Skip to content

Commit 8e97e02

Browse files
ManexAlkortaManexAlkorta
committed
Minimal changes in the logout
1 parent d310125 commit 8e97e02

1 file changed

Lines changed: 44 additions & 91 deletions

File tree

‎SCHAModules/module_hess.f90‎

Lines changed: 44 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -225,7 +225,7 @@ subroutine get_ref4fc(orbit4t, indep_4fc_elem, n_indep_4fc_elem, kernel_4fc, rot
225225

226226
if (verbose) then
227227
print *, "======================= get_ref4fc() ======================= "
228-
print *, "Computing the independent elements for the 4th order FCs."
228+
print *, " Computing the independent elements for the 4th order FCs"
229229
end if
230230

231231
!$omp parallel private (i, nat1, nat2, nat3, nat4, alpha, beta, gamma, delta)
@@ -297,14 +297,15 @@ subroutine get_ref_vsq(refq2,pol_vecs,rot_3fc,ref_3fc,mapping_triplet,verbose,v1
297297
tstart = omp_get_wtime()
298298
print*, "======================= get_ref_vsq() ======================"
299299
print*, ""
300-
print*, " Computing all reference V_{a}^{\alpha}(q1,q2)."
300+
print*, " Computing all reference V_{a}^{\alpha}(q1,q2)"
301301
print*, ""
302302
end if
303303
v1 = 0
304304
!$omp parallel private (qmu,qnu,at3,a,alpha,mu,laux1,lres1,ref3,equiv,iperm,isym,i,j)
305305
!$omp do schedule (dynamic, 1) private (b,c,aux_3fc,beta,gamma,index,nu,laux2)
306306
do rq2 = 1, nrefq2
307-
print*, " Calculating V(rq2=", rq2,") out of", nrefq2
307+
print "(A, I8, A, I8)", &
308+
" Calculating V(rq2=", rq2,") out of", nrefq2
308309
qmu = refq2(rq2,1,1)
309310
qnu = refq2(rq2,1,2)
310311
do at3 = 1, n_mode
@@ -348,7 +349,7 @@ subroutine get_ref_vsq(refq2,pol_vecs,rot_3fc,ref_3fc,mapping_triplet,verbose,v1
348349
if (verbose) then
349350
tend = omp_get_wtime()
350351
print*, ""
351-
print*, "Elapsed time inside get_ref_vsq():", tend-tstart, "seconds."
352+
print "(A, ES16.6,A)", " Elapsed time inside get_ref_vsq():", tend-tstart, " seconds"
352353
print*, "======================= DONE ======================="
353354
print*, ""
354355
endif
@@ -431,7 +432,8 @@ subroutine get_ref_wsq(refq4,pol_vecs,rot_4fc,ref_4fc,mapping_quadruplet,verbose
431432
tstart = omp_get_wtime()
432433
print*, "======================= get_ref_wsq() ======================"
433434
print*, ""
434-
print*, " Computing the", nrefq4, "reference W(-q1,q2,q3,-q4)."
435+
print "(A, I8, A)", &
436+
" Computing the", nrefq4, " reference W^{(1)}(q1,q2,q3,q4)"
435437
print*, ""
436438
print*, " This might take some time..."
437439
print*, ""
@@ -442,7 +444,7 @@ subroutine get_ref_wsq(refq4,pol_vecs,rot_4fc,ref_4fc,mapping_quadruplet,verbose
442444
v1 = 0.0d0
443445
tstart = omp_get_wtime()
444446
do rq4 = 1, nrefq4
445-
print*, " Calculating W(rq4=", rq4,") out of", nrefq4
447+
print "(A, I8, A, I8)", " Calculating W(rq4=", rq4,") out of", nrefq4
446448
q1 = refq4(rq4,1,1)+1
447449
q2 = refq4(rq4,1,2)+1
448450
q3 = refq4(rq4,1,3)+1
@@ -533,7 +535,7 @@ subroutine get_ref_wsq(refq4,pol_vecs,rot_4fc,ref_4fc,mapping_quadruplet,verbose
533535
if (verbose) then
534536
tend = omp_get_wtime()
535537
print*, ""
536-
print*, "Elapsed time inside get_ref_wsq():", tend-tstart, "seconds."
538+
print "(A, ES16.6,A)", " Elapsed time inside get_ref_wsq():", tend-tstart, " seconds"
537539
print*, "======================= DONE ======================="
538540
print*, ""
539541
endif
@@ -565,7 +567,7 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps,
565567

566568
double precision :: tstart, tend
567569

568-
integer, dimension(12,4) :: perms
570+
integer, dimension(8,4) :: perms
569571
integer, dimension(4) :: nunu, lala, mumu1, mumu2
570572
integer :: nrefq4, n_mode, iq, dimq4, nsym
571573
integer :: mu1,mu2,mu3,mu4,mu5,mu6,iter
@@ -574,29 +576,31 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps,
574576
integer :: nu1,nu2,nu5,nu6,lambda3,lambda4,lambda5,lambda6
575577
integer :: q10,q20,q30,q40,q5,q6,iperm_1,iperm_2,isym_1,isym_2
576578

579+
logical :: its_conv
577580
complex(8) :: ktea1, ktea2, W1, W2
578-
integer, parameter :: maxiter = 50
581+
integer, parameter :: maxiter = 10
579582

580583
perms = reshape([ &
581584
1,2,3,4, 2,1,4,3, 3,4,1,2, 4,3,2,1, &
582-
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1, &
583-
4,2,3,1, 2,4,1,3, 3,1,4,2, 1,3,2,4], &
584-
[12,4], order=[2,1])
585+
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1], &
586+
[8,4], order=[2,1])
585587

586588
wsq_scf = wsq1
587589

590+
its_conv = .true.
588591
if (verbose) then
589592
tstart = omp_get_wtime()
590593
print*, "======================= get_scf_wsq() ======================="
591594
print*, ""
592-
print*, " Self-consintent loop to compute \Theta(-q1,q2,q3,-q4)."
595+
print*, " Self-consintent loop to compute W(q1,q2,q3,q4)"
593596
print*, ""
594-
print*, " eps=", eps, " alpha_mix=", alpha_mix
597+
print "(A, ES16.6, A, ES16.6)", &
598+
" eps=", eps, "; alpha_mix=", alpha_mix
595599
print*, ""
596600
end if
597601
do1 : do iter = 1, maxiter
598602
if (verbose) then
599-
print*, " Iteration", iter
603+
print "(A, I4)", " Iteration", iter
600604
print*, " ---------------------"
601605
end if
602606
if (iter == 1) then
@@ -683,18 +687,6 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps,
683687
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))))
684688
elseif (iperm_1 == 7) then
685689
W1 = &
686-
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)))
687-
elseif (iperm_1 == 8) then
688-
W1 = &
689-
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)))
690-
elseif (iperm_1 == 9) then
691-
W1 = &
692-
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))))
693-
elseif (iperm_1 == 10) then
694-
W1 = &
695-
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))))
696-
elseif (iperm_1 == 11) then
697-
W1 = &
698690
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)))
699691
end if
700692
do mu3 = 1, n_mode
@@ -737,18 +729,6 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps,
737729
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))))
738730
elseif (iperm_2 == 7) then
739731
W2 = &
740-
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)))
741-
elseif (iperm_2 == 8) then
742-
W2 = &
743-
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)))
744-
elseif (iperm_2 == 9) then
745-
W2 = &
746-
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))))
747-
elseif (iperm_2 == 10) then
748-
W2 = &
749-
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))))
750-
elseif (iperm_2 == 11) then
751-
W2 = &
752732
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)))
753733
end if
754734
wsqn(rq4_0,mu1,mu2,mu3,mu4) = &
@@ -787,16 +767,28 @@ subroutine get_scf_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, eps,
787767
end do
788768
wsq_scf = wsq_scf*(1-alpha_mix) +alpha_mix*(wsq1+wsqn)
789769
if (verbose) then
790-
print*, " MAXVAL(|W(n)-W(n-1)|)", MAXVAL(ABS(wsq_scf-wsq2)), ". Threshold (weighted by alpha_mix)=", eps*alpha_mix
770+
print "(A, ES16.6)", &
771+
" MAXVAL(|W(n)-W(n-1)|) =", MAXVAL(ABS(wsq_scf-wsq2))
772+
print "(A, ES16.6)", &
773+
" Threshold (weighted by alpha_mix) =", eps*alpha_mix
791774
print*, ""
792775
end if
793776
if (MAXVAL(ABS(wsq_scf-wsq2)) < eps*alpha_mix) exit do1
777+
if (iter == maxiter) then
778+
wsq_scf = wsq1
779+
its_conv = .false.
780+
end if
794781
end do do1
795782
if (verbose) then
796783
tend = omp_get_wtime()
797-
print*, "Convergence found with", iter, "iterations"
784+
if (its_conv) then
785+
print "(A, I4, A)", " Convergence found with", iter, " iterations"
786+
else
787+
print "(A, I4, A)", " Convergence not found with", iter, " iterations"
788+
print "(A)", " Returning first order W^{(1)}(q1,q2,q3,q4)"
789+
end if
798790
print*, ""
799-
print*, "Elapsed time inside get_scf_wsq():", tend-tstart, "seconds."
791+
print "(A, ES16.6,A)", " Elapsed time inside get_scf_wsq():", tend-tstart, " seconds"
800792
print*, "======================= DONE ======================="
801793
print*, ""
802794
endif
@@ -1059,7 +1051,7 @@ subroutine get_sum_wsq(wsq1, F, refq4, refq4o, norbitq4, P, degs, verbose, &
10591051
tend = omp_get_wtime()
10601052
print*, "Convergence found with", iter, "iterations"
10611053
print*, ""
1062-
print*, "Elapsed time inside get_scf_wsq():", tend-tstart, "seconds."
1054+
print "(A, ES16.6,A)", " Elapsed time inside get_scf_wsq():", tend-tstart, " seconds"
10631055
print*, "======================= DONE ======================="
10641056
print*, ""
10651057
endif
@@ -1299,43 +1291,6 @@ subroutine get_indep2fc_v4( &
12991291
F(q10,q20,mu1,mu2) * &
13001292
F(q30,q40,mu3,mu4) * &
13011293
vsa*W*CONJG(vsb))
1302-
elseif (iperm == 8) then
1303-
vsa= &
1304-
dot_product(rot_cart(isym+1,alpha,:),vs(3*nat1p+1:3*(nat1p+1)+1,q40,q20,mu4,mu2))
1305-
vsb= &
1306-
dot_product(rot_cart(isym+1,beta,:),vs(3*nat2p+1:3*(nat2p+1)+1,q30,q10,mu3,mu1))
1307-
indep_fc(ref2,i) = indep_fc(ref2,i) + real(0.25d0* &
1308-
F(q10,q20,mu1,mu2) * &
1309-
F(q30,q40,mu3,mu4) * &
1310-
vsa*W*CONJG(vsb))
1311-
elseif (iperm == 9) then
1312-
vsa= &
1313-
dot_product(rot_cart(isym+1,alpha,:),vs(3*nat1p+1:3*(nat1p+1)+1,q20,q40,mu2,mu4))
1314-
vsb= &
1315-
dot_product(rot_cart(isym+1,beta,:),vs(3*nat2p+1:3*(nat2p+1)+1,q10,q30,mu1,mu3))
1316-
indep_fc(ref2,i) = indep_fc(ref2,i) + real(0.25d0* &
1317-
F(q10,q20,mu1,mu2) * &
1318-
F(q30,q40,mu3,mu4) * &
1319-
vsa*W*CONJG(vsb))
1320-
elseif (iperm == 10) then
1321-
vsa= &
1322-
dot_product(rot_cart(isym+1,alpha,:),vs(3*nat1p+1:3*(nat1p+1)+1,q30,q10,mu3,mu1))
1323-
vsb= &
1324-
dot_product(rot_cart(isym+1,beta,:),vs(3*nat2p+1:3*(nat2p+1)+1,q40,q20,mu4,mu2))
1325-
indep_fc(ref2,i) = indep_fc(ref2,i) + real(0.25d0* &
1326-
F(q10,q20,mu1,mu2) * &
1327-
F(q30,q40,mu3,mu4) * &
1328-
vsa*W*CONJG(vsb))
1329-
elseif (iperm == 11) then
1330-
vsa= &
1331-
dot_product(rot_cart(isym+1,alpha,:),vs(3*nat1p+1:3*(nat1p+1)+1,q10,q30,mu1,mu3))
1332-
vsb= &
1333-
dot_product(rot_cart(isym+1,beta,:),vs(3*nat2p+1:3*(nat2p+1)+1,q20,q40,mu2,mu4))
1334-
indep_fc(ref2,i) = indep_fc(ref2,i) + real(0.25d0* &
1335-
F(q10,q20,mu1,mu2) * &
1336-
F(q30,q40,mu3,mu4) * &
1337-
vsa*W*CONJG(vsb))
1338-
13391294
end if
13401295
end do
13411296
end do
@@ -2088,7 +2043,7 @@ subroutine get_q_nref4(mappings, nref4, norbitq4, iq, nsym)
20882043
integer, allocatable, dimension(:) :: norbit4
20892044
integer, dimension(12*nsym, 4) :: equilist
20902045

2091-
integer, dimension(12,4) :: permutations
2046+
integer, dimension(8,4) :: permutations
20922047
integer, dimension(4) :: qplet, qplet_perm, qplet_sym
20932048
integer :: q1, q2, q3, q4, nall4, equiv, iperm, isym, iq, nsym
20942049
logical :: its_in_list
@@ -2098,9 +2053,8 @@ subroutine get_q_nref4(mappings, nref4, norbitq4, iq, nsym)
20982053

20992054
permutations = reshape([ &
21002055
1,2,3,4, 2,1,4,3, 3,4,1,2, 4,3,2,1, &
2101-
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1, &
2102-
4,2,3,1, 2,4,1,3, 3,1,4,2, 1,3,2,4], &
2103-
[12,4], order=[2,1])
2056+
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1], &
2057+
[8,4], order=[2,1])
21042058

