From b57b512e77d0d102a372384acfe3a0b73509f458 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 16 Jul 2019 16:18:52 +0200 Subject: [PATCH 01/41] Added the q vector for each polarization mode in real space supercell in the Lanczos. --- Modules/DynamicalLanczos.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index e9d643d51..0e7ba1748 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -115,6 +115,7 @@ def __init__(self, ensemble = None, mode = 1): self.N_degeneracy = None self.initialized = False self.perturbation_modulus = 1 + self.q_vectors = None # The q vectors of each mode # Perform a bare initialization if the ensemble is not provided if ensemble is None: @@ -130,6 +131,7 @@ def __init__(self, ensemble = None, mode = 1): ws, pols = self.dyn.DiagonalizeSupercell() + self.nat = superdyn.structure.N_atoms self.qe_sym = CC.symmetries.QE_Symmetry(self.dyn.structure) @@ -148,6 +150,17 @@ def __init__(self, ensemble = None, mode = 1): self.n_modes = len(self.w) + + # Prepare the list of q point starting from the polarization vectors + q_list = CC.symmetries.GetQForEachMode(self.pols, self.uci_structure, self.super_structure, self.dyn.GetSupercell()) + # Store the q vectors in crystal space + bg = self.uci_structure.get_reciprocal_vectors() / 2* np.pi + self.q_vectors = np.zeros((self.n_modes, 3), dtype = np.double, order = "C") + for iq, q in enumerate(q_list): + self.q_vectors[iq, :] = CC.Methods.covariant_coordinate(bg, q) + + + # Ignore v3 or v4. You can set them for testing self.ignore_v3 = False self.ignore_v4 = False @@ -645,7 +658,8 @@ def save_status(self, file): N_degeneracy = self.N_degeneracy, initialized = self.initialized, degenerate_space = self.degenerate_space, - perturbation_modulus = self.perturbation_modulus) + perturbation_modulus = self.perturbation_modulus, + q_vectors = self.q_vectors) def load_status(self, file): """ @@ -691,6 +705,7 @@ def load_status(self, file): if "perturbation_modulus" in data.keys(): self.perturbation_modulus = data["perturbation_modulus"] + self.q_vectors = data["q_vectors"] # Rebuild the Linear operator self.L_linop = scipy.sparse.linalg.LinearOperator(shape = (len(self.psi), len(self.psi)), matvec = self.apply_full_L, dtype = TYPE_DP) From 6057961c1dc8ea739dfbb4b44d54823829e9962f Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 12:06:26 +0200 Subject: [PATCH 02/41] Added the way of doing many Lanczos perturbations to get the curvature --- Modules/DynamicalLanczos.py | 259 +++++++++++++++++++++++++++++++++++- 1 file changed, 256 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 0e7ba1748..133fbabfc 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -203,6 +203,28 @@ def __init__(self, ensemble = None, mode = 1): self.shift_value = 0 self.symmetrize = False + def reset(self): + """ + RESET THE LANCZOS + ================= + + This function reset the Lanczos algorithm, allowing for a new responce function calculation + with the same ensemble and the same settings. + """ + + + # Prepare the solution of the Lanczos algorithm + self.eigvals = None + self.eigvects = None + + # Store the basis and the coefficients of the Lanczos procedure + # In the custom lanczos mode + self.a_coeffs = [] #Coefficients on the diagonal + self.b_coeffs = [] # Coefficients close to the diagonal + self.krilov_basis = [] # The basis of the krilov subspace + self.arnoldi_matrix = [] # If requested, the upper triangular arnoldi matrix + + def prepare_symmetrization(self): @@ -263,11 +285,11 @@ def prepare_symmetrization(self): def prepare_ir(self, effective_charges = None, pol_vec = np.array([1,0,0])): """ - PREPARE LANCZOS FOR IR COMPUTATION - ================================== + PREPARE LANCZOS FOR INFRARED SPECTRUM COMPUTATION + ================================================= In this subroutine we prepare the lanczos algorithm for the computation of the - IR signal. + signal. Parameters ---------- @@ -1170,6 +1192,147 @@ def get_static_odd_fc(self, use_arnoldi = False): return fc_matrix + def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_dir = None, verbose = True): + """ + GET ALL THE GREEN FUNCTIONS + =========================== + + This will compute a set of lanczos coefficients for each element of the odd matrix. + a_n and b_n. + We will run lanczos for all the elements and all the crosses. + In this way we have the whole evolution with frequency of the matrix. + + NOTE: This can be a very intensive computation. + + Parameters + ---------- + N_steps : int + The number of Lanczos iteration for each green function + mode_mixing : bool + If True also non diagonal elements are computed, otherwise the + SSCHA eigenvector are supposed to be conserved, and only diagonal + green functions are considered. + If False the computation is much less expensive (a factor nat_sc), + but it is approximated. + save_step_dir : string + If not None, the path to the directory in which you want to save + each step. So even if stopped the calculation can restart. + verbose : bool + If true print all the progress to standard output + + Results + ------- + a_ns : ndarray( (n_modes, n_modes, N_steps)) + The a coefficients for each element in the mode x mode space + b_ns : ndarray( (n_modes, n_modes, N_steps-1)) + The b_n coefficients for each mode in the space. + """ + + # Check if the save directory exists + # Otherwise we create it + if not save_step_dir is None: + if not os.path.exists(save_step_dir): + os.makedirs(save_step_dir) + + # Load all the data + a_ns = np.zeros( (self.n_modes, self.n_modes, N_steps), dtype = np.double) + b_ns = np.zeros( (self.n_modes, self.n_modes, N_steps), dtype = np.double) + + # Incompatible with shift for now + self.shift_value = 0 + + # Compute the diagonal parts + for i in range(self.n_modes): + if verbose: + print("\n") + print(" ========================== ") + print(" | | ") + print(" | DIAGONAL ELEMENTS | ") + print(" | STEP {:5d} | ".format(i)) + print(" | | ") + print(" ========================== ") + print() + + # Setup the Lanczos + self.reset() + + # Prepare the perturbation + self.psi[:] = 0 + self.psi[i] = 1 + + # Run the Lanczos perturbation + self.run(N_steps, save_dir = save_step_dir, verbose = verbose) + + if verbose: + print() + print(" ---- > LANCZOS RUN COMPLEATED < ---- ") + print() + + # Save the status + if save_step_dir: + self.save_status("full_lanczos_diagonal_{}".format(i)) + + # Fill the a_n and b_n + a_ns[i, i, :] = self.a_coeffs + b_ns[i, i, :] = self.b_coeffs + + # If we must compute the mode mixing + if mode_mixing: + for i in range(self.n_modes): + for j in range(i+1, self.n_modes): + # TODO: Neglect (i,j) forbidden by symmetries + + if verbose: + print("\n") + print(" ============================ ") + print(" | | ") + print(" | NON DIAGONAL ELEMENT | ") + print(" | STEP ({:5d},{:5d}) | ".format(i, j)) + print(" | | ") + print(" ============================ ") + print() + + # Setup the Lanczos + self.reset() + + # Prepare the perturbation + self.psi[:] = 0 + self.psi[i] = 1 + self.psi[j] = 1 + + # Run the Lanczos perturbation + self.run(N_steps, save_dir = save_step_dir, verbose = verbose) + + if verbose: + print() + print(" ---- > LANCZOS RUN COMPLEATED < ---- ") + print() + + # Save the status + if save_step_dir: + self.save_status("full_lanczos_off_diagonal_{}_{}".format(i, j)) + + # Fill the a_n and b_n + a_ns[i, j, :] = self.a_coeffs + b_ns[i, j, :] = self.b_coeffs + a_ns[j, i, :] = self.a_coeffs + b_ns[j, i, :] = self.b_coeffs + + + if verbose: + print() + print() + print(" ================= ") + print(" | | ") + print(" | DONE | ") + print(" | | ") + print(" ================= ") + print() + print() + + return a_ns, b_ns + + def get_spectral_function_from_Lenmann(self, w_array, smearing, use_arnoldi=True): """ @@ -1778,3 +1941,93 @@ def FastApplyD4ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, deg + + +# Here some functions to analyze the data that comes out by a Lanczos +def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mode_mixing = False,\ + use_terminator = True, last_average = 5, smearing = 0): + """ + GET THE FREE ENERGY CURVATURE FROM MANY LANCZOS + =============================================== + + This function computes the free energy curvature from the result + of a full Lanczos computation between all possible perturbations. + + Parameters + ---------- + a_ns : ndarray(size = (n_modes, n_modes, N_steps)) + The a_n coefficients for each Lanczos perturbation + b_ns : ndarray(size = (n_modes, n_modes, N_steps-1)) + The b_n coefficients for each Lanczos perturbation + pols_sc : ndarray(size = (3*nat_sc, n_modes)) + The polarization vectors in the supercell + masses : ndarray(size = (3*nat_sc)) + The mass associated to each component of pols_sc + use_terminator : bool + If true the infinite volume interpolation is performed trought the + terminator trick + last_average : int + Used in combination with the terminator, average the last 'last_average' + coefficients and replicate them. + smearing : float + The smearing for the green function calculation. + Usually not needed for this kind of calculation. + + Results + ------- + odd_fc : ndarray( (3*nat_sc, 3*nat_sc)) + The free energy curvature in the supercell + + """ + + n_modes = np.shape(pols_sc)[1] + nat_sc = int(np.shape(pols_sc)[0] / 3) + N_steps = np.shape(a_ns)[2] + + assert N_steps == np.shape(b_ns)[2], "Error, an and bn has an incompatible size:\n a_n = {}, b_n = {}".format(np.shape(a_ns), np.shape(b_ns)) + + + mat_pol = np.zeros( (n_modes, n_modes), dtype = np.double) + for i in range(n_modes): + + # Create the Lanczos class + lanc = Lanczos(None) + lanc.a_coeffs = a_ns[i, i, :] + lanc.b_coeffs = b_ns[i, i, :] + lanc.perturbation_modulus = 1 + + # get the green function from continued fraction + gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ + smearing = smearing, last_average = last_average)[0] + + mat_pol[i,i] = np.real(gf) + + # If there is the mode-mixing compute also the off-diagonal terms + if mode_mixing: + for i in range(n_modes): + for j in range(i+1, n_modes): + + # Create the Lanczos class + lanc = Lanczos(None) + lanc.a_coeffs = a_ns[i, j, :] + lanc.b_coeffs = b_ns[i, j, :] + lanc.perturbation_modulus = 2 + + # get the green function from continued fraction + gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ + smearing = smearing, last_average = last_average)[0] + + # Lanczos can compute only diagonal green functions + # Therefore we need to trick it to get the off-diagonal elements + # <1|L|2> = 1/2*( <1+2|L|1+2> - <1|L|1> - <2|L|2>) + mat_pol[i,j] = (np.real(gf) - mat_pol[i,i] - mat_pol[j,j]) / 2 + mat_pol[j,i] = (np.real(gf) - mat_pol[i,i] - mat_pol[j,j]) / 2 + + # The green function is the inverse of the free energy curvature + fc_pols = np.linalg.inv(mat_pol) + + # Get back into real space + epols_m = np.einsum("ab, a->ab", pols_sc, np.sqrt(masses)) + fc_odd = np.einsum("ab, ca, da ->cd", fc_pols, epols_m, epols_m) + + return fc_odd From 646bc903555005333b254592fa952ae5f186e64f Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 12:23:56 +0200 Subject: [PATCH 03/41] Added the ability to use Lanczos to compute static free energy curvature --- Modules/DynamicalLanczos.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 133fbabfc..20f2843ba 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -156,8 +156,8 @@ def __init__(self, ensemble = None, mode = 1): # Store the q vectors in crystal space bg = self.uci_structure.get_reciprocal_vectors() / 2* np.pi self.q_vectors = np.zeros((self.n_modes, 3), dtype = np.double, order = "C") - for iq, q in enumerate(q_list): - self.q_vectors[iq, :] = CC.Methods.covariant_coordinate(bg, q) + #for iq, q in enumerate(q_list): + # self.q_vectors[iq, :] = CC.Methods.covariant_coordinate(bg, q) From 76e87995e1cf89c780a8b6133bbc5c20be92358a Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 12:25:40 +0200 Subject: [PATCH 04/41] bugfix --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 20f2843ba..18a480cb0 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -152,7 +152,7 @@ def __init__(self, ensemble = None, mode = 1): # Prepare the list of q point starting from the polarization vectors - q_list = CC.symmetries.GetQForEachMode(self.pols, self.uci_structure, self.super_structure, self.dyn.GetSupercell()) + #q_list = CC.symmetries.GetQForEachMode(self.pols, self.uci_structure, self.super_structure, self.dyn.GetSupercell()) # Store the q vectors in crystal space bg = self.uci_structure.get_reciprocal_vectors() / 2* np.pi self.q_vectors = np.zeros((self.n_modes, 3), dtype = np.double, order = "C") From 358c3e32c664344e43e9284f66e564ee1274de28 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 16:54:31 +0200 Subject: [PATCH 05/41] Bugfix on lanczos --- Modules/DynamicalLanczos.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 18a480cb0..f47804385 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1273,8 +1273,12 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d self.save_status("full_lanczos_diagonal_{}".format(i)) # Fill the a_n and b_n - a_ns[i, i, :] = self.a_coeffs - b_ns[i, i, :] = self.b_coeffs + a_tmp = np.zeros(N_steps, dtype = np.double) + a_tmp[len(self.a_coeffs)] = self.a_coeffs + b_tmp = np.zeros(N_steps-1, dtype = np.double) + b_tmp[len(self.b_coeffs)] = self.b_coeffs + a_ns[i, i, :] = a_tmp + b_ns[i, i, :] = b_tmp # If we must compute the mode mixing if mode_mixing: From 7040d06fa3dea18053f6836bcb031075d6d8fd56 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 16:55:40 +0200 Subject: [PATCH 06/41] Other bugfix on lanczos --- Modules/DynamicalLanczos.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index f47804385..40e157441 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1274,9 +1274,9 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d # Fill the a_n and b_n a_tmp = np.zeros(N_steps, dtype = np.double) - a_tmp[len(self.a_coeffs)] = self.a_coeffs + a_tmp[:len(self.a_coeffs)] = self.a_coeffs b_tmp = np.zeros(N_steps-1, dtype = np.double) - b_tmp[len(self.b_coeffs)] = self.b_coeffs + b_tmp[:len(self.b_coeffs)] = self.b_coeffs a_ns[i, i, :] = a_tmp b_ns[i, i, :] = b_tmp From 57d698a51dbad476ea2ed02f840a73469e697367 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:00:29 +0200 Subject: [PATCH 07/41] Added a fix if lanczos converges easily --- Modules/DynamicalLanczos.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 40e157441..955caace8 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1236,7 +1236,7 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d # Load all the data a_ns = np.zeros( (self.n_modes, self.n_modes, N_steps), dtype = np.double) - b_ns = np.zeros( (self.n_modes, self.n_modes, N_steps), dtype = np.double) + b_ns = np.zeros( (self.n_modes, self.n_modes, N_steps-1), dtype = np.double) # Incompatible with shift for now self.shift_value = 0 @@ -2010,11 +2010,13 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod if mode_mixing: for i in range(n_modes): for j in range(i+1, n_modes): + # Get the number of steps + n_steps = np.arange(N_steps)[b_ns[i, j, :] == 0][0] + 1 - # Create the Lanczos class + # Create the Lanczos class) lanc = Lanczos(None) - lanc.a_coeffs = a_ns[i, j, :] - lanc.b_coeffs = b_ns[i, j, :] + lanc.a_coeffs = a_ns[i, j, :n_steps] + lanc.b_coeffs = b_ns[i, j, :n_steps] lanc.perturbation_modulus = 2 # get the green function from continued fraction From 707515c1e2978fb0e0ab2a5ef793656b4db3d63c Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:02:47 +0200 Subject: [PATCH 08/41] Another bugfix on lanczos --- Modules/DynamicalLanczos.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 955caace8..72b8b4f40 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1317,10 +1317,14 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d self.save_status("full_lanczos_off_diagonal_{}_{}".format(i, j)) # Fill the a_n and b_n - a_ns[i, j, :] = self.a_coeffs - b_ns[i, j, :] = self.b_coeffs - a_ns[j, i, :] = self.a_coeffs - b_ns[j, i, :] = self.b_coeffs + a_tmp = np.zeros(N_steps, dtype = np.double) + a_tmp[:len(self.a_coeffs)] = self.a_coeffs + b_tmp = np.zeros(N_steps-1, dtype = np.double) + b_tmp[:len(self.b_coeffs)] = self.b_coeffs + a_ns[i, j, :] = a_tmp + b_ns[i, j, :] = b_tmp + a_ns[j, i, :] = a_tmp + b_ns[j, i, :] = b_tmp if verbose: From cdb5b973bdfb8c7b02f22118140f21a52d9b6d1d Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:04:45 +0200 Subject: [PATCH 09/41] Another bugfix --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 72b8b4f40..9b5f45c8b 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1992,7 +1992,7 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod nat_sc = int(np.shape(pols_sc)[0] / 3) N_steps = np.shape(a_ns)[2] - assert N_steps == np.shape(b_ns)[2], "Error, an and bn has an incompatible size:\n a_n = {}, b_n = {}".format(np.shape(a_ns), np.shape(b_ns)) + assert N_steps -1 == np.shape(b_ns)[2], "Error, an and bn has an incompatible size:\n a_n = {}, b_n = {}".format(np.shape(a_ns), np.shape(b_ns)) mat_pol = np.zeros( (n_modes, n_modes), dtype = np.double) From 26053b5229a89132554ae9602cedc21a254255d7 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:07:35 +0200 Subject: [PATCH 10/41] Aded debug --- Modules/DynamicalLanczos.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 9b5f45c8b..20e77e56e 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2020,9 +2020,11 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod # Create the Lanczos class) lanc = Lanczos(None) lanc.a_coeffs = a_ns[i, j, :n_steps] - lanc.b_coeffs = b_ns[i, j, :n_steps] + lanc.b_coeffs = b_ns[i, j, :n_steps-1] lanc.perturbation_modulus = 2 + print("Computing ({},{}) ... n_steps = {}".format(i, j, n_steps)) + # get the green function from continued fraction gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ smearing = smearing, last_average = last_average)[0] From 325d94b8a1d9bd2116e0b83c69da17659dc14ab7 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:09:25 +0200 Subject: [PATCH 11/41] Other bugfix --- Modules/DynamicalLanczos.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 20e77e56e..b55fb6fdf 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1997,13 +1997,20 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod mat_pol = np.zeros( (n_modes, n_modes), dtype = np.double) for i in range(n_modes): + + # Get the number of steps + n_steps = np.arange(N_steps)[b_ns[i, i, :] == 0][0] + 1 + # Create the Lanczos class lanc = Lanczos(None) - lanc.a_coeffs = a_ns[i, i, :] - lanc.b_coeffs = b_ns[i, i, :] + lanc.a_coeffs = a_ns[i, i, :n_steps] + lanc.b_coeffs = b_ns[i, i, :n_steps - 1] lanc.perturbation_modulus = 1 + + print("Computing ({},{}) ... n_steps = {}".format(i, j, n_steps)) + # get the green function from continued fraction gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ smearing = smearing, last_average = last_average)[0] From 26508a2468350228fafeab477f215f0e2cd0c197 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:09:40 +0200 Subject: [PATCH 12/41] Changed default mode_mixing --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index b55fb6fdf..809d45389 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1952,7 +1952,7 @@ def FastApplyD4ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, deg # Here some functions to analyze the data that comes out by a Lanczos -def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mode_mixing = False,\ +def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mode_mixing = True,\ use_terminator = True, last_average = 5, smearing = 0): """ GET THE FREE ENERGY CURVATURE FROM MANY LANCZOS From 44676cbdaa3ae3854b16ead09a2647ca59618a33 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:12:08 +0200 Subject: [PATCH 13/41] Bugfix on the lanczos --- Modules/DynamicalLanczos.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 809d45389..3d1a430b1 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1999,7 +1999,11 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod for i in range(n_modes): # Get the number of steps - n_steps = np.arange(N_steps)[b_ns[i, i, :] == 0][0] + 1 + n_steps = np.arange(N_steps-1)[b_ns[i, i, :] == 0] + if len(n_steps) == 0: + n_steps = N_steps + else: + n_steps = n_steps[0] + 1 # Create the Lanczos class @@ -2022,7 +2026,12 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod for i in range(n_modes): for j in range(i+1, n_modes): # Get the number of steps - n_steps = np.arange(N_steps)[b_ns[i, j, :] == 0][0] + 1 + n_steps = np.arange(N_steps-1)[b_ns[i, i, :] == 0] + if len(n_steps) == 0: + n_steps = N_steps + else: + n_steps = n_steps[0] + 1 + # Create the Lanczos class) lanc = Lanczos(None) From 5f90de99cba1d0d008fcea66506e7ed817a35e3a Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:13:07 +0200 Subject: [PATCH 14/41] Stupid bug corrected --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 3d1a430b1..26d774054 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2013,7 +2013,7 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod lanc.perturbation_modulus = 1 - print("Computing ({},{}) ... n_steps = {}".format(i, j, n_steps)) + print("Computing ({},{}) ... n_steps = {}".format(i, i, n_steps)) # get the green function from continued fraction gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ From c77192f9eceb5e59ab83c22d63ea9fc5bf2e9343 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:14:55 +0200 Subject: [PATCH 15/41] Added some extra debug --- Modules/DynamicalLanczos.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 26d774054..d136cf858 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2052,7 +2052,9 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod mat_pol[j,i] = (np.real(gf) - mat_pol[i,i] - mat_pol[j,j]) / 2 # The green function is the inverse of the free energy curvature + np.savetxt("gf_mat.dat", mat_pol) fc_pols = np.linalg.inv(mat_pol) + np.savetxt("fc_pols.dat", fc_pols) # Get back into real space epols_m = np.einsum("ab, a->ab", pols_sc, np.sqrt(masses)) From 2f6ea964b9c721c226a9c1f149a319116252581c Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:27:37 +0200 Subject: [PATCH 16/41] Added a bugfix on lanczos --- Modules/DynamicalLanczos.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index d136cf858..af82b8b18 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -975,6 +975,7 @@ def run(self, n_iter, save_dir = ".", verbose = True): # Lets repeat twice the orthogonalization + converged = False for k_orth in range(N_REP_ORTH): for j in range(len(self.krilov_basis)): coeff = new_vect.dot(self.krilov_basis[j]) @@ -992,15 +993,17 @@ def run(self, n_iter, save_dir = ".", verbose = True): # Check the normalization (If zero the algorithm converged) if norm < __EPSILON__: + converged = True if verbose: print("Obtained a linear dependent vector.") print("The algorithm converged.") - return + break new_vect /= norm - self.krilov_basis.append(new_vect) - self.psi = new_vect + if not converged: + self.krilov_basis.append(new_vect) + self.psi = new_vect t2 = time.time() # Add the coefficients to the variables @@ -1026,6 +1029,10 @@ def run(self, n_iter, save_dir = ".", verbose = True): if verbose: print("Lanczos step %d ultimated." % i) + + + if converged: + return def build_lanczos_matrix_from_coeffs(self, use_arnoldi=True): From 32bf4fe15c8c43a8be8f2532efa10b2210cd2336 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 17 Jul 2019 17:48:44 +0200 Subject: [PATCH 17/41] Bugfix on lanczos odd --- Modules/DynamicalLanczos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index af82b8b18..bc9b987dc 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1030,7 +1030,7 @@ def run(self, n_iter, save_dir = ".", verbose = True): if verbose: print("Lanczos step %d ultimated." % i) - + if converged: return @@ -2033,7 +2033,7 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod for i in range(n_modes): for j in range(i+1, n_modes): # Get the number of steps - n_steps = np.arange(N_steps-1)[b_ns[i, i, :] == 0] + n_steps = np.arange(N_steps-1)[b_ns[i, j, :] == 0] if len(n_steps) == 0: n_steps = N_steps else: @@ -2046,7 +2046,7 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod lanc.b_coeffs = b_ns[i, j, :n_steps-1] lanc.perturbation_modulus = 2 - print("Computing ({},{}) ... n_steps = {}".format(i, j, n_steps)) + print("Computing ({},{}) ..., n_steps = {}".format(i, j, n_steps)) # get the green function from continued fraction gf = lanc.get_green_function_continued_fraction(np.array([0]), use_terminator = use_terminator, \ From bac5651497d6de07bc1ab8380108367dd86a5278 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 10:47:47 +0200 Subject: [PATCH 18/41] Added a total timing function in the super Lanczos --- Modules/DynamicalLanczos.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index bc9b987dc..6d9c7f692 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1235,6 +1235,9 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d The b_n coefficients for each mode in the space. """ + # Time the function + t_start = time.time() + # Check if the save directory exists # Otherwise we create it if not save_step_dir is None: @@ -1333,15 +1336,22 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d a_ns[j, i, :] = a_tmp b_ns[j, i, :] = b_tmp + t_end = time.time() + + total_time = t_end - t_start + minutes = int(total_time / 60) + hours = int(minutes / 60) + minutes -= hours * 60 + seconds = total_time - hours*3600 - minutes * 60 if verbose: print() print() - print(" ================= ") - print(" | | ") - print(" | DONE | ") - print(" | | ") - print(" ================= ") + print(" ====================== ") + print(" | | ") + print(" | DONE | ") + print(" | In {:3d}:{:02d}:{:02d}s | ".format(hours, minutes, seconds)) + print(" ====================== ") print() print() From 1edfe5c256d39353b57fcd9de07e74e2953f2611 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 11:13:46 +0200 Subject: [PATCH 19/41] Added debugging prints for the Lanczos with degeneracies --- Modules/DynamicalLanczos.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 6d9c7f692..1ab0bf380 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1817,13 +1817,15 @@ def FastApplyD3ToVector(X, Y, rho, w, T, input_vector, symmetries, n_degeneracie """ n_modes = len(w) output_dyn = np.zeros(n_modes*n_modes, dtype = TYPE_DP) - #print( "Apply to vector, nmodes:", n_modes, "shape:", np.shape(output_dyn)) + print( "Apply to vector, nmodes:", n_modes, "shape:", np.shape(output_dyn)) deg_space_new = np.zeros(np.sum(n_degeneracies), dtype = np.intc) i = 0 i_mode = 0 j_mode = 0 + print("Mapping degeneracies:", np.sum(n_degeneracies)) while i_mode < len(deg_space_new): + print("cross_modes: ({}, {}) | deg_i = {}".format(i_mode, j_mode, n_degeneracies[i_mode])) deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 i += 1 From 7163b2a40b4dea8dea1a4f21f2b951749e25fa5f Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 11:17:26 +0200 Subject: [PATCH 20/41] Added other debug prints --- Modules/DynamicalLanczos.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 1ab0bf380..42f6a4ea1 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1749,16 +1749,19 @@ def FastApplyD3ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, de n_modes = len(w) output_vector = np.zeros(n_modes, dtype = TYPE_DP) - #print( "Apply to dyn, nmodes:", n_modes, "shape:", np.shape(output_vector)) + print( "Apply to dyn, nmodes:", n_modes, "shape:", np.shape(output_vector)) deg_space_new = np.zeros(np.sum(n_degeneracies), dtype = np.intc) i = 0 i_mode = 0 j_mode = 0 - #print("len1 = ", len(deg_space_new), "len2 = ", len(n_degeneracies)) + print("len1 = ", len(deg_space_new), "len2 = ", len(n_degeneracies)) + print("Mapping degeneracies:", np.sum(n_degeneracies)) while i_mode < len(deg_space_new): #print("i= ", i_mode, "Ndeg:", n_degeneracies[i_mode], "j = ", j_mode, "len = ", len(degenerate_space[i_mode])) #print("new_i = ", i, "tot = ", np.sum(n_degeneracies)) + print("cross_modes: ({}, {}) | deg_imu = {} | i = {}".format(i_mode, j_mode, n_degeneracies[i_mode], i)) + deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 i+=1 From f6eaf3ebec380fac32d61c47382512b784dfb35b Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 11:23:23 +0200 Subject: [PATCH 21/41] Bugfix in assigning degeneracies --- Modules/DynamicalLanczos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 42f6a4ea1..151cfa012 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1755,9 +1755,9 @@ def FastApplyD3ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, de i = 0 i_mode = 0 j_mode = 0 - print("len1 = ", len(deg_space_new), "len2 = ", len(n_degeneracies)) + print("len1 = ", len(deg_space_new), "len2 = ", n_modes) print("Mapping degeneracies:", np.sum(n_degeneracies)) - while i_mode < len(deg_space_new): + while i_mode < n_modes: #print("i= ", i_mode, "Ndeg:", n_degeneracies[i_mode], "j = ", j_mode, "len = ", len(degenerate_space[i_mode])) #print("new_i = ", i, "tot = ", np.sum(n_degeneracies)) print("cross_modes: ({}, {}) | deg_imu = {} | i = {}".format(i_mode, j_mode, n_degeneracies[i_mode], i)) @@ -1827,7 +1827,7 @@ def FastApplyD3ToVector(X, Y, rho, w, T, input_vector, symmetries, n_degeneracie i_mode = 0 j_mode = 0 print("Mapping degeneracies:", np.sum(n_degeneracies)) - while i_mode < len(deg_space_new): + while i_mode < n_modes: print("cross_modes: ({}, {}) | deg_i = {}".format(i_mode, j_mode, n_degeneracies[i_mode])) deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 From bf292c12fa159c55e166aaf78e5794f5fc1c70f7 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 12:55:53 +0200 Subject: [PATCH 22/41] Bugfix in timing --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 151cfa012..e3b8d652a 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1342,7 +1342,7 @@ def get_all_green_functions(self, N_steps = 100, mode_mixing = True, save_step_d minutes = int(total_time / 60) hours = int(minutes / 60) minutes -= hours * 60 - seconds = total_time - hours*3600 - minutes * 60 + seconds = int(total_time - hours*3600 - minutes * 60) if verbose: print() From 034c636482d8b9a328f4b391ef15befeca0b2679 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 18 Jul 2019 12:56:25 +0200 Subject: [PATCH 23/41] Deleted some debuggings --- Modules/DynamicalLanczos.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index e3b8d652a..cb5d6ba0c 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1749,18 +1749,18 @@ def FastApplyD3ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, de n_modes = len(w) output_vector = np.zeros(n_modes, dtype = TYPE_DP) - print( "Apply to dyn, nmodes:", n_modes, "shape:", np.shape(output_vector)) + #print( "Apply to dyn, nmodes:", n_modes, "shape:", np.shape(output_vector)) deg_space_new = np.zeros(np.sum(n_degeneracies), dtype = np.intc) i = 0 i_mode = 0 j_mode = 0 - print("len1 = ", len(deg_space_new), "len2 = ", n_modes) - print("Mapping degeneracies:", np.sum(n_degeneracies)) + #print("len1 = ", len(deg_space_new), "len2 = ", n_modes) + #print("Mapping degeneracies:", np.sum(n_degeneracies)) while i_mode < n_modes: #print("i= ", i_mode, "Ndeg:", n_degeneracies[i_mode], "j = ", j_mode, "len = ", len(degenerate_space[i_mode])) #print("new_i = ", i, "tot = ", np.sum(n_degeneracies)) - print("cross_modes: ({}, {}) | deg_imu = {} | i = {}".format(i_mode, j_mode, n_degeneracies[i_mode], i)) + #print("cross_modes: ({}, {}) | deg_imu = {} | i = {}".format(i_mode, j_mode, n_degeneracies[i_mode], i)) deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 @@ -1820,15 +1820,15 @@ def FastApplyD3ToVector(X, Y, rho, w, T, input_vector, symmetries, n_degeneracie """ n_modes = len(w) output_dyn = np.zeros(n_modes*n_modes, dtype = TYPE_DP) - print( "Apply to vector, nmodes:", n_modes, "shape:", np.shape(output_dyn)) + #print( "Apply to vector, nmodes:", n_modes, "shape:", np.shape(output_dyn)) deg_space_new = np.zeros(np.sum(n_degeneracies), dtype = np.intc) i = 0 i_mode = 0 j_mode = 0 - print("Mapping degeneracies:", np.sum(n_degeneracies)) + #print("Mapping degeneracies:", np.sum(n_degeneracies)) while i_mode < n_modes: - print("cross_modes: ({}, {}) | deg_i = {}".format(i_mode, j_mode, n_degeneracies[i_mode])) + #print("cross_modes: ({}, {}) | deg_i = {}".format(i_mode, j_mode, n_degeneracies[i_mode])) deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 i += 1 From a592437f1a02f8fbf37c9e7136243529054ba841 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 11:51:23 +0200 Subject: [PATCH 24/41] edited a print --- Modules/Ensemble.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index 1219b3743..7bff30691 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -2063,7 +2063,7 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v Returns ------- - phi_sc : Phonons() + phi_hessian : Phonons() The dynamical matrix of the free energy hessian in (Ry/bohr^2) """ # For now the v4 is not implemented From 54a3c6bba2412d8992aa86981a72798ba640b605 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 12:38:50 +0200 Subject: [PATCH 25/41] Try to use the stored d3 to compute the full L in this way we can compare the C implementation of symmetries with the fortran implementation of d3 --- Modules/DynamicalLanczos.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index cb5d6ba0c..57a39d76d 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -962,7 +962,8 @@ def run(self, n_iter, save_dir = ".", verbose = True): # Apply the matrix L t1 = time.time() - self.psi = self.apply_full_L() + #self.psi = self.apply_full_L() + self.psi = self.L_linop.dot(self.psi) t2 = time.time() if verbose: @@ -1496,6 +1497,10 @@ def get_full_L_operator(self, verbose = False, only_pert=False): d3 += np.einsum("ai,bi,ci", new_Y, new_X, new_X) d3 /= - 3 * N_eff + if verbose: + np.save("d3_modes_nosym.npy", d3) + + # Reshape the d3 d3_reshaped = d3.reshape((self.n_modes, self.n_modes * self.n_modes)) @@ -1513,6 +1518,9 @@ def get_full_L_operator(self, verbose = False, only_pert=False): d4 += np.einsum("ai,bi,ci,di", new_Y, new_X, new_X, new_X) d4 /= - 4 * N_eff + if verbose: + np.save("d4_modes_nosym.npy", d4) + # Reshape the d4 d4_reshaped = d4.reshape((self.n_modes*self.n_modes, self.n_modes * self.n_modes)) From 153b5fb18cfde76be07c533e6a69f7adecc3e4f6 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 15:08:11 +0200 Subject: [PATCH 26/41] Added the d3 calculation directly in mode space --- Modules/Ensemble.py | 60 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index 7bff30691..376594d7a 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -1768,7 +1768,65 @@ def get_average_stress(self): # #TODO: apply symmetries # return df_dfc, err_df_dfc - + + def get_d3_muspace(self): + r""" + GET V3 IN MODE SPACE + ==================== + + This subroutine gets the d3 directly in the space of the modes. + + ..math:: + + D^{(3)}_{abc} = \sum_{xyz} \frac{\Phi^{(3)}_{xyz} e_a^x e_b^y e_c^z}{\sqrt{m_x m_y m_z}} + + + """ + + # Be shure to have the correct units + self.convert_units(UNITS_DEFAULT) + + # Convert from A to Bohr the space + u_disps = self.u_disps * __A_TO_BOHR__ + n_rand, n_modes = np.shape(u_disps) + forces = (self.forces - self.sscha_forces).reshape(self.N, n_modes) / __A_TO_BOHR__ + + Ups = self.current_dyn.GetUpsilonMatrix(self.current_T) + v_disp = u_disps.dot(Ups) + + # pass in the polarization space + w, pols = self.current_dyn.DiagonalizeSupercell() + + # Discard translations + trans = CC.Methods.get_translations(pols, self.structures[0].get_masses_array()) + pols = pols[:, ~trans] + + m = np.tile(self.structures[0].get_masses_array(), (3,1)).T.ravel() + + pol_vec = np.einsum("ab, a->ab", pols, 1 / np.sqrt(m)) + + v_mode = v_disp.dot(pol_vec) + f_mode = forces.dot(pol_vec) + + # Now compute the d3 as + N_eff = np.sum(self.rho) + f_mode = np.einsum("ia, i->ia", f_mode, self.rho) + d3_noperm = np.einsum("ia,ib,ic->abc", v_mode, v_mode, f_mode) + d3_noperm /= N_eff + + # Apply the permuatations + d3 = d3_noperm.copy() + d3 += np.einsum("abc->acb", d3_noperm) + d3 += np.einsum("abc->bac", d3_noperm) + d3 += np.einsum("abc->bca", d3_noperm) + d3 += np.einsum("abc->cab", d3_noperm) + d3 += np.einsum("abc->cba", d3_noperm) + d3 /= 6 + + # TODO: symmetrize + + return d3 + def get_v3_realspace(self): """ This is a testing function that computes the V3 matrix in real space: From f022cba214c23eaf3122bc539e9645405eaf9408 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 15:13:40 +0200 Subject: [PATCH 27/41] Bugfix in the new function --- Modules/Ensemble.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index 376594d7a..5c43df8fc 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -1786,6 +1786,8 @@ def get_d3_muspace(self): # Be shure to have the correct units self.convert_units(UNITS_DEFAULT) + supersturct = self.current_dyn.structure.generate_supercell(self.supercell) + # Convert from A to Bohr the space u_disps = self.u_disps * __A_TO_BOHR__ n_rand, n_modes = np.shape(u_disps) @@ -1798,10 +1800,10 @@ def get_d3_muspace(self): w, pols = self.current_dyn.DiagonalizeSupercell() # Discard translations - trans = CC.Methods.get_translations(pols, self.structures[0].get_masses_array()) + trans = CC.Methods.get_translations(pols, supersturct.get_masses_array()) pols = pols[:, ~trans] - m = np.tile(self.structures[0].get_masses_array(), (3,1)).T.ravel() + m = np.tile(supersturct.get_masses_array(), (3,1)).T.ravel() pol_vec = np.einsum("ab, a->ab", pols, 1 / np.sqrt(m)) From 02d746471fed78506bdb8c48686f8b689cc9cd8c Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 15:16:52 +0200 Subject: [PATCH 28/41] Added a new function to compute the d3 --- Modules/Ensemble.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index 5c43df8fc..fffedbfc8 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -1814,7 +1814,7 @@ def get_d3_muspace(self): N_eff = np.sum(self.rho) f_mode = np.einsum("ia, i->ia", f_mode, self.rho) d3_noperm = np.einsum("ia,ib,ic->abc", v_mode, v_mode, f_mode) - d3_noperm /= N_eff + d3_noperm /= -N_eff # there is a minus # Apply the permuatations d3 = d3_noperm.copy() From ee4ba27681b6b60ad4c8d92f9d7d0e92b00bd080 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 15:35:05 +0200 Subject: [PATCH 29/41] Checking the new d3 computation --- Modules/DynamicalLanczos.py | 38 ++++++++++++++++++++++++++++--------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 57a39d76d..94a800088 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -170,15 +170,24 @@ def __init__(self, ensemble = None, mode = 1): self.rho = ensemble.rho.copy() self.N_eff = np.sum(self.rho) - u = ensemble.u_disps + u = ensemble.u_disps / Ensemble.Bohr f = ensemble.forces.reshape(self.N, 3 * self.nat).copy() f -= ensemble.sscha_forces.reshape(self.N, 3 * self.nat) + f *= Ensemble.Bohr self.X = np.zeros((self.N, self.n_modes), order = order, dtype = TYPE_DP) - self.X[:,:] = np.einsum("a,ia, ab->ib", np.sqrt(self.m), u, self.pols) / Ensemble.Bohr + #self.X[:,:] = np.einsum("a,ia, ab->ib", np.sqrt(self.m), u, self.pols) / Ensemble.Bohr self.Y = np.zeros((self.N, self.n_modes), order = order, dtype = TYPE_DP) - self.Y[:,:] = np.einsum("a,ia, ab->ib", 1/np.sqrt(self.m), f, self.pols) * Ensemble.Bohr + #self.Y[:,:] = np.einsum("a,ia, ab->ib", 1/np.sqrt(self.m), f, self.pols) * Ensemble.Bohr + + ups_mat = self.dyn.GetUpsilonMatrix(self.T) + v_vec = u.dot(ups_mat) + + # Convert in the polarization basis + pol_mat = np.einsum("ab, a->ab", self.pols, 1 / np.sqrt(self.m)) + self.X[:,:] = v_vec.dot(pol_mat) + self.Y[:,:] = f.dot(pol_mat) # Prepare the variable used for the working self.psi = np.zeros(self.n_modes + self.n_modes*self.n_modes, dtype = TYPE_DP) @@ -1485,17 +1494,28 @@ def get_full_L_operator(self, verbose = False, only_pert=False): # Compute the d3 operator - new_X = np.einsum("ia,a->ai", self.X, f_ups(self.w, self.T)) - new_Y = np.einsum("ia,i->ai", self.Y, self.rho) + #new_X = np.einsum("ia,a->ai", self.X, f_ups(self.w, self.T)) + N_eff = np.sum(self.rho) + Y_weighted = np.einsum("ia, i->ia", self.Y, self.rho) + #new_Y = np.einsum("ia,i->ai", self.Y, self.rho) N_eff = np.sum(self.rho) if not self.ignore_v3: if verbose: print("Computing d3...") - d3 = np.einsum("ai,bi,ci", new_X, new_X, new_Y) - d3 += np.einsum("ai,bi,ci", new_X, new_Y, new_X) - d3 += np.einsum("ai,bi,ci", new_Y, new_X, new_X) - d3 /= - 3 * N_eff + d3_noperm = np.einsum("ia,ib,ic->abc", self.X, self.X, Y_weighted) + d3_noperm /= -N_eff + + # Apply the permuatations + d3 = d3_noperm.copy() + d3 += np.einsum("abc->acb", d3_noperm) + d3 += np.einsum("abc->bac", d3_noperm) + d3 += np.einsum("abc->bca", d3_noperm) + d3 += np.einsum("abc->cab", d3_noperm) + d3 += np.einsum("abc->cba", d3_noperm) + d3 /= 6 + + # TODO, symmetrization if verbose: np.save("d3_modes_nosym.npy", d3) From f441bf0ca8ca71f6ed4ae328a30762b9b2c5b102 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 15:47:06 +0200 Subject: [PATCH 30/41] Added the symmetrization --- Modules/DynamicalLanczos.py | 41 ++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 94a800088..3d84400b3 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1515,7 +1515,8 @@ def get_full_L_operator(self, verbose = False, only_pert=False): d3 += np.einsum("abc->cba", d3_noperm) d3 /= 6 - # TODO, symmetrization + # Perform the standard symmetrization + d3 = symmetrize_d3_muspace(d3, self.symmetries) if verbose: np.save("d3_modes_nosym.npy", d3) @@ -2111,3 +2112,41 @@ def GetFreeEnergyCurvatureFromContinuedFraction(a_ns, b_ns, pols_sc, masses, mod fc_odd = np.einsum("ab, ca, da ->cd", fc_pols, epols_m, epols_m) return fc_odd + + +def symmetrize_d3_muspace(d3, symmetries): + """ + SYMMETRIZE D3 IN MODE SPACE + =========================== + + This function symmetrizes the d3 in the mu space. + It is quite fast. + + Parameters + ---------- + d3 : ndarray(n_modes, n_modes, n_modes) + The d3 tensor to be symmetrized + symmetries : ndarray(N_sym, n_modes, n_modes) + The full symmetry matrix + + Results + ------- + new_d3 : ndarray(n_modes, n_modes, n_modes) + The d3 tensor symmetrized + """ + + new_d3 = d3.copy() + + N_sym, nmode, dumb = np.shape(symmetries) + + for i in range(N_sym): + symmat = symmetries[i, :, :] + + ap = np.einsum("abc, lc ->abl", d3, symmat) + ap = np.einsum("abc, lb ->alc", ap, symmat) + ap = np.einsum("abc, la ->lbc", ap, symmat) + + new_d3 += ap + + new_d3 /= N_sym + return new_d3 \ No newline at end of file From 31a0191aabfd972fbf281d411398fb504d21818f Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 16:00:12 +0200 Subject: [PATCH 31/41] Added other debugging prints in the Lanczos --- Modules/DynamicalLanczos.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 3d84400b3..f65cfadfd 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1515,11 +1515,15 @@ def get_full_L_operator(self, verbose = False, only_pert=False): d3 += np.einsum("abc->cba", d3_noperm) d3 /= 6 + if verbose: + np.save("d3_modes_nosym.npy", d3) + # Perform the standard symmetrization d3 = symmetrize_d3_muspace(d3, self.symmetries) if verbose: - np.save("d3_modes_nosym.npy", d3) + np.save("d3_modes_sym.npy", d3) + np.save("symmetries_modes.npy", self.symmetries) # Reshape the d3 From 325b791ed70933df6fe1e91842043b41fff30cbe Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Fri, 19 Jul 2019 16:11:33 +0200 Subject: [PATCH 32/41] bugfix in mode symmetrization --- Modules/DynamicalLanczos.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index f65cfadfd..ebb8a419e 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2139,7 +2139,7 @@ def symmetrize_d3_muspace(d3, symmetries): The d3 tensor symmetrized """ - new_d3 = d3.copy() + new_d3 = np.zeros(np.shape(d3), dtype = np.double) N_sym, nmode, dumb = np.shape(symmetries) @@ -2149,6 +2149,7 @@ def symmetrize_d3_muspace(d3, symmetries): ap = np.einsum("abc, lc ->abl", d3, symmat) ap = np.einsum("abc, lb ->alc", ap, symmat) ap = np.einsum("abc, la ->lbc", ap, symmat) + #ap = np.einsum("abc, aa, bb, cc->abc", d3, symmat, symmat, symmat) new_d3 += ap From 5b1b98812932f3a8412a459e4dc32fe2e20cfd24 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 23 Jul 2019 15:44:25 +0200 Subject: [PATCH 33/41] Debugging print --- Modules/DynamicalLanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index ebb8a419e..f63cc6909 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1485,7 +1485,7 @@ def get_full_L_operator(self, verbose = False, only_pert=False): w_a = np.tile(self.w, (self.n_modes,1)).ravel() w_b = np.tile(self.w, (self.n_modes,1)).T.ravel() - chi_beta = -.5 * np.sqrt((w_a + w_b)/(w_a*w_b)) + chi_beta = -.5 * np.sqrt(w_a + w_b)/(np.sqrt(w_a)*np.sqrt(w_b)) B_mat = (w_a + w_b)**2 From 0355d3e58ffa4cc20d7adb273f27b4cca7e7a6b2 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 23 Jul 2019 16:54:32 +0200 Subject: [PATCH 34/41] Reverted symmetries --- Modules/DynamicalLanczos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index f63cc6909..b7516ff8c 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2146,9 +2146,9 @@ def symmetrize_d3_muspace(d3, symmetries): for i in range(N_sym): symmat = symmetries[i, :, :] - ap = np.einsum("abc, lc ->abl", d3, symmat) - ap = np.einsum("abc, lb ->alc", ap, symmat) - ap = np.einsum("abc, la ->lbc", ap, symmat) + ap = np.einsum("abc, cl ->abl", d3, symmat) + ap = np.einsum("abc, bl ->alc", ap, symmat) + ap = np.einsum("abc, al ->lbc", ap, symmat) #ap = np.einsum("abc, aa, bb, cc->abc", d3, symmat, symmat, symmat) new_d3 += ap From 7db512c7a149114797774851948e763c4b21b1cf Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Tue, 23 Jul 2019 16:55:49 +0200 Subject: [PATCH 35/41] Back to good symmetries --- Modules/DynamicalLanczos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index b7516ff8c..f63cc6909 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -2146,9 +2146,9 @@ def symmetrize_d3_muspace(d3, symmetries): for i in range(N_sym): symmat = symmetries[i, :, :] - ap = np.einsum("abc, cl ->abl", d3, symmat) - ap = np.einsum("abc, bl ->alc", ap, symmat) - ap = np.einsum("abc, al ->lbc", ap, symmat) + ap = np.einsum("abc, lc ->abl", d3, symmat) + ap = np.einsum("abc, lb ->alc", ap, symmat) + ap = np.einsum("abc, la ->lbc", ap, symmat) #ap = np.einsum("abc, aa, bb, cc->abc", d3, symmat, symmat, symmat) new_d3 += ap From 6be004ddc61c6ceea8cddacc080a30a64cc0652e Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 24 Jul 2019 12:22:15 +0200 Subject: [PATCH 36/41] Edited the fast lanczos functions for the clever usage of symmetries --- CModules/LanczosFunctions.c | 413 ++++++++++++++++++++++++++++++++++-- 1 file changed, 395 insertions(+), 18 deletions(-) diff --git a/CModules/LanczosFunctions.c b/CModules/LanczosFunctions.c index bace2fbe1..ebbb2df09 100644 --- a/CModules/LanczosFunctions.c +++ b/CModules/LanczosFunctions.c @@ -13,6 +13,140 @@ double f_ups(double w, double T) { +// void OMP_ApplyD3ToVector(const double * X, const double * Y, const double * rho, const double * w, double T, int N_modes, +// int N_configs, const double * input_vector, double * output_dyn, +// double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space ) { + +// // Compute the N_eff +// double N_eff = 0; +// int i; + +// //#pragma omp parallel for private(i) reduction(+:N_eff) +// for (i = 0; i < N_configs; ++i) +// N_eff += rho[i]; + +// if (DEB) { +// printf("File %s, Line %d: N_eff = %.4f\n", __FILE__, __LINE__, N_eff); +// fflush(stdout); +// } + +// // Prepare the new modified X +// double * new_X = malloc(sizeof(double) * N_configs* N_modes); + +// //#pragma omp parallel for private(i) +// for (i = 0; i < N_configs*N_modes; ++i) { +// new_X[i] = X[i] * f_ups(w[i / N_configs], T); +// } + +// if (DEB) { +// printf("File %s, Line %d: Got the new X\n", __FILE__, __LINE__, N_eff); +// fflush(stdout); +// } + + +// // Initialize the output +// for (i = 0; i < N_modes*N_modes; ++i) +// output_dyn[i] = 0; + +// if (DEB) printf("Applying the d3 to vector!\n"); + +// // Perform the application +// int a, b, c, new_a, new_b, new_c; +// int j, k, i_sym, N_sym_tmp; +// double sym_coeff = 0; +// //#pragma omp parallel for collapse(3) private(a,b,c) +// for (a = 0; a < N_modes; ++a) { +// for (b = 0; b < N_modes; ++b) { +// for (c = 0; c < N_modes; ++c) { +// // Check if this element is zero by symmetry +// int stop= 0; + +// if (DEB) printf("Element a=%d, b=%d, c=%d ... \n", a, b, c); +// for (i = 0; i < N_sym; ++i) { +// if (DEB) printf("%d) s_aa = %.2f, s_bb = %.2f, s_cc = %.2f\n", i, +// symmetries[i * N_modes*N_modes + a*N_modes + a], +// symmetries[i * N_modes*N_modes + b*N_modes + b], +// symmetries[i * N_modes*N_modes + c*N_modes + c]); +// if (fabs(symmetries[i * N_modes*N_modes + a*N_modes + a] * +// symmetries[i * N_modes*N_modes + b*N_modes + b] * +// symmetries[i * N_modes*N_modes + c*N_modes + c] + 1) < __EPSILON__) { +// stop = 1; +// break; +// } +// } +// if (stop == 1) continue; +// if (DEB)printf("I'm computing this element.\n"); + +// double tmp = 0; + +// //#pragma omp parallel for private(i) reduce(+:tmp) +// for (i = 0; i < N_configs; ++i) { +// tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; +// //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i]; +// //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i]; + +// } + + +// // Apply all the symmetries in the degenerate subspace +// for (i = 0; i < N_degeneracy[a]; ++i) { +// new_a = degenerate_space[a][i]; +// for (j = 0; j < N_degeneracy[b]; ++j) { +// new_b = degenerate_space[b][j]; +// for (k = 0; k < N_degeneracy[c]; ++k) { +// new_c = degenerate_space[c][k]; + +// // Check if there are degeneracies +// // If not, symmetries are useless, apply only the identity +// N_sym_tmp = N_sym; +// if (N_degeneracy[a] * N_degeneracy[b] * N_degeneracy[c] == 1) +// N_sym_tmp = 1; + +// if (DEB) +// printf("Deg space = %d | new_a = %d, new_b = %d, new_c = %d\n", N_degeneracy[a] * N_degeneracy[b] * N_degeneracy[c], +// new_a, new_b, new_c); + +// for (i_sym = 0; i_sym < N_sym_tmp; ++i_sym) { +// sym_coeff = symmetries[i_sym * N_modes * N_modes + a * N_modes + new_a] * +// symmetries[i_sym * N_modes * N_modes + b * N_modes + new_b] * +// symmetries[i_sym * N_modes * N_modes + c * N_modes + new_c]; + +// if (DEB) +// printf("IN_VEC_OUT_DYN: symfactor = %.2f | d3[%d, %d, %d] = %.6e\n", sym_coeff, a, b, c, -tmp / (N_eff)); + + + +// output_dyn[new_a * N_modes + new_b] += -tmp * input_vector[new_c] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_dyn[new_b * N_modes + new_a] += -tmp * input_vector[new_c] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_dyn[new_a * N_modes + new_c] += -tmp * input_vector[new_b] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_dyn[new_c * N_modes + new_a] += -tmp * input_vector[new_b] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_dyn[new_c * N_modes + new_b] += -tmp * input_vector[new_a] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_dyn[new_b * N_modes + new_c] += -tmp * input_vector[new_a] * sym_coeff / (6 * N_eff * N_sym_tmp); +// } +// } +// } +// } +// if (DEB) +// printf("\n"); + +// // Check if the symmetries are mixing something + +// if (DEB && c == N_modes -1) { +// printf("a = %d, b = %d: output = %.8e\n", a, b, output_dyn[a * N_modes + b]); +// } +// } + +// } +// } + +// // Free memory +// free(new_X); +// } + + +/* + * This function uses the new defintions of the X and Y matrices + */ void OMP_ApplyD3ToVector(const double * X, const double * Y, const double * rho, const double * w, double T, int N_modes, int N_configs, const double * input_vector, double * output_dyn, double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space ) { @@ -31,12 +165,12 @@ void OMP_ApplyD3ToVector(const double * X, const double * Y, const double * rho, } // Prepare the new modified X - double * new_X = malloc(sizeof(double) * N_configs* N_modes); + //double * new_X = malloc(sizeof(double) * N_configs* N_modes); //#pragma omp parallel for private(i) - for (i = 0; i < N_configs*N_modes; ++i) { - new_X[i] = X[i] * f_ups(w[i / N_configs], T); - } + //for (i = 0; i < N_configs*N_modes; ++i) { + // new_X[i] = X[i] * f_ups(w[i / N_configs], T); + //} if (DEB) { printf("File %s, Line %d: Got the new X\n", __FILE__, __LINE__, N_eff); @@ -81,7 +215,7 @@ void OMP_ApplyD3ToVector(const double * X, const double * Y, const double * rho, //#pragma omp parallel for private(i) reduce(+:tmp) for (i = 0; i < N_configs; ++i) { - tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; + tmp += X[N_configs*a + i] * X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i]; //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i]; @@ -140,11 +274,117 @@ void OMP_ApplyD3ToVector(const double * X, const double * Y, const double * rho, } // Free memory - free(new_X); + //free(new_X); } + +// void OMP_ApplyD3ToDyn(const double * X, const double * Y, const double * rho, const double * w, +// double T, int N_modes, int N_configs, const double * input_dyn, double * output_vector, +// double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space){ + +// // Compute the N_eff +// double N_eff = 0; +// int i; + +// //#pragma omp parallel for private(i) reduction(+:N_eff) +// for (i = 0; i < N_configs; ++i) +// N_eff += rho[i]; + + +// // Prepare the new modified X +// double * new_X = malloc(sizeof(double) * N_configs* N_modes); + +// //#pragma omp parallel for private(i) +// for (i = 0; i < N_configs*N_modes; ++i) { +// new_X[i] = X[i] * f_ups(w[i / N_configs], T); +// } + +// // Initialize the output +// for (i = 0; i < N_modes; ++i) +// output_vector[i] = 0; + +// // Perform the application +// int a, b, c, new_a, new_b, new_c; +// int j, k, i_sym, N_sym_tmp; +// double sym_coeff = 0; +// //#pragma omp parallel for collapse(3) private(a,b,c) +// for (a = 0; a < N_modes; ++a) { +// for (b = 0; b < N_modes; ++b) { +// for (c = 0; c < N_modes; ++c) { +// double tmp = 0; + + +// // Check if this element is zero by symmetry +// int stop= 0; +// for (i = 0; i < N_sym; ++i) { +// if (fabs(symmetries[i * N_modes*N_modes + a*N_modes + a] * +// symmetries[i * N_modes*N_modes + b*N_modes + b] * +// symmetries[i * N_modes*N_modes + c*N_modes + c] + 1) < __EPSILON__) { +// stop = 1; +// break; +// } +// } + +// if (stop == 1) continue; + +// for (i = 0; i < N_configs; ++i) { +// tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; +// //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i]; +// //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i]; +// //tmp += tmp1 * rho[i]; +// } + +// // Apply all the symmetries in the degenerate subspace +// for (i = 0; i < N_degeneracy[a]; ++i) { +// new_a = degenerate_space[a][i]; +// for (j = 0; j < N_degeneracy[b]; ++j) { +// new_b = degenerate_space[b][j]; +// for (k = 0; k < N_degeneracy[c]; ++k) { +// new_c = degenerate_space[c][k]; + +// // Check if there are degeneracies +// // If not, symmetries are useless, apply only the identity +// N_sym_tmp = N_sym; +// if (N_degeneracy[a] * N_degeneracy[b] * N_degeneracy[c] == 1) +// N_sym_tmp = 1; + +// for (i_sym = 0; i_sym < N_sym_tmp; ++i_sym) { +// sym_coeff = symmetries[i_sym * N_modes * N_modes + a * N_modes + new_a] * +// symmetries[i_sym * N_modes * N_modes + b * N_modes + new_b] * +// symmetries[i_sym * N_modes * N_modes + c * N_modes + new_c]; + + +// if (DEB) +// printf("IN_DYN_OUT_VEC: symfactor = %.2f | d3[%d, %d, %d] = %.6e\n", sym_coeff, a, b, c, -tmp / (N_eff)); + + +// output_vector[new_a] += -tmp * input_dyn[N_modes * new_b + new_c] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_vector[new_a] += -tmp * input_dyn[N_modes * new_c + new_b] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_vector[new_b] += -tmp * input_dyn[N_modes * new_c + new_a] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_vector[new_b] += -tmp * input_dyn[N_modes * new_a + new_c] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_vector[new_c] += -tmp * input_dyn[N_modes * new_a + new_b] * sym_coeff / (6 * N_eff * N_sym_tmp); +// output_vector[new_c] += -tmp * input_dyn[N_modes * new_b + new_a] * sym_coeff / (6 * N_eff * N_sym_tmp); + +// } +// } +// } +// } + + +// if (DEB && b == N_modes - 1 && c == N_modes - 1) { +// printf("a = %d, output = %.8e\n", a, output_vector[a]); +// } + +// } +// } +// } + +// // Free memory +// free(new_X); +// } + void OMP_ApplyD3ToDyn(const double * X, const double * Y, const double * rho, const double * w, double T, int N_modes, int N_configs, const double * input_dyn, double * output_vector, double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space){ @@ -159,12 +399,12 @@ void OMP_ApplyD3ToDyn(const double * X, const double * Y, const double * rho, co // Prepare the new modified X - double * new_X = malloc(sizeof(double) * N_configs* N_modes); + //double * new_X = malloc(sizeof(double) * N_configs* N_modes); //#pragma omp parallel for private(i) - for (i = 0; i < N_configs*N_modes; ++i) { - new_X[i] = X[i] * f_ups(w[i / N_configs], T); - } + //for (i = 0; i < N_configs*N_modes; ++i) { + // new_X[i] = X[i] * f_ups(w[i / N_configs], T); + //} // Initialize the output for (i = 0; i < N_modes; ++i) @@ -195,7 +435,7 @@ void OMP_ApplyD3ToDyn(const double * X, const double * Y, const double * rho, co if (stop == 1) continue; for (i = 0; i < N_configs; ++i) { - tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; + tmp += X[N_configs*a + i] * X[N_configs*b + i] * Y[N_configs*c +i] * rho[i]; //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i]; //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i]; //tmp += tmp1 * rho[i]; @@ -247,11 +487,148 @@ void OMP_ApplyD3ToDyn(const double * X, const double * Y, const double * rho, co } // Free memory - free(new_X); + //free(new_X); } +// void OMP_ApplyD4ToDyn(const double * X, const double * Y, const double * rho, const double * w, +// double T, int N_modes, int N_configs, const double * input_dyn, double * output_dyn, +// double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space){ + +// // Compute the N_eff +// double N_eff = 0; +// int i; + +// #pragma omp parallel for private(i) reduction(+:N_eff) +// for (i = 0; i < N_configs; ++i) +// N_eff += rho[i]; + + +// // Prepare the new modified X +// double * new_X = malloc(sizeof(double) * N_configs* N_modes); + +// //#pragma omp parallel for private(i) +// for (i = 0; i < N_configs*N_modes; ++i) { +// new_X[i] = X[i] * f_ups(w[i / N_configs], T); +// } + +// // Initialize the output +// for (i = 0; i < N_modes*N_modes; ++i) +// output_dyn[i] = 0; + +// // Perform the application +// int a, b, c, d; +// int new_a, new_b, new_c, new_d; +// int j, k, h, i_sym, N_sym_tmp; +// double sym_coeff = 0; + +// //#pragma omp parallel for collapse(4) private(a,b,c, d) +// for (a = 0; a < N_modes; ++a) { +// for (b = 0; b < N_modes; ++b) { +// for (c = 0; c < N_modes; ++c) { +// for (d = 0; d < N_modes; ++d) { + +// // Check if the operation is allowed by symmetries +// int stop = 0; + +// for (i = 0; i < N_sym; ++i) { +// if (fabs(symmetries[i * N_modes*N_modes + a*N_modes + a] * +// symmetries[i * N_modes*N_modes + b*N_modes + b] * +// symmetries[i * N_modes*N_modes + c*N_modes + c] * +// symmetries[i * N_modes*N_modes + d*N_modes + d] + 1) < __EPSILON__) { +// stop = 1; +// break; +// } +// } + +// if (stop == 1) continue; + + +// double tmp = 0; +// for (i = 0; i < N_configs; ++i) { +// tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i] * Y[N_configs*d +i] * rho[i]; +// //tmp1 += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * new_X[N_configs*d +i]; +// //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i] * new_X[N_configs*d +i]; +// //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i] * new_X[N_configs*d +i]; +// //tmp += tmp1 * rho[i]; +// } + + +// // Apply all the symmetries in the degenerate subspace +// for (i = 0; i < N_degeneracy[a]; ++i) { +// new_a = degenerate_space[a][i]; +// for (j = 0; j < N_degeneracy[b]; ++j) { +// new_b = degenerate_space[b][j]; +// for (k = 0; k < N_degeneracy[c]; ++k) { +// new_c = degenerate_space[c][k]; +// for (h = 0; h < N_degeneracy[d]; ++h) { +// new_d = degenerate_space[d][h]; + +// // Check if there are degeneracies +// // If not, symmetries are useless, apply only the identity +// N_sym_tmp = N_sym; +// if (N_degeneracy[a] * N_degeneracy[b] * N_degeneracy[c] * N_degeneracy[d] == 1) +// N_sym_tmp = 1; + +// for (i_sym = 0; i_sym < N_sym_tmp; ++i_sym) { +// sym_coeff = symmetries[i_sym * N_modes * N_modes + a * N_modes + new_a] * +// symmetries[i_sym * N_modes * N_modes + b * N_modes + new_b] * +// symmetries[i_sym * N_modes * N_modes + c * N_modes + new_c] * +// symmetries[i_sym * N_modes * N_modes + d * N_modes + new_d]; + + +// if (DEB) +// printf("IN_DYN_OUT_DYN: symfactor = %.2f | d4[%d, %d, %d, %d] = %.6e\n", sym_coeff, a, b, c, d, -tmp / (N_eff)); + + +// output_dyn[N_modes * new_a + new_b] += -tmp * input_dyn[N_modes * new_c + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_b + new_a] += -tmp * input_dyn[N_modes * new_c + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_a + new_b] += -tmp * input_dyn[N_modes * new_d + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_b + new_a] += -tmp * input_dyn[N_modes * new_d + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); + +// output_dyn[N_modes * new_a + new_c] += -tmp * input_dyn[N_modes * new_b + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_c + new_a] += -tmp * input_dyn[N_modes * new_b + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_a + new_c] += -tmp * input_dyn[N_modes * new_d + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_c + new_a] += -tmp * input_dyn[N_modes * new_d + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); + +// output_dyn[N_modes * new_a + new_d] += -tmp * input_dyn[N_modes * new_b + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_a] += -tmp * input_dyn[N_modes * new_b + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_a + new_d] += -tmp * input_dyn[N_modes * new_c + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_a] += -tmp * input_dyn[N_modes * new_c + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); + +// output_dyn[N_modes * new_b + new_c] += -tmp * input_dyn[N_modes * new_a + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_c + new_b] += -tmp * input_dyn[N_modes * new_a + new_d] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_b + new_c] += -tmp * input_dyn[N_modes * new_d + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_c + new_b] += -tmp * input_dyn[N_modes * new_d + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); + +// output_dyn[N_modes * new_b + new_d] += -tmp * input_dyn[N_modes * new_a + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_b] += -tmp * input_dyn[N_modes * new_a + new_c] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_b + new_d] += -tmp * input_dyn[N_modes * new_c + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_b] += -tmp * input_dyn[N_modes * new_c + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); + +// output_dyn[N_modes * new_c + new_d] += -tmp * input_dyn[N_modes * new_a + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_c] += -tmp * input_dyn[N_modes * new_a + new_b] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_c + new_d] += -tmp * input_dyn[N_modes * new_b + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); +// output_dyn[N_modes * new_d + new_c] += -tmp * input_dyn[N_modes * new_b + new_a] * sym_coeff / (24 * N_eff * N_sym_tmp); +// } +// } +// } +// } +// } +// //output_dyn[a*N_modes + b] += -tmp * input_dyn[N_modes * b + c] / (4 * N_eff); +// } +// } +// } +// } + +// // Free memory +// free(new_X); +// } + + + + void OMP_ApplyD4ToDyn(const double * X, const double * Y, const double * rho, const double * w, double T, int N_modes, int N_configs, const double * input_dyn, double * output_dyn, double * symmetries, int N_sym, int * N_degeneracy, int ** degenerate_space){ @@ -266,12 +643,12 @@ void OMP_ApplyD4ToDyn(const double * X, const double * Y, const double * rho, co // Prepare the new modified X - double * new_X = malloc(sizeof(double) * N_configs* N_modes); + //double * new_X = malloc(sizeof(double) * N_configs* N_modes); //#pragma omp parallel for private(i) - for (i = 0; i < N_configs*N_modes; ++i) { - new_X[i] = X[i] * f_ups(w[i / N_configs], T); - } + //for (i = 0; i < N_configs*N_modes; ++i) { + // new_X[i] = X[i] * f_ups(w[i / N_configs], T); + //} // Initialize the output for (i = 0; i < N_modes*N_modes; ++i) @@ -307,7 +684,7 @@ void OMP_ApplyD4ToDyn(const double * X, const double * Y, const double * rho, co double tmp = 0; for (i = 0; i < N_configs; ++i) { - tmp += new_X[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i] * Y[N_configs*d +i] * rho[i]; + tmp += X[N_configs*a + i] * X[N_configs*b + i] * X[N_configs*c +i] * Y[N_configs*d +i] * rho[i]; //tmp1 += new_X[N_configs*a + i] * new_X[N_configs*b + i] * Y[N_configs*c +i] * new_X[N_configs*d +i]; //tmp1 += new_X[N_configs*a + i] * Y[N_configs*b + i] * new_X[N_configs*c +i] * new_X[N_configs*d +i]; //tmp1 += Y[N_configs*a + i] * new_X[N_configs*b + i] * new_X[N_configs*c +i] * new_X[N_configs*d +i]; @@ -383,7 +760,7 @@ void OMP_ApplyD4ToDyn(const double * X, const double * Y, const double * rho, co } // Free memory - free(new_X); + //free(new_X); } From be79ccb66de6a7d13bb0a5ac92a8d45e40840465 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 24 Jul 2019 12:33:40 +0200 Subject: [PATCH 37/41] Bugfix in the symmetrization preparation for the v4 computation in Lanczos --- Modules/DynamicalLanczos.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index f63cc6909..943d29c19 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -1989,7 +1989,9 @@ def FastApplyD4ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, deg i = 0 i_mode = 0 j_mode = 0 - while i_mode < len(deg_space_new): + #print("Mapping degeneracies:", np.sum(n_degeneracies)) + while i_mode < n_modes: + #print("cross_modes: ({}, {}) | deg_i = {}".format(i_mode, j_mode, n_degeneracies[i_mode])) deg_space_new[i] = degenerate_space[i_mode][j_mode] j_mode += 1 i += 1 @@ -1997,6 +1999,7 @@ def FastApplyD4ToDyn(X, Y, rho, w, T, input_dyn, symmetries, n_degeneracies, deg i_mode += 1 j_mode = 0 + #print( "Apply to vector, nmodes:", n_modes, "shape:", np.shape(output_dyn)) sscha_HP_odd.ApplyV4ToDyn(X, Y, rho, w, T, input_dyn, output_dyn, mode, symmetries, n_degeneracies, deg_space_new) From 1a99645d412b1d4eb95341799903611f581998bf Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 24 Jul 2019 22:12:55 +0200 Subject: [PATCH 38/41] Added a semplification for the prepare_from_namelist methods --- Modules/Calculator.py | 7 ++++++- Modules/Cluster.py | 3 +++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/Modules/Calculator.py b/Modules/Calculator.py index 18a770afc..a624a54b9 100644 --- a/Modules/Calculator.py +++ b/Modules/Calculator.py @@ -46,14 +46,19 @@ def prepare_calculator_from_namelist(namelist): Parameters ---------- - namelist : dict + namelist : dict or string The parsed namelist. + If a string is passed, it will be parsed Returns ------- ase_calc : The ASE calculator. """ + + # Parse the namelist if needed + if isinstance(namelist, str): + namelist = CC.Methods.read_namelist(namelist) # Check if the namelist has the correct keys if not __CALCULATOR_HEAD__ in namelist.keys(): diff --git a/Modules/Cluster.py b/Modules/Cluster.py index 4ea615208..8206227fb 100644 --- a/Modules/Cluster.py +++ b/Modules/Cluster.py @@ -625,6 +625,9 @@ def setup_from_namelist(self, namelist): namelist: The parsed namelist dictionary. """ + # Parse the namelist if needed + if isinstance(namelist, str): + namelist = CC.Methods.read_namelist(namelist) # Check if the cluster namelist is present if not __CLUSTER_NAMELIST__ in namelist.keys(): From 6765ae84d00c3fb87589870433b41364f1e4dc0c Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Wed, 24 Jul 2019 22:13:27 +0200 Subject: [PATCH 39/41] Added a missing import in Calculator --- Modules/Calculator.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Modules/Calculator.py b/Modules/Calculator.py index a624a54b9..b98082249 100644 --- a/Modules/Calculator.py +++ b/Modules/Calculator.py @@ -1,4 +1,6 @@ # -*- coding: utf-8 -*- +import cellconstructor as CC +import cellconstructor.Methods import difflib __ASE_CALC__ = False From 8fbd0be970f30ed8e6c499ed0d5c1ef79b3d0443 Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 25 Jul 2019 10:59:39 +0200 Subject: [PATCH 40/41] Added the possibility to symmetrize the stress tensor using spglib --- Modules/Ensemble.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/Modules/Ensemble.py b/Modules/Ensemble.py index fffedbfc8..cc5c39d43 100644 --- a/Modules/Ensemble.py +++ b/Modules/Ensemble.py @@ -1411,7 +1411,8 @@ def get_covmat_from_ensemble(self): return cov_mat - def get_stress_tensor(self, offset_stress = None, add_centroid_contrib = False): + def get_stress_tensor(self, offset_stress = None, add_centroid_contrib = False, \ + symmetrize_with_spglib = False): """ GET STRESS TENSOR ================= @@ -1430,6 +1431,8 @@ def get_stress_tensor(self, offset_stress = None, add_centroid_contrib = False): add_centroid_contrib : bool, optional If true the contribution of the centroid is added. This is always zero when the system is relaxed. + symmetrize_with_spglib : bool, optional + If true the symmetrization is performed using spglib Results ------- @@ -1532,6 +1535,10 @@ def get_stress_tensor(self, offset_stress = None, add_centroid_contrib = False): # Symmetrize the stress tensor qe_sym = CC.symmetries.QE_Symmetry(self.current_dyn.structure) + if symmetrize_with_spglib: + qe_sym.SetupFromSPGLIB() + else: + qe_sym.SetupQPoint() qe_sym.ApplySymmetryToMatrix(stress, err_stress) return stress, err_stress From c4bdb3c37989668d1315d43738e6b6654c75f29c Mon Sep 17 00:00:00 2001 From: Lorenzo Monacelli Date: Thu, 25 Jul 2019 12:07:52 +0200 Subject: [PATCH 41/41] Added the automatic symmetrization with spglib of the stress tensor when the used of the spglib library is requested for symmetrization --- Modules/SchaMinimizer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/SchaMinimizer.py b/Modules/SchaMinimizer.py index d571ceb70..f928acc3d 100644 --- a/Modules/SchaMinimizer.py +++ b/Modules/SchaMinimizer.py @@ -1268,7 +1268,7 @@ def get_stress_tensor(self): """ - return self.ensemble.get_stress_tensor(self.stress_offset) + return self.ensemble.get_stress_tensor(self.stress_offset, symmetrize_with_spglib = self.use_spglib) def get_root_dyn(dyn_fc, root_representation):