Skip to content

Commit ce95fd7

Browse files
mesonepigrecomesonepigreco
committed
Now the report seems to work properly!
1 parent 897f857 commit ce95fd7

5 files changed

Lines changed: 261 additions & 96 deletions

File tree

‎Modules/Ensemble.py‎

Lines changed: 74 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1342,7 +1342,7 @@ def _unwrap_symmetries_(self):
13421342

13431343

13441344

1345-
def update_weights(self, new_dynamical_matrix, newT, update_q = False):
1345+
def update_weights(self, new_dynamical_matrix, newT, update_q = False, timer=None):
13461346
"""
13471347
IMPORTANCE SAMPLING
13481348
===================
@@ -1381,7 +1381,10 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False):
13811381
old_disps[:,:] = self.xats.reshape((self.N, 3*Nat_sc)) - np.tile(super_struct0.coords.ravel(), (self.N,1))
13821382

13831383
if changed_dyn:
1384-
w_new, pols = new_dynamical_matrix.DiagonalizeSupercell()#new_super_dyn.DyagDinQ(0)
1384+
if timer:
1385+
w_new, pols = timer.execute_timed_function(new_dynamical_matrix.DiagonalizeSupercell)
1386+
else:
1387+
w_new, pols = new_dynamical_matrix.DiagonalizeSupercell()#new_super_dyn.DyagDinQ(0)
13851388
self.current_w = w_new.copy()
13861389
self.current_pols = pols.copy()
13871390
else:
@@ -1391,7 +1394,11 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False):
13911394
#self.current_w = w_new.copy()
13921395
#self.current_pols = pols.copy()
13931396
# Update sscha energies and forces
1394-
self.sscha_energies[:], self.sscha_forces[:,:,:] = new_dynamical_matrix.get_energy_forces(None, displacement = self.u_disps, w_pols = (w_new, pols))
1397+
if timer:
1398+
self.sscha_energies[:], self.sscha_forces[:,:,:] = timer.execute_timed_function(new_dynamical_matrix.get_energy_forces,
1399+
None, displacement = self.u_disps, w_pols = (w_new, pols))
1400+
else:
1401+
self.sscha_energies[:], self.sscha_forces[:,:,:] = new_dynamical_matrix.get_energy_forces(None, displacement = self.u_disps, w_pols = (w_new, pols))
13951402

13961403

13971404
t1 = time.time()
@@ -1455,15 +1462,16 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False):
14551462

14561463

14571464
# Get the new displacements in the supercell
1458-
t3 = time.time()
1465+
t2 = time.time()
1466+
if timer:
1467+
timer.add_timer("update preparation", t2-t1)
14591468
# for i in range(self.N):
14601469
# self.u_disps[i, :] = (self.xats[i, :, :] - super_structure.coords).reshape( 3*Nat_sc )
14611470

14621471
# old_disps[i,:] = (self.xats[i, :, :] - super_dyn.structure.coords).reshape( 3*Nat_sc )
14631472

14641473
# # TODO: this method recomputes the displacements, it is useless since we already have them in self.u_disps
14651474

1466-
t4 = time.time()
14671475

14681476

14691477
# Convert the q vectors in the Hartree units
@@ -1483,22 +1491,23 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False):
14831491

14841492

14851493
# Get the covariance matrices of the ensemble
1486-
ups_new = np.real(new_dynamical_matrix.GetUpsilonMatrix(self.current_T, w_pols = (w_new, pols)))
1487-
ups_old = np.real(self.dyn_0.GetUpsilonMatrix(self.T0, w_pols = (w_original, pols_original)))
1494+
if timer:
1495+
ups_new = np.real(timer.execute_timed_function(new_dynamical_matrix.GetUpsilonMatrix, self.current_T, w_pols = (w_new, pols)))
1496+
ups_old = np.real(timer.execute_timed_function(self.dyn_0.GetUpsilonMatrix, self.T0, w_pols = (w_original, pols_original)))
1497+
else:
1498+
ups_new = np.real(new_dynamical_matrix.GetUpsilonMatrix(self.current_T, w_pols = (w_new, pols)))
1499+
ups_old = np.real(self.dyn_0.GetUpsilonMatrix(self.T0, w_pols = (w_original, pols_original)))
14881500

14891501
# Get the normalization ratio
14901502
#norm = np.sqrt(np.abs(np.linalg.det(ups_new) / np.linalg.det(ups_old)))
14911503
norm = np.prod( old_a / new_a)
14921504

14931505
t2 = time.time()
1494-
print( "Time elapsed to prepare the rho update:", t2 - t1, "s")
1495-
print (" (of which to diagonalize and prepare the structure: %.4f s)" % (t3-t1))
1496-
print ( "(of which to update sscha energies and forces: %.4f s)" % (t4-t3))
1497-
print ( "(of which computing the Upsilon matrix: %.4f s)" % (t2 - t4))
14981506

14991507
rho_tmp = np.ones( self.N, dtype = np.float64) * norm
15001508
if __DEBUG_RHO__:
15011509
print("Norm factor:", norm)
1510+
15021511
for i in range(self.N):
15031512
v_new = self.u_disps[i, :].dot(ups_new.dot(self.u_disps[i, :])) * __A_TO_BOHR__**2
15041513
v_old = old_disps[i, :].dot(ups_old.dot(old_disps[i, :])) * __A_TO_BOHR__**2
@@ -1524,11 +1533,12 @@ def update_weights(self, new_dynamical_matrix, newT, update_q = False):
15241533
#self.u_disps[i, :] = self.structures[i].get_displacement(new_super_dyn.structure).reshape(3 * new_super_dyn.structure.N_atoms)
15251534
#self.u_disps[i, :] = (self.xats[i, :, :] - super_structure.coords).reshape( 3*Nat_sc )
15261535
t1 = time.time()
1536+
if timer:
1537+
timer.add_timer("<u| Y |u> and weights update (TODO: parallelizable)", t1-t2)
1538+
self.current_dyn = timer.execute_timed_function(new_dynamical_matrix.Copy)
15271539
#print( "Time elapsed to update weights the sscha energies, forces and displacements:", t1 - t3, "s")
1528-
print( "(of which to update the weights):", t1 - t2, "s")
1529-
self.current_dyn = new_dynamical_matrix.Copy()
1530-
t2 = time.time()
1531-
print(" | to copy the dynamical matrix: {} s".format(t2-t1))
1540+
else:
1541+
self.current_dyn = new_dynamical_matrix.Copy()
15321542

15331543

15341544
if __DEBUG_RHO__:
@@ -1880,7 +1890,7 @@ def get_fc_from_self_consistency(self, subtract_sscha = False, return_error = Fa
18801890

18811891
def get_preconditioned_gradient(self, subtract_sscha = True, return_error = False,
18821892
use_ups_supercell = True, preconditioned = 1,
1883-
fast_grad = False, verbose = True):
1893+
fast_grad = False, verbose = True, timer=None):
18841894
"""
18851895
SELF CONSISTENT SCHA EQUATION
18861896
=============================
@@ -1926,12 +1936,14 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals
19261936
self-consistent equation.
19271937
"""
19281938

1929-
t1 = time.time()
19301939
super_struct = self.current_dyn.structure.generate_supercell(self.supercell)
19311940
#supercell_dyn = self.current_dyn.GenerateSupercellDyn(self.supercell)
19321941

19331942
# Dyagonalize
1934-
w, pols = self.current_dyn.DiagonalizeSupercell()#supercell_dyn.DyagDinQ(0)
1943+
if timer:
1944+
w, pols = timer.execute_timed_function(self.current_dyn.DiagonalizeSupercell)
1945+
else:
1946+
w, pols = self.current_dyn.DiagonalizeSupercell()#supercell_dyn.DyagDinQ(0)
19351947

19361948
if not self.ignore_small_w:
19371949
trans = CC.Methods.get_translations(pols, super_struct.get_masses_array())
@@ -1949,6 +1961,8 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals
19491961
nat = super_struct.N_atoms
19501962
eforces = np.zeros((self.N, nat, 3), dtype = np.float64, order = "F")
19511963
u_disp = np.zeros((self.N, nat, 3), dtype = np.float64, order = "F")
1964+
1965+
t1 = time.time()
19521966
#print nat
19531967
if subtract_sscha:
19541968
eforces[:,:,:] = self.forces - self.sscha_forces
@@ -1962,25 +1976,34 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals
19621976
pols = np.real(pols)
19631977

19641978
t2 = time.time()
1965-
if verbose:
1966-
print( " [GRADIENT] Time to prepare the gradient calculation:", t2 -t1,"s")
1979+
if timer:
1980+
timer.add_timer("Preparation of the gradient", t2 - t1)
19671981

19681982

1969-
t1 = time.time()
19701983
if fast_grad or not preconditioned:
1971-
grad, grad_err = SCHAModules.get_gradient_supercell(self.rho, u_disp, eforces, w, pols, trans,
1984+
if timer:
1985+
grad, grad_err = timer.execute_timed_function(SCHAModules.get_gradient_supercell,
1986+
self.rho, u_disp, eforces, w, pols, trans,
1987+
self.current_T, mass, ityp, log_err, self.N,
1988+
nat, 3*nat, len(mass), preconditioned,
1989+
override_name = "get_gradient_supercell")
1990+
else:
1991+
grad, grad_err = SCHAModules.get_gradient_supercell(self.rho, u_disp, eforces, w, pols, trans,
19721992
self.current_T, mass, ityp, log_err, self.N,
19731993
nat, 3*nat, len(mass), preconditioned)
19741994
else:
1975-
grad, grad_err = SCHAModules.get_gradient_supercell_new(self.rho, u_disp, eforces, w, pols, trans,
1995+
if timer:
1996+
grad, grad_err = timer.execute_timed_function(SCHAModules.get_gradient_supercell_new,
1997+
self.rho, u_disp, eforces, w, pols, trans,
1998+
self.current_T, mass, ityp, log_err, self.N,
1999+
nat, 3*nat, len(mass),
2000+
override_name = "get_gradient_supercell_new")
2001+
else:
2002+
grad, grad_err = SCHAModules.get_gradient_supercell_new(self.rho, u_disp, eforces, w, pols, trans,
19762003
self.current_T, mass, ityp, log_err, self.N,
19772004
nat, 3*nat, len(mass))
19782005

19792006

1980-
t2 = time.time()
1981-
if verbose:
1982-
print (" [GRADIENT] Time to call the fortran code:", t2 - t1, "s")
1983-
19842007
# If we are at gamma, we can skip this part
19852008
# Which makes the code faster
19862009
if np.prod(self.dyn_0.GetSupercell()) > 1:
@@ -1989,17 +2012,37 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals
19892012
if return_error:
19902013
# Check if a multiprocessing function can be exploited
19912014
if hasattr(CC.Phonons, 'GetDynQFromFCSupercell_parallel') and CC.Settings.GetNProc() > 1:
1992-
q_grad,q_grad_err = CC.Phonons.GetDynQFromFCSupercell_parallel(grad, np.array(self.current_dyn.q_tot),
2015+
if timer:
2016+
q_grad,q_grad_err = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell_parallel,
2017+
grad, np.array(self.current_dyn.q_tot),
2018+
self.current_dyn.structure, super_struct,fc2=grad_err)
2019+
else:
2020+
q_grad,q_grad_err = CC.Phonons.GetDynQFromFCSupercell_parallel(grad, np.array(self.current_dyn.q_tot),
19932021
self.current_dyn.structure, super_struct,fc2=grad_err)
19942022
else:
1995-
q_grad,q_grad_err = CC.Phonons.GetDynQFromFCSupercell(grad, np.array(self.current_dyn.q_tot),
2023+
if timer:
2024+
q_grad,q_grad_err = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell,
2025+
grad, np.array(self.current_dyn.q_tot),
2026+
self.current_dyn.structure, super_struct,fc2=grad_err)
2027+
else:
2028+
q_grad,q_grad_err = CC.Phonons.GetDynQFromFCSupercell(grad, np.array(self.current_dyn.q_tot),
19962029
self.current_dyn.structure, super_struct,fc2=grad_err)
19972030
else:
19982031
if hasattr(CC.Phonons, 'GetDynQFromFCSupercell_parallel'):
1999-
q_grad = CC.Phonons.GetDynQFromFCSupercell_parallel(grad, np.array(self.current_dyn.q_tot),
2032+
if timer:
2033+
q_grad = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell_parallel,
2034+
grad, np.array(self.current_dyn.q_tot),
2035+
self.current_dyn.structure, super_struct)
2036+
else:
2037+
q_grad = CC.Phonons.GetDynQFromFCSupercell_parallel(grad, np.array(self.current_dyn.q_tot),
20002038
self.current_dyn.structure, super_struct)
20012039
else:
2002-
q_grad = CC.Phonons.GetDynQFromFCSupercell(grad, np.array(self.current_dyn.q_tot),
2040+
if timer:
2041+
q_grad = timer.execute_timed_function(CC.Phonons.GetDynQFromFCSupercell,
2042+
grad, np.array(self.current_dyn.q_tot),
2043+
self.current_dyn.structure, super_struct)
2044+
else:
2045+
q_grad = CC.Phonons.GetDynQFromFCSupercell(grad, np.array(self.current_dyn.q_tot),
20032046
self.current_dyn.structure, super_struct)
20042047
#q_grad_err = CC.Phonons.GetDynQFromFCSupercell(grad_err, np.array(self.current_dyn.q_tot),
20052048
# self.current_dyn.structure, supercell_dyn.structure)
@@ -2010,11 +2053,6 @@ def get_preconditioned_gradient(self, subtract_sscha = True, return_error = Fals
20102053
q_grad[0, :, :] = grad
20112054
q_grad_err[0, :, :] = grad_err
20122055

2013-
t1 = time.time()
2014-
if verbose:
2015-
print (" [GRADIENT] Time to get back in fourier space:", t1 - t2, "s")
2016-
2017-
20182056
if return_error:
20192057
return q_grad, q_grad_err
20202058
else:

0 commit comments

Comments
 (0)