Skip to content

Commit b78fd8e

Browse files
ManexAlkortaManexAlkorta
committed
Differentiate between get_scf_wsq and get_sum_wsq
1 parent 9701230 commit b78fd8e

2 files changed

Lines changed: 273 additions & 13 deletions

File tree

‎Modules/Ensemble.py‎

Lines changed: 7 additions & 5 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, get_full_hessian = True, verbose = False):
3849+
def get_free_energy_hessian_dev(self, include_v4 = False, do_scf = True, get_full_hessian = True, verbose = False):
38503850
"""
38513851
Dev function.
38523852
@@ -4059,12 +4059,14 @@ def get_free_energy_hessian_dev(self, include_v4 = False, get_full_hessian = Tru
40594059
degs = qClassify.find_degeneracies(trs_wq)
40604060
Pmn = qClassify.construct_Pmn(mapping, orbitq1a, orbitq1s, trs_polvecs, rot_cart)
40614061

4062-
ws_red_scf = SCHAModules.module_hess.get_scf_wsq(ws_red, trs_gq, refq4, refq4o, norbitq4, Pmn, degs, verbose)
4063-
40644062
v_red, ref_3fc = SCHAModules.module_hess.get_v3_red(nat, norbit3, orbit3a, orbit3s, 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)
40654063
vs = SCHAModules.module_hess.get_all_vsq(trs_l, v_red, map_uc)
4066-
4067-
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)
4064+
4065+
if do_scf:
4066+
ws_red_scf = SCHAModules.module_hess.get_scf_wsq(ws_red, trs_gq, refq4, refq4o, norbitq4, Pmn, degs, verbose)
4067+
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)
4068+
else:
4069+
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)
40684070
indep_fc += indep_fc4
40694071
phi_sc_odd = np.zeros((n_modes, n_modes), dtype = np.double)
40704072
for ref2 in range(nref2):

‎SCHAModules/module_hess.f90‎

Lines changed: 266 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -160,12 +160,6 @@ subroutine get_ref3fc(nat, orbit3t, indep_3fc_elem, n_indep_3fc_elem, kernel_3fc
160160
!$omp parallel private (i, nat1, nat2, nat3, alpha, beta, gamma, v3, indep_3fc, tmp_ref_3fc)
161161
!$omp do schedule (dynamic, 1) private (v3, indep_3fc, tmp_ref_3fc)
162162
do ref3 = 1, nref3
163-
!debugging parallel installation
164-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166-
print *, "Iteration ", ref3, "executed by thread ", omp_get_thread_num()
167-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
168-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
169163
nat1 = orbit3t(ref3,1,1)
170164
nat2 = orbit3t(ref3,1,2)
171165
nat3 = orbit3t(ref3,1,3)
@@ -567,6 +561,7 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
567561

568562
complex(8) :: ktea1, ktea2, W1, W2
569563
double precision, parameter :: eps = 1e-6
564+
double precision, parameter :: alpha_mix = 0.5d0
570565
integer, parameter :: maxiter = 50
571566

572567
perms = reshape([ &
@@ -575,11 +570,270 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
575570
4,2,3,1, 2,4,1,3, 3,1,4,2, 1,3,2,4], &
576571
[12,4], order=[2,1])
577572

573+
wsq_scf = wsq1
574+
578575
if (verbose) then
579576
tstart = omp_get_wtime()
580577
print*, "======================= get_scf_wsq() ======================="
581578
print*, "Self-consintent loop to compute \Theta(-q1,q2,q3,-q4)."
582579
print*, ""
580+
end if
581+
do1 : do iter = 1, maxiter
582+
if (verbose) then
583+
print*, " Iteration", iter
584+
end if
585+
if (iter == 1) then
586+
wsq2 = wsq1
587+
else
588+
wsq2 = wsq_scf
589+
end if
590+
wsqn = (0.0d0,0.0d0)
591+
do rq4_0 = 1, nrefq4
592+
q10 = refq4(rq4_0,1,1)+1
593+
q20 = refq4(rq4_0,1,2)+1
594+
q30 = refq4(rq4_0,1,3)+1
595+
q40 = refq4(rq4_0,1,4)+1
596+
do rq4_1 = 1, nrefq4
597+
q1_ref = refq4(rq4_1,1,1)+1
598+
q2_ref = refq4(rq4_1,1,2)+1
599+
q5_1ref = refq4(rq4_1,1,3)+1
600+
q6_1ref = refq4(rq4_1,1,4)+1
601+
do q4i_1 = 1, norbitq4(rq4_1)
602+
q1 = refq4(rq4_1,q4i_1,1)+1
603+
q2 = refq4(rq4_1,q4i_1,2)+1
604+
q5_1 = refq4(rq4_1,q4i_1,3)+1
605+
q6_1 = refq4(rq4_1,q4i_1,4)+1
606+
iperm_1 = refq4o(rq4_1,q4i_1,1)
607+
isym_1 = refq4o(rq4_1,q4i_1,2)+1
608+
if ((q10 == q1) .and. (q20 == q2)) then
609+
do rq4_2 = 1, nrefq4
610+
q5_2ref = refq4(rq4_2,1,1)+1
611+
q6_2ref = refq4(rq4_2,1,2)+1
612+
q3_ref = refq4(rq4_2,1,3)+1
613+
q4_ref = refq4(rq4_2,1,4)+1
614+
do q4i_2 = 1, norbitq4(rq4_2)
615+
q5_2 = refq4(rq4_2,q4i_2,1)+1
616+
q6_2 = refq4(rq4_2,q4i_2,2)+1
617+
q3 = refq4(rq4_2,q4i_2,3)+1
618+
q4 = refq4(rq4_2,q4i_2,4)+1
619+
iperm_2 = refq4o(rq4_2,q4i_2,1)
620+
isym_2 = refq4o(rq4_2,q4i_2,2)+1
621+
if ((q30 == q3) .and. (q40 == q4)) then
622+
if ((q5_1 == q5_2) .and. (q6_1 == q6_2)) then
623+
!$omp parallel &
624+
!$omp private (mu2,mu3,mu4,mu5,mu6) &
625+
!$omp private (nu1,nu2,nu5,nu6,lambda5,lambda6,lambda3,lambda4,nunu,lala,mumu1,mumu2)
626+
!$omp do schedule (dynamic, 1) private (ktea1,ktea2,W1,W2)
627+
do mu1 = 1, n_mode
628+
do nu1 = 1, n_mode
629+
if (degs(q1,nu1,mu1)) then
630+
do mu2 = 1, n_mode
631+
do nu2 = 1, n_mode
632+
if (degs(q2,nu2,mu2)) then
633+
do mu5 = 1, n_mode
634+
do nu5 = 1, n_mode
635+
if (degs(q5_1,nu5,mu5)) then
636+
do mu6 = 1, n_mode
637+
mumu1 = [mu1,mu2,mu5,mu6]
638+
do nu6 = 1, n_mode
639+
nunu = [nu1,nu2,nu5,nu6]
640+
if (degs(q6_1,nu6,mu6)) then
641+
ktea1 = &
642+
P(q1_ref,mumu1(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,1)),isym_1)*&
643+
CONJG(P(q2_ref,mumu1(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,2)),isym_1))*&
644+
CONJG(P(q5_1ref,mumu1(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,3)),isym_1))*&
645+
P(q6_1ref,mumu1(perms(iperm_1+1,4)),nunu(perms(iperm_1+1,4)),isym_1)
646+
if (iperm_1 == 0) then
647+
W1 = &
648+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
649+
elseif (iperm_1 == 1) then
650+
W1 = &
651+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
652+
elseif (iperm_1 == 2) then
653+
W1 = &
654+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
655+
elseif (iperm_1 == 3) then
656+
W1 = &
657+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
658+
elseif (iperm_1 == 4) then
659+
W1 = &
660+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
661+
elseif (iperm_1 == 5) then
662+
W1 = &
663+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
664+
elseif (iperm_1 == 6) then
665+
W1 = &
666+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
667+
elseif (iperm_1 == 7) then
668+
W1 = &
669+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
670+
elseif (iperm_1 == 8) then
671+
W1 = &
672+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
673+
elseif (iperm_1 == 9) then
674+
W1 = &
675+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
676+
elseif (iperm_1 == 10) then
677+
W1 = &
678+
CONJG(ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4))))
679+
elseif (iperm_1 == 11) then
680+
W1 = &
681+
ktea1*wsq2(rq4_1,nunu(perms(iperm_1+1,1)),nunu(perms(iperm_1+1,2)),nunu(perms(iperm_1+1,3)),nunu(perms(iperm_1+1,4)))
682+
end if
683+
do mu3 = 1, n_mode
684+
do lambda3 = 1, n_mode
685+
if (degs(q3,lambda3,mu3)) then
686+
do mu4 = 1, n_mode
687+
do lambda4 = 1, n_mode
688+
mumu2 = [mu5,mu6,mu3,mu4]
689+
if (degs(q4,lambda4,mu4)) then
690+
do lambda5 = 1, n_mode
691+
if (degs(q5_2,lambda5,mu5)) then
692+
do lambda6 = 1, n_mode
693+
lala = [lambda5,lambda6,lambda3,lambda4]
694+
if (degs(q6_2,lambda6,mu6)) then
695+
ktea2 = &
696+
P(q5_2ref,mumu2(perms(iperm_2+1,1)),lala(perms(iperm_2+1,1)),isym_2)*&
697+
CONJG(P(q6_2ref,mumu2(perms(iperm_2+1,2)),lala(perms(iperm_2+1,2)),isym_2))*&
698+
CONJG(P(q3_ref,mumu2(perms(iperm_2+1,3)),lala(perms(iperm_2+1,3)),isym_2))*&
699+
P(q4_ref,mumu2(perms(iperm_2+1,4)),lala(perms(iperm_2+1,4)),isym_2)
700+
if (iperm_2 == 0) then
701+
W2 = &
702+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
703+
elseif (iperm_2 == 1) then
704+
W2 = &
705+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
706+
elseif (iperm_2 == 2) then
707+
W2 = &
708+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
709+
elseif (iperm_2 == 3) then
710+
W2 = &
711+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
712+
elseif (iperm_2 == 4) then
713+
W2 = &
714+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
715+
elseif (iperm_2 == 5) then
716+
W2 = &
717+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
718+
elseif (iperm_2 == 6) then
719+
W2 = &
720+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
721+
elseif (iperm_2 == 7) then
722+
W2 = &
723+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
724+
elseif (iperm_2 == 8) then
725+
W2 = &
726+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
727+
elseif (iperm_2 == 9) then
728+
W2 = &
729+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
730+
elseif (iperm_2 == 10) then
731+
W2 = &
732+
CONJG(ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4))))
733+
elseif (iperm_2 == 11) then
734+
W2 = &
735+
ktea2*wsq1(rq4_2,lala(perms(iperm_2+1,1)),lala(perms(iperm_2+1,2)),lala(perms(iperm_2+1,3)),lala(perms(iperm_2+1,4)))
736+
end if
737+
wsqn(rq4_0,mu1,mu2,mu3,mu4) = &
738+
wsqn(rq4_0,mu1,mu2,mu3,mu4) + 0.5d0*F(q5_1,q6_1,mu5,mu6)*W1*W2
739+
end if
740+
end do
741+
end if
742+
end do
743+
end if
744+
end do
745+
end do
746+
end if
747+
end do
748+
end do
749+
end if
750+
end do
751+
end do
752+
end if
753+
end do
754+
end do
755+
end if
756+
end do
757+
end do
758+
end if
759+
end do
760+
end do
761+
!$omp end do
762+
!$omp end parallel
763+
end if
764+
end if
765+
end do
766+
end do
767+
end if
768+
end do
769+
end do
770+
end do
771+
wsq_scf = wsq_scf*(1-alpha_mix) +alpha_mix*(wsq1+wsqn)
772+
if (verbose) then
773+
print*, " MAXVAL(|W(n)-W(n-1)|)", MAXVAL(ABS(wsq_scf-wsq2)), ". Threshold (weighted by alpha_mix)=", eps*alpha_mix
774+
print*, ""
775+
end if
776+
if (MAXVAL(ABS(wsq_scf-wsq2)) < eps*alpha_mix) exit do1
777+
end do do1
778+
if (verbose) then
779+
tend = omp_get_wtime()
780+
print*, "Convergence found with", iter, "iterations"
781+
print*, ""
782+
print*, "Elapsed time inside get_scf_wsq():", tend-tstart
783+
print*, "======================= DONE ======================="
784+
print*, ""
785+
endif
786+
787+
end subroutine get_scf_wsq
788+
789+
subroutine get_sum_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
790+
wsq_scf, nrefq4, n_mode, iq, dimq4, nsym)
791+
792+
implicit none
793+
794+
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode), intent(in) :: wsq1
795+
796+
double precision, dimension(iq, iq, n_mode, n_mode), intent(in) :: F
797+
798+
integer, dimension(nrefq4,dimq4,4), intent(in) :: refq4
799+
integer, dimension(nrefq4,dimq4,2), intent(in) :: refq4o
800+
integer, dimension(nrefq4), intent(in) :: norbitq4
801+
802+
complex(8), dimension(iq, n_mode, n_mode, nsym), intent(in) :: P
803+
logical, dimension(iq, n_mode, n_mode), intent(in) :: degs
804+
logical, intent(in) :: verbose
805+
806+
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode), intent(out) :: wsq_scf
807+
808+
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode) :: wsq2
809+
complex(8), dimension(nrefq4,n_mode,n_mode,n_mode,n_mode) :: wsqn
810+
811+
double precision :: tstart, tend
812+
813+
integer, dimension(12,4) :: perms
814+
integer, dimension(4) :: nunu, lala, mumu1, mumu2
815+
integer :: nrefq4, n_mode, iq, dimq4, nsym
816+
integer :: mu1,mu2,mu3,mu4,mu5,mu6,iter
817+
integer :: rq4_0,rq4_1,rq4_2,q4i_1,q4i_2,q1,q2,q3,q4,q5_1,q6_1,q5_2,q6_2
818+
integer :: q1_ref, q2_ref, q3_ref, q4_ref, q5_1ref, q5_2ref, q6_1ref, q6_2ref
819+
integer :: nu1,nu2,nu5,nu6,lambda3,lambda4,lambda5,lambda6
820+
integer :: q10,q20,q30,q40,q5,q6,iperm_1,iperm_2,isym_1,isym_2
821+
822+
complex(8) :: ktea1, ktea2, W1, W2
823+
double precision, parameter :: eps = 1e-6
824+
integer, parameter :: maxiter = 50
825+
826+
perms = reshape([ &
827+
1,2,3,4, 2,1,4,3, 3,4,1,2, 4,3,2,1, &
828+
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1, &
829+
4,2,3,1, 2,4,1,3, 3,1,4,2, 1,3,2,4], &
830+
[12,4], order=[2,1])
831+
832+
if (verbose) then
833+
tstart = omp_get_wtime()
834+
print*, "======================= get_sum_wsq() ======================="
835+
print*, "Looking where to truncate the analytical sum of \Theta(-q1,q2,q3,-q4)."
836+
print*, ""
583837
end if
584838
wsq_scf = wsq1
585839
do1 : do iter = 1, maxiter
@@ -773,11 +1027,15 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
7731027
end do
7741028
end do
7751029
if (verbose) then
776-
print*, " MAXVAL(W^{n})=", MAXVAL(ABS(wsqn)), ". Threshold=", eps
1030+
print*, " MAXVAL(|W^{n}|)=", MAXVAL(ABS(wsqn)), ". Threshold=", eps
7771031
print*, ""
7781032
end if
7791033
wsq_scf = wsq_scf + wsqn
7801034
if (MAXVAL(ABS(wsqn)) < eps) exit do1
1035+
if (MAXVAL(ABS(wsqn)) > MAXVAL(ABS(wsq2))) then
1036+
print*, "W^{(n)} is not decaying. Returning first order correction..."
1037+
wsqn=wsq1
1038+
end if
7811039
end do do1
7821040
if (verbose) then
7831041
tend = omp_get_wtime()
@@ -788,7 +1046,7 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
7881046
print*, ""
7891047
endif
7901048

791-
end subroutine get_scf_wsq
1049+
end subroutine get_sum_wsq
7921050

7931051
subroutine get_indep2fc( &
7941052
vs_red, refq2, refq2o, norbitq2, orbit2t, n_indep_elem, indep_elem, rot_cart, mapping, map_uc, map_tr, T_list, q_list, F, &

0 commit comments

Comments
 (0)