21052059
!permutations = reshape([1,2,3,4], [1,4])
21062060

@@ -2116,7 +2070,7 @@ subroutine get_q_nref4(mappings, nref4, norbitq4, iq, nsym)
21162070
if (its_in_list) cycle do4
21172071
nref4 = nref4 + 1
21182072
equiv = 0
2119-
do iperm = 1, 12
2073+
do iperm = 1, 8
21202074
qplet_perm(1) = qplet(permutations(iperm,1))
21212075
qplet_perm(2) = qplet(permutations(iperm,2))
21222076
qplet_perm(3) = qplet(permutations(iperm,3))
@@ -2158,7 +2112,7 @@ subroutine recognize_q_quadruplet(nref4, norbitq4, q_list, mappings, verbose, &
21582112
integer, dimension(iq**4,4) :: all4
21592113
integer, dimension(norbitq4, 4) :: equilist
21602114

2161-
integer, dimension(12,4) :: permutations
2115+
integer, dimension(8,4) :: permutations
21622116
integer, dimension(4) :: qplet, qplet_perm, qplet_sym
21632117
integer :: q1, q2, q3, q4, nall4, equiv, iperm, isym
21642118
integer :: iq, nsym
@@ -2171,9 +2125,8 @@ subroutine recognize_q_quadruplet(nref4, norbitq4, q_list, mappings, verbose, &
21712125

21722126
permutations = reshape([ &
21732127
1,2,3,4, 2,1,4,3, 3,4,1,2, 4,3,2,1, &
2174-
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1, &
2175-
4,2,3,1, 2,4,1,3, 3,1,4,2, 1,3,2,4], &
2176-
[12,4], order=[2,1])
2128+
1,3,2,4, 3,1,4,2, 2,4,1,3, 4,2,3,1], &
2129+
[8,4], order=[2,1])
21772130

21782131
ref4 = 0
21792132
nall4 = 0
@@ -2190,7 +2143,7 @@ subroutine recognize_q_quadruplet(nref4, norbitq4, q_list, mappings, verbose, &
21902143
if (its_cons) then
21912144
ref4 = ref4 + 1
21922145
equiv = 0
2193-
do iperm = 1, 12
2146+
do iperm = 1, 8
21942147
qplet_perm(1) = qplet(permutations(iperm,1))
21952148
qplet_perm(2) = qplet(permutations(iperm,2))
21962149
qplet_perm(3) = qplet(permutations(iperm,3))
@@ -2212,7 +2165,7 @@ subroutine recognize_q_quadruplet(nref4, norbitq4, q_list, mappings, verbose, &
22122165
norbit(ref4) = equiv
22132166
else
22142167
equiv = 0
2215-
do iperm = 1, 12
2168+
do iperm = 1, 8
22162169
qplet_perm(1) = qplet(permutations(iperm,1))
22172170
qplet_perm(2) = qplet(permutations(iperm,2))
22182171
qplet_perm(3) = qplet(permutations(iperm,3))

0 commit comments

Comments
 (0)