Skip to content

Commit d310125

Browse files
ManexAlkortaManexAlkorta
committed
Fix in V permutations when summing in get_indep_2fc(), and improved outlogs
1 parent b78fd8e commit d310125

3 files changed

Lines changed: 53 additions & 21 deletions

File tree

‎Modules/Ensemble.py‎

Lines changed: 15 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3846,7 +3846,7 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v
38463846
return dyn_hessian, d3* 2.0 # Ha to Ry
38473847
return dyn_hessian
38483848

3849-
def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_full_hessian = True, verbose = False):
3849+
def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, eps = 1e-6, alpha_mix=0.3, get_full_hessian = True, verbose = False):
38503850
"""
38513851
Dev function.
38523852
@@ -3865,6 +3865,15 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_ful
38653865
include_v4 : bool
38663866
If True we include the fourth order force constant matrix.
38673867
This requires a lot of memory
3868+
do_scf : bool
3869+
If True, the W matrix is self-consitently converged. Otherwise, just the first
3870+
correction of the fourth-order is considered (False often enough in practice).
3871+
eps : float
3872+
Precision to reach for the max value of |W|. It is later weighted by alpha-mix
3873+
in practice to ensure convergence. Defaults to 1e-6.
3874+
alpha_mix : float
3875+
Mixing parameter used in the scf loop of W, where
3876+
W(n+1)=W(n)(1-alpha_mix)+alpha_mix*W(n+1). Default to 0.3.
38683877
get_full_hessian : bool
38693878
If True the full hessian matrix is returned, if false, only the correction to
38703879
the SSCHA dynamical matrix is returned.
@@ -4039,10 +4048,10 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_ful
40394048
nat = self.current_dyn.structure.N_atoms
40404049
n_modes = self.current_dyn.structure.N_atoms*mod[0]*mod[1]*mod[2]*3
40414050

4042-
ref_3fc = SCHAModules.module_hess.get_ref3fc(nat, orbit3a, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc, rot_3fc, self.ur, self.upsilon, f, self.rho, log_err, self.s_inv_cart, self.irt, self.translations_irt, verbose)
4051+
ref_3fc = SCHAModules.module_hess.get_ref3fc(nat, orbit3a, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc, rot_3fc, self.ur, self.upsilon, f, self.rho, log_err, self.s_inv_cart, self.irt, self.translations_irt, True)
40434052

40444053
vs_red = np.empty([nrefq2,nat*3,nat*3,nat*3], dtype=np.complex128)
4045-
vs_red = SCHAModules.module_hess.get_ref_vsq(refq2,trs_l,rot_3fc,ref_3fc,mapping_triplet,verbose)
4054+
vs_red = SCHAModules.module_hess.get_ref_vsq(refq2,trs_l,rot_3fc,ref_3fc,mapping_triplet,True)
40464055
trs_gq, daq = SCHAModules.get_gq(trs_aq, trs_wq, transq, self.current_T)
40474056
indep_fc = SCHAModules.module_hess.get_indep2fc(vs_red, refq2, refq2o, norbitq2, orbit2a, n_indep_elem, indep_elem, rot_cart, mapping, map_uc, map_tr, T_list, trs_qlist, trs_gq, verbose)
40484057

@@ -4051,10 +4060,10 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_ful
40514060

40524061
orbit4t, orbit4o, norbit_4, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot_4fc, mapping_quadruplet = Classify.recognize_quadruplet(self.current_dyn, mapping, map_uc, verbose)
40534062

4054-
ref_4fc = SCHAModules.module_hess.get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot_4fc, self.ur, self.upsilon, f, self.rho, log_err, self.s_inv_cart, self.irt, self.translations_irt, verbose)
4063+
ref_4fc = SCHAModules.module_hess.get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot_4fc, self.ur, self.upsilon, f, self.rho, log_err, self.s_inv_cart, self.irt, self.translations_irt, True)
40554064

40564065
ws_red = np.zeros([nrefq4,nat*3,nat*3,nat*3,nat*3], dtype=np.complex128)
4057-
ws_red = SCHAModules.module_hess.get_ref_wsq(refq4,trs_l,rot_4fc,ref_4fc,mapping_quadruplet,verbose)
4066+
ws_red = SCHAModules.module_hess.get_ref_wsq(refq4,trs_l,rot_4fc,ref_4fc,mapping_quadruplet,True)
40584067

40594068
degs = qClassify.find_degeneracies(trs_wq)
40604069
Pmn = qClassify.construct_Pmn(mapping, orbitq1a, orbitq1s, trs_polvecs, rot_cart)
@@ -4063,7 +4072,7 @@ def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_ful
40634072
vs = SCHAModules.module_hess.get_all_vsq(trs_l, v_red, map_uc)
40644073

40654074
if do_scf:
4066-
ws_red_scf = SCHAModules.module_hess.get_scf_wsq(ws_red, trs_gq, refq4, refq4o, norbitq4, Pmn, degs, verbose)
4075+
ws_red_scf = SCHAModules.module_hess.get_scf_wsq(ws_red, trs_gq, refq4, refq4o, norbitq4, Pmn, degs, True, eps, alpha_mix)
40674076
indep_fc4 = SCHAModules.module_hess.get_indep2fc_v4(vs, ws_red_scf, refq4, refq4o, norbitq4, orbit2a, n_indep_elem, indep_elem, trs_gq, Pmn, degs, mapping, rot_cart, verbose)
40684077
else:
40694078
indep_fc4 = SCHAModules.module_hess.get_indep2fc_v4(vs, ws_red, refq4, refq4o, norbitq4, orbit2a, n_indep_elem, indep_elem, trs_gq, Pmn, degs, mapping, rot_cart, verbose)

‎Modules/qClassify.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ def map_singlet(q_list, q_list_frac, rcell, rot_cart):
5252
q_list_frac_fixed = np.empty(q_list_frac.shape)
5353
its_zb = np.empty(q_list_frac.shape[0], dtype=np.int8)
5454
for qi,q in enumerate(q_list):
55-
q_list_frac_fixed[qi] = CC.Methods.cart_to_cryst(rcell, _map_q_to_1st_bz(rcell, q))
55+
q_list_frac_fixed[qi] = np.round(CC.Methods.cart_to_cryst(rcell, _map_q_to_1st_bz(rcell, q)),6)
5656
if np.all((np.abs(np.round(2*q_list_frac_fixed[qi], 6))%1)<1e-6):
5757
its_zb[qi] = 0
5858
else:
@@ -70,7 +70,7 @@ def map_singlet(q_list, q_list_frac, rcell, rot_cart):
7070
for isym in range(rot_cart.shape[0]):
7171
q_sym_cart = rot_cart[isym] @ q
7272
q_sym_cart_1bz = _map_q_to_1st_bz(rcell, q_sym_cart)
73-
q_sym = CC.Methods.cart_to_cryst(rcell, q_sym_cart_1bz)
73+
q_sym = np.round(CC.Methods.cart_to_cryst(rcell, q_sym_cart_1bz),6)
7474
match = np.all(np.abs(q_list_frac_fixed-q_sym)<1e-3, axis=1)
7575
qii = np.where(match)
7676
mapping[qi,isym] = qii[0][0]

‎SCHAModules/module_hess.f90‎

Lines changed: 36 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -277,6 +277,7 @@ subroutine get_ref_vsq(refq2,pol_vecs,rot_3fc,ref_3fc,mapping_triplet,verbose,v1
277277
complex, dimension(nrefq2, n_mode, n_mode, n_mode), intent(out) :: v1
278278

279279
double precision, dimension(27) :: aux_3fc
280+
double precision :: tstart, tend
280281
integer :: n_mode_sc, n_mode, ns, qmu, qnu, mu, nu, at3, iq, index, nrefq2, dimq2, i,j, rq2
281282
integer :: nperm,nsym,nref3,ref3,equiv,iperm,isym,a,b,c,alpha,beta,gamma,nat,nat_sc
282283
complex, dimension(:), allocatable :: laux1, laux2
@@ -293,13 +294,17 @@ subroutine get_ref_vsq(refq2,pol_vecs,rot_3fc,ref_3fc,mapping_triplet,verbose,v1
293294
allocate(laux2(n_mode_sc))
294295

295296
if (verbose) then
296-
print*, "======================= get_ref_vsq() ======================="
297-
print*, "Computing all reference V_{a}^{\alpha}(q1,q2)."
297+
tstart = omp_get_wtime()
298+
print*, "======================= get_ref_vsq() ======================"
299+
print*, ""
300+
print*, " Computing all reference V_{a}^{\alpha}(q1,q2)."
301+
print*, ""
298302
end if
299303
v1 = 0
300304
!$omp parallel private (qmu,qnu,at3,a,alpha,mu,laux1,lres1,ref3,equiv,iperm,isym,i,j)
301305
!$omp do schedule (dynamic, 1) private (b,c,aux_3fc,beta,gamma,index,nu,laux2)
302306
do rq2 = 1, nrefq2
307+
print*, " Calculating V(rq2=", rq2,") out of", nrefq2
303308
qmu = refq2(rq2,1,1)
304309
qnu = refq2(rq2,1,2)
305310
do at3 = 1, n_mode
@@ -341,6 +346,9 @@ subroutine get_ref_vsq(refq2,pol_vecs,rot_3fc,ref_3fc,mapping_triplet,verbose,v1
341346
!$omp end do
342347
!$omp end parallel
343348
if (verbose) then
349+
tend = omp_get_wtime()
350+
print*, ""
351+
print*, "Elapsed time inside get_ref_vsq():", tend-tstart, "seconds."
344352
print*, "======================= DONE ======================="
345353
print*, ""
346354
endif
@@ -421,16 +429,20 @@ subroutine get_ref_wsq(refq4,pol_vecs,rot_4fc,ref_4fc,mapping_quadruplet,verbose
421429

422430
if (verbose) then
423431
tstart = omp_get_wtime()
424-
print*, "======================= get_ref_wsq() ======================="
425-
print*, "Computing the", nrefq4, "reference W(-q1,q2,q3,-q4)."
426-
print*, "This might take some time..."
432+
print*, "======================= get_ref_wsq() ======================"
433+
print*, ""
434+
print*, " Computing the", nrefq4, "reference W(-q1,q2,q3,-q4)."
435+
print*, ""
436+
print*, " This might take some time..."
437+
print*, ""
427438
end if
428439

429440
ns = n_mode_sc
430441

431442
v1 = 0.0d0
432443
tstart = omp_get_wtime()
433444
do rq4 = 1, nrefq4
445+
print*, " Calculating W(rq4=", rq4,") out of", nrefq4
434446
q1 = refq4(rq4,1,1)+1
435447
q2 = refq4(rq4,1,2)+1
436448
q3 = refq4(rq4,1,3)+1
@@ -520,13 +532,14 @@ subroutine get_ref_wsq(refq4,pol_vecs,rot_4fc,ref_4fc,mapping_quadruplet,verbose
520532
end do
521533
if (verbose) then
522534
tend = omp_get_wtime()
523-
print*, "Elapsed time inside get_ref_wsq():", tend-tstart
535+
print*, ""
536+
print*, "Elapsed time inside get_ref_wsq():", tend-tstart, "seconds."
524537
print*, "======================= DONE ======================="
525538
print*, ""
526539
endif
527540
end subroutine get_ref_wsq
528541

529-
subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
542+
subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps, alpha_mix, &
530543
wsq_scf, nrefq4, n_mode, iq, dimq4, nsym)
531544

532545
implicit none
@@ -543,6 +556,8 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
543556
logical, dimension(iq, n_mode, n_mode), intent(in) :: degs
544557
logical, intent(in) :: verbose
545558

559+
double precision, intent(in) :: eps, alpha_mix
560+
546561
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode), intent(out) :: wsq_scf
547562

548563
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode) :: wsq2
@@ -560,8 +575,6 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
560575
integer :: q10,q20,q30,q40,q5,q6,iperm_1,iperm_2,isym_1,isym_2
561576

562577
complex(8) :: ktea1, ktea2, W1, W2
563-
double precision, parameter :: eps = 1e-6
564-
double precision, parameter :: alpha_mix = 0.5d0
565578
integer, parameter :: maxiter = 50
566579

567580
perms = reshape([ &
@@ -575,12 +588,16 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
575588
if (verbose) then
576589
tstart = omp_get_wtime()
577590
print*, "======================= get_scf_wsq() ======================="
578-
print*, "Self-consintent loop to compute \Theta(-q1,q2,q3,-q4)."
591+
print*, ""
592+
print*, " Self-consintent loop to compute \Theta(-q1,q2,q3,-q4)."
593+
print*, ""
594+
print*, " eps=", eps, " alpha_mix=", alpha_mix
579595
print*, ""
580596
end if
581597
do1 : do iter = 1, maxiter
582598
if (verbose) then
583599
print*, " Iteration", iter
600+
print*, " ---------------------"
584601
end if
585602
if (iter == 1) then
586603
wsq2 = wsq1
@@ -779,7 +796,7 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
779796
tend = omp_get_wtime()
780797
print*, "Convergence found with", iter, "iterations"
781798
print*, ""
782-
print*, "Elapsed time inside get_scf_wsq():", tend-tstart
799+
print*, "Elapsed time inside get_scf_wsq():", tend-tstart, "seconds."
783800
print*, "======================= DONE ======================="
784801
print*, ""
785802
endif
@@ -839,6 +856,7 @@ subroutine get_sum_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
839856
do1 : do iter = 1, maxiter
840857
if (verbose) then
841858
print*, " Iteration", iter
859+
print*, " -----------------------"
842860
end if
843861
if (iter == 1) then
844862
wsq2 = wsq1
@@ -1041,7 +1059,7 @@ subroutine get_sum_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
10411059
tend = omp_get_wtime()
10421060
print*, "Convergence found with", iter, "iterations"
10431061
print*, ""
1044-
print*, "Elapsed time inside get_scf_wsq():", tend-tstart
1062+
print*, "Elapsed time inside get_scf_wsq():", tend-tstart, "seconds."
10451063
print*, "======================= DONE ======================="
10461064
print*, ""
10471065
endif
@@ -1119,8 +1137,13 @@ subroutine get_indep2fc( &
11191137
Vsb = &
11201138
dot_product(rot_cart(isym+1,beta+1,:),vs_q(nat2p_uc*3+1:(nat2p_uc+1)*3+1)) * &
11211139
exp(-2*j*pi*dot_product(q_list(refq2(rq2,1,1)+1,:)-q_list(refq2(rq2,1,2)+1,:),T_list(map_tr(nat2p+1)+1,:)))
1122-
indep_fc(ref2,i) = indep_fc(ref2,i) + &
1140+
if (iperm == 0) then
1141+
indep_fc(ref2,i) = indep_fc(ref2,i) + &
11231142
real(0.5d0*F(refq2(rq2,1,1)+1,refq2(rq2,1,2)+1,munu(permutations(iperm+1,1)), munu(permutations(iperm+1,2)))*Vsa*CONJG(Vsb))
1143+
else
1144+
indep_fc(ref2,i) = indep_fc(ref2,i) + &
1145+
real(0.5d0*F(refq2(rq2,1,1)+1,refq2(rq2,1,2)+1,munu(permutations(iperm+1,1)), munu(permutations(iperm+1,2)))*CONJG(Vsa)*Vsb)
1146+
end if
11241147
end do
11251148
end do
11261149
end do

0 commit comments

Comments
 (0)