Add the latest change from the master

.github/workflows/python-app.yml

@@ -24,6 +24,8 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install flake8 pytest
+        pip install ase spglib==2.2 julia
+        python -c 'import julia; julia.install()'
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
@@ -37,15 +39,16 @@ jobs:
         sudo apt-get install gfortran libblas-dev liblapack-dev build-essential \
                              mpich libopenmpi-dev libgl1-mesa-glx  libglib2.0-0 \
                              libboost-all-dev
+
+        pip install mpi4py
 
-        pip install mpi4py numpy matplotlib scipy spglib --user
 
         git clone https://github.com/SSCHAcode/CellConstructor.git
-        cd cellconstructor && sudo python setup.py install && cd ..
-        rm -r cellconstructor
+        cd CellConstructor && python setup.py install --user && cd ..
+        rm -r CellConstructor
 
         git clone https://github.com/SSCHAcode/python-sscha.git
-        cd python-sscha && sudo python setup.py install && cd ..
+        cd python-sscha && python setup.py install --user && cd ..
         rm -r python-sscha
 
         MPICC=mpicc python setup.py install --user

CModules/Lanczos.cpp

@@ -58,6 +58,18 @@ void Lanczos::setup_from_input(string rootname) {
     ignore_v3 = root.get<bool>("ignore_v3");
     ignore_v4 = root.get<bool>("ignore_v4");
 
+    // Check the optional key restart (default false)
+    int step_done = 0;
+    restart = false;
+    if (root.count("restart")) {
+        restart = root.get<bool>("restart");
+        step_done = root.get<int>("step_done");
+    }
+
+    if (step_done >= n_steps) {
+        throw invalid_argument("Error, n_steps must be greather than step_done.\n");
+    }
+
     // Now read the more complex values
     pt::ptree &data = root.get_child("data");
     N = data.get<int>("n_configs");
@@ -143,11 +155,12 @@ void Lanczos::setup_from_input(string rootname) {
     }
     file.close();
 
-    // Now read the 2D arrays
+
+    // Now read the 2D arrays (TODO: read also previous steps)
     file.open(rootname + ".X.dat");
     if (file.is_open()) {
         string line;
-        for (int k = 0; k < N*n_modes; ++k) {
+        for (int k = 0; k < N*n_modes ; ++k) {
             file >> X[k];
         }
     }
@@ -157,13 +170,58 @@ void Lanczos::setup_from_input(string rootname) {
     file.open(rootname + ".Y.dat");
     if (file.is_open()) {
         string line;
-        for (int k = 0; k < N*n_modes; ++k) {
+        for (int k = 0; k < N*n_modes ; ++k) {
             file >> Y[k];
         }
     }
     file.close();
 
 
+    // Check if we need a restart
+    // And load P and Q
+    if (restart) {
+        file.open(rootname + ".abc");
+        if (file.is_open()) {
+            string line;
+            for (int k = 0; k < step_done ; ++k) {
+                file >> a[k];
+                file >> b[k];
+                file >> c[k];
+            }
+        }
+        file.close();
+
+        file.open(rootname + ".qbasis.out");
+        if (file.is_open()) {
+            string line;
+            for (int k = 0; k < n_psi * (step_done+1) ; ++k) {
+                file >> Qbasis[k];
+            }
+        }
+        file.close();
+
+        file.open(rootname + ".pbasis.out");
+        if (file.is_open()) {
+            string line;
+            for (int k = 0; k < n_psi * (step_done+1) ; ++k) {
+                file >> Pbasis[k];
+            }
+        }
+        file.close();
+
+        file.open(rootname + ".snorm.out");
+        if (file.is_open()) {
+            string line;
+            for (int k = 0; k < step_done + 1; ++k) {
+                file >> snorm[k];
+            }
+        }
+        file.close();
+
+        i_step = step_done;
+    }
+
+
     file.open(rootname + ".blockid");
     if (file.is_open()) {
         string line;
@@ -480,8 +538,8 @@ void Lanczos::run() {
     double * p_L = (double*) malloc(sizeof(double) * n_psi);
 
     // Store the lenght of the variables
-    int lena = 0, lenb = 0, lenc = 0, lens = 1;
-    int len_bq = 1, len_bp = 1;
+    int lena = i_step, lenb = i_step, lenc = i_step, lens = i_step + 1;
+    int len_bq = i_step + 1, len_bp = i_step + 1;
 
     // Prepare the first vector computing the norm of psi
     double psi_norm = 0;
@@ -490,26 +548,23 @@ void Lanczos::run() {
 
 
     // Fill the first vector of P and Q basis
+
+    psi_q = (double *) malloc(sizeof(double) * n_psi);
+    psi_p = (double *) malloc(sizeof(double) * n_psi);
+    double * sk_tilde = (double*) malloc(sizeof(double) * n_psi);
+
     if (i_step == 0) {
         for (int i = 0; i < n_psi; ++i) {
             Qbasis[i] = psi[i] / psi_norm;
             Pbasis[i] = psi[i] / psi_norm;
         }
         snorm[0] = 1;
-    } else {
-        cerr << "Error, starting from a non zero step is not implemented." << endl;
-        exit(EXIT_FAILURE);
-    }
-
-
-    psi_q = (double *) malloc(sizeof(double) * n_psi);
-    psi_p = (double *) malloc(sizeof(double) * n_psi);
-    double * sk_tilde = (double*) malloc(sizeof(double) * n_psi);
+    } 
 
 
     for (int j = 0; j < n_psi; ++j) {
-        psi_q[j] = Qbasis[j];
-        psi_p[j] = Pbasis[j];
+        psi_q[j] = Qbasis[i_step * n_psi + j];
+        psi_p[j] = Pbasis[i_step * n_psi + j];
     }
 
     double a_coeff, b_coeff, c_coeff;
@@ -521,30 +576,49 @@ void Lanczos::run() {
     fstream file_snorm;
 
     if (am_i_the_master()) {
-        file_abc.open(rootname + ".abc", fstream::out);
-        file_qbasis.open(rootname + ".qbasis.out", fstream::out);
-        file_pbasis.open(rootname + ".pbasis.out", fstream::out);
-        file_snorm.open(rootname + ".snorm.out", fstream::out);
+        if (restart) {
+            file_abc.open(rootname + ".abc", fstream::out | fstream::app);
+            file_qbasis.open(rootname + ".qbasis.out", fstream::out | fstream::app);
+            file_pbasis.open(rootname + ".pbasis.out", fstream::out | fstream::app);
+            file_snorm.open(rootname + ".snorm.out", fstream::out | fstream::app);
+        } else {
+            file_abc.open(rootname + ".abc", fstream::out);
+            file_qbasis.open(rootname + ".qbasis.out", fstream::out);
+            file_pbasis.open(rootname + ".pbasis.out", fstream::out);
+            file_snorm.open(rootname + ".snorm.out", fstream::out);
+        }
     }
 
-    if (file_qbasis.is_open()) {
+    // Setup the required precision after the restart
+    if (file_qbasis.is_open()) 
         file_qbasis << scientific << setprecision(16);
-        for (int j = 0; j < n_psi; ++j) {
-            file_qbasis << psi_q[j] << " ";
-        }
-        file_qbasis << endl;
 
-    } 
-    if (file_pbasis.is_open()) { 
+    if (file_pbasis.is_open())
         file_pbasis << scientific << setprecision(16);
-        for (int j = 0; j < n_psi; ++j) {
-            file_pbasis << psi_p[j] << " ";
+
+    if (file_snorm.is_open())
+        file_snorm << scientific << setprecision(16);
+
+    if (! restart) {
+        if (file_qbasis.is_open()) {
+            file_qbasis << scientific << setprecision(16);
+            for (int j = 0; j < n_psi; ++j) {
+                file_qbasis << psi_q[j] << " ";
+            }
+            file_qbasis << endl;
+
+        } 
+        if (file_pbasis.is_open()) { 
+            file_pbasis << scientific << setprecision(16);
+            for (int j = 0; j < n_psi; ++j) {
+                file_pbasis << psi_p[j] << " ";
+            }
+            file_pbasis << endl;
         }
-        file_pbasis << endl;
-    }
 
-    if (file_snorm.is_open()) {
-        file_snorm << scientific << setprecision(16) << snorm[0] << endl << flush;
+        if (file_snorm.is_open()) {
+            file_snorm << scientific << setprecision(16) << snorm[0] << endl << flush;
+        }
     }
 
     // Here the run
@@ -555,7 +629,7 @@ void Lanczos::run() {
     if (DEBUG_LANC && am_i_the_master()) 
         cout << "PSI POINTER (LINE " << __LINE__ << "): " << psi << endl;
 
-    for (int i = i_step; i < i_step + n_steps; ++i) {
+    for (int i = i_step; i < n_steps; ++i) {
         if (am_i_the_master()) {
             cout << endl;
             cout << "===== NEW STEP " << i + 1 <<" =====" <<endl << endl;

CModules/Lanczos.hpp

@@ -16,7 +16,7 @@ using namespace std;
 class Lanczos {
     int N, n_syms, n_modes, n_steps, n_blocks;
     int i_step;
-    bool ignore_v2, ignore_v3, ignore_v4, reverse_L;
+    bool ignore_v2, ignore_v3, ignore_v4, reverse_L, restart;
 
     double T, shift_value;
 

CModules/LanczosFunctions.c

@@ -3,6 +3,8 @@
 
 #define DEB 0
 #define DEB_L 0
+#define USEBLAS 
+
 
 // These are used for debugging
 #define X_VAL 1534
@@ -2191,6 +2193,7 @@ void get_f_average_from_Y_pert_sym_fast( double * X,  double * Y,  double * w,
 
 			// Compute the standard weight Y1
 			weight = 0;
+			// TODO, explot BLAS
 			for (mu = 0; mu < n_modes; ++mu) {
 				for (nu = 0; nu < n_modes; ++nu) {
 					weight -= displacement[mu] * displacement[nu] * Y1[mu * n_modes +  nu] / 2;
@@ -2261,6 +2264,7 @@ void get_d2v_dR2_from_R_pert_sym_fast( double * X,  double * Y,  double * w,  do
 	double N_eff = 0;
 	double u_mu, u_nu;
 
+	double sum = 0;
 
 	// Prepare the temporaney force and displacement (after symmetry applicaiton)
 	double * force = (double*) calloc(sizeof(double), n_modes);
@@ -2374,6 +2378,8 @@ void get_d2v_dR2_from_R_pert_sym_fast( double * X,  double * Y,  double * w,  do
 				weight += f_ups(w[mu], T) * R1[mu] * displacement[mu];
 			}
 
+			//printf("i = %d, j = %d, weight = %e\n", i, j, weight);
+
 			// Compute the d2V_dR2
 			for (mu = 0; mu < n_modes; ++mu) {
 				u_mu = f_ups(w[mu], T) * displacement[mu];
@@ -2399,6 +2405,10 @@ void get_d2v_dR2_from_R_pert_sym_fast( double * X,  double * Y,  double * w,  do
 					d2v_tmp[mu * n_modes + nu] -= u_mu * u_nu * weight * w_is[i] / 3;
 				}
 			}	
+
+			/*sum = 0;
+			for (mu = 0; mu < n_modes * n_modes; ++mu) sum += d2v_tmp[mu];
+			printf("Total sum: %e\n", sum);*/
 		}
 	}
 
@@ -2678,16 +2688,17 @@ void get_d2v_dR2_from_Y_pert_sym_fast( double * X,  double * Y,  double * w,  do
 			// First the standard weight
 			for (mu = 0; mu < n_modes; ++mu) {
 				for(nu = 0; nu < n_modes; ++nu) {
-					weight -= displacement[mu] * displacement[nu] * Y1[mu * n_modes + nu] / 2;
+					weight -= displacement[mu] * displacement[nu] * Y1[mu * n_modes + nu];
 				}
 			}
+			weight *= w_is[i] / 8;
 
 			// Now the first part of the potential
 			for (mu = 0; mu < n_modes; ++mu) {
 				u_mu = f_ups(w[mu], T) * displacement[mu];
 				for(nu = 0; nu < n_modes; ++nu) {
-					d2v_tmp[mu * n_modes + nu] -= force[nu] * u_mu * weight * w_is[i] / 4; // Permutation symmetry
-					d2v_tmp[nu * n_modes + mu] -= force[nu] * u_mu * weight * w_is[i] / 4; // Permutation symmetry
+					d2v_tmp[mu * n_modes + nu] -= force[nu] * u_mu * weight; // Permutation symmetry
+					d2v_tmp[nu * n_modes + mu] -= force[nu] * u_mu * weight; // Permutation symmetry
 				}
 			}
 
@@ -2698,17 +2709,18 @@ void get_d2v_dR2_from_Y_pert_sym_fast( double * X,  double * Y,  double * w,  do
 				f_mu = f_psi(w[mu], T) * force[mu];
 				for(nu = 0; nu < n_modes; ++nu) {
 					f_nu = f_psi(w[nu], T) * force[nu];
-					weight -= f_mu * displacement[nu] * Y1[mu * n_modes + nu] / 4;
-					weight -= f_nu * displacement[mu] * Y1[mu * n_modes + nu] / 4;
+					weight -= f_mu * displacement[nu] * Y1[mu * n_modes + nu];
+					weight -= f_nu * displacement[mu] * Y1[mu * n_modes + nu];
 				}
 			}
+			weight *= w_is[i] / 8;
 
 			// Now the first part of the potential
 			for (mu = 0; mu < n_modes; ++mu) {
 				u_mu = f_ups(w[mu], T) * displacement[mu];
 				for(nu = 0; nu < n_modes; ++nu) {
 					u_nu = f_ups(w[nu], T) * displacement[nu];
-					d2v_tmp[mu * n_modes + nu] -= u_nu * u_mu * weight * w_is[i] / 2; // Permutation symmetry
+					d2v_tmp[mu * n_modes + nu] -= u_nu * u_mu * weight; // Permutation symmetry
 				}
 			}
 		}

CModules/Utils.cpp

@@ -1,10 +1,18 @@
 #include "Utils.hpp"
 
+using namespace std;
 bool am_i_the_master(void) {
     int size, rank;
 	MPI_Comm_size(MPI_COMM_WORLD, &size);
 	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
     if (rank == 0) return true;
     return false;
+}
+
+
+bool is_file_exist(const char *fileName)
+{
+    std::ifstream infile(fileName);
+    return infile.good();
 }

CModules/Utils.hpp

@@ -1,6 +1,9 @@
 #ifndef LANC_UTILS
 #define LANC_UTILS
 #include <mpi.h>
+#include <iostream>
+#include <fstream>
 
 bool am_i_the_master(void);
+bool is_file_exist(const char *fileName);
 #endif