• R/O
  • HTTP
  • SSH
  • HTTPS

pwdft: Commit

Source code and sample for Educational-PWDFT


Commit MetaInfo

Revisionbde80125bdf9c5d66be234a908dc7dfb6c983cb0 (tree)
Zeit2018-12-07 18:43:57
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Bag fix: LOBPCG eval exceeded the array size in diag_ovlp
In the 2nd SCF, the previous WF was not used.

Ändern Zusammenfassung

Diff

--- a/sample/Al/bands.in
+++ b/sample/Al/bands.in
@@ -19,7 +19,6 @@ ATOMIC_SPECIES
1919 ATOMIC_POSITIONS
2020 Al 0.000000 0.000000 0.000000
2121 K_POINTS
22-3
22+2
2323 0.0 0.0 0.0
2424 0.5 0.0 0.0
25-0.0 0.0 0.0
--- a/sample/Al/scf.in
+++ b/sample/Al/scf.in
@@ -2,14 +2,16 @@
22 calculation = 'scf'
33 /
44 &SYSTEM
5+nbnd = 5
56 nat = 1
67 ntyp = 1
7- ecutwfc = 30.000000
8+ ecutwfc =30.000000
89 ecutrho = 120.000000
910 /
1011 &ELECTRONS
1112 mixing_beta = 0.3
12- conv_thr = 1.000000e-5
13+conv_thr = 1.000000e-5
14+electron_maxstep = 100
1315 /
1416 CELL_PARAMETERS
1517 0.000000 2.024700 2.024700
--- a/src/lobpcg.F90
+++ b/src/lobpcg.F90
@@ -31,11 +31,11 @@ contains
3131 !
3232 end subroutine initialize
3333 !
34- subroutine diag_ovrp(nsub,hsub,ovlp,eval)
34+ subroutine diag_ovrp(nsub,hsub,ovlp,esub)
3535 !
3636 integer,intent(in) :: nsub
3737 complex(8),intent(inout) :: hsub(nsub,nsub), ovlp(nsub,nsub)
38- real(8),intent(out) :: eval(nsub)
38+ real(8),intent(out) :: esub(nsub)
3939 !
4040 integer :: lwork, isub, nsub2, info
4141 real(8) :: rwork(3*nsub-2)
@@ -43,18 +43,18 @@ contains
4343 !
4444 lwork = -1
4545 allocate(work(1))
46- call zheev('V', 'U', nsub, ovlp, nsub, eval, work, lwork, rwork, info)
46+ call zheev('V', 'U', nsub, ovlp, nsub, esub, work, lwork, rwork, info)
4747 lwork = nint(dble(work(1)))
4848 deallocate(work)
4949 allocate(work(lwork))
50- call zheev('V', 'U', nsub, ovlp, nsub, eval, work, lwork, rwork, info)
50+ call zheev('V', 'U', nsub, ovlp, nsub, esub, work, lwork, rwork, info)
5151 deallocate(work)
5252 !
5353 nsub2 = 0
5454 do isub = 1, nsub
55- if(eval(isub) > 1.0d-14) then
55+ if(esub(isub) > 1.0d-14) then
5656 nsub2 = nsub2 + 1
57- ovlp(1:nsub,nsub2) = ovlp(1:nsub,isub) / sqrt(eval(isub))
57+ ovlp(1:nsub,nsub2) = ovlp(1:nsub,isub) / sqrt(esub(isub))
5858 end if
5959 end do
6060 ovlp(1:nsub,nsub2+1:nsub) = 0.0d0
@@ -64,11 +64,11 @@ contains
6464 !
6565 lwork = -1
6666 allocate(work(1))
67- call zheev('V', 'U', nsub2, hsub, nsub, eval, work, lwork, rwork, info)
67+ call zheev('V', 'U', nsub2, hsub, nsub, esub, work, lwork, rwork, info)
6868 lwork = nint(dble(work(1)))
6969 deallocate(work)
7070 allocate(work(lwork))
71- call zheev('V', 'U', nsub2, hsub, nsub, eval, work, lwork, rwork, info)
71+ call zheev('V', 'U', nsub2, hsub, nsub, esub, work, lwork, rwork, info)
7272 deallocate(work)
7373 !
7474 hsub(1:nsub, 1:nsub2) = matmul(ovlp(1:nsub,1:nsub2), hsub(1:nsub2,1:nsub2))
@@ -76,7 +76,7 @@ contains
7676 !
7777 end subroutine diag_ovrp
7878 !
79- subroutine lobpcg_main(linit,npw,kvec,evec,eval)
79+ subroutine lobpcg_main(linit,npw,kvec,evec,eval,istep)
8080 !
8181 use kohn_sham, only : nbnd, calculation
8282 use hamiltonian, only : h_psi
@@ -88,9 +88,11 @@ contains
8888 real(8),intent(in) :: kvec(3)
8989 real(8),intent(out) :: eval(nbnd)
9090 complex(8),intent(out) :: evec(npw,nbnd)
91+ integer,intent(out) :: istep
9192 !
92- integer :: ii, ibnd, istep, nsub, ipw, cg_maxstep = 100
93- real(8) :: norm, maxnorm, ekin(npw), pre(npw), gv(3), ekin0, cg_thr = 1.0d-4
93+ integer :: ii, ibnd, nsub, ipw, cg_maxstep = 100
94+ real(8) :: norm, maxnorm, ekin(npw), pre(npw), gv(3), ekin0, &
95+ & cg_thr = 1.0d-4, esub(3*nbnd)
9496 complex(8) :: wxp(npw,nbnd,3), hwxp(npw,nbnd,3), xp(npw,nbnd,2:3), &
9597 & hsub(nbnd,3,3*nbnd), ovlp(3*nbnd,3*nbnd), rotmat(nbnd,3,nbnd,2:3)
9698 !
@@ -104,12 +106,13 @@ contains
104106 wxp( 1:npw,1:nbnd,1:3) = 0.0d0
105107 hwxp(1:npw,1:nbnd,1:3) = 0.0d0
106108 !
107- if(linit) call initialize(wxp(1:npw,1:nbnd,2))
109+ if(linit) call initialize(evec(1:npw,1:nbnd))
110+ wxp(1:npw,1:nbnd,2) = evec(1:npw,1:nbnd)
108111 !
109112 call h_psi(kvec,wxp(1:npw,1:nbnd,2), hwxp(1:npw,1:nbnd,2))
110113 !
111114 do ibnd = 1, nbnd
112- eval(ibnd) = dble(dot_product(wxp(1:npw,ibnd,2), hwxp(1:npw,ibnd,2)))
115+ esub(ibnd) = dble(dot_product(wxp(1:npw,ibnd,2), hwxp(1:npw,ibnd,2)))
113116 end do
114117 !
115118 if(calculation=="iterative") write(*,*) "Step Residual"
@@ -119,7 +122,7 @@ contains
119122 maxnorm = 0.0d0
120123 do ibnd = 1, nbnd
121124 !
122- wxp(1:npw,ibnd,1) = hwxp(1:npw,ibnd,2) - eval(ibnd) * wxp(1:npw,ibnd,2)
125+ wxp(1:npw,ibnd,1) = hwxp(1:npw,ibnd,2) - esub(ibnd) * wxp(1:npw,ibnd,2)
123126 norm = sqrt(dble(dot_product(wxp(1:npw,ibnd,1), wxp(1:npw,ibnd,1))))
124127 maxnorm = max(norm, maxnorm)
125128 !
@@ -148,7 +151,7 @@ contains
148151 call zgemm("C", "N", 3*nbnd, 3*nbnd, npw, &
149152 & (1.0d0,0.0d0), wxp, npw, wxp, npw, (0.0d0,0.0d0), ovlp, 3*nbnd)
150153 !
151- call diag_ovrp(nsub,hsub,ovlp,eval)
154+ call diag_ovrp(nsub,hsub,ovlp,esub)
152155 !
153156 rotmat(1:nbnd,1:3,1:nbnd,2) = hsub(1:nbnd,1:3,1:nbnd)
154157 rotmat(1:nbnd, 1,1:nbnd,3) = hsub(1:nbnd, 1,1:nbnd)
@@ -173,10 +176,11 @@ contains
173176 end do
174177 !
175178 if(istep >= cg_maxstep) then
176- write(*,*) "Not converged at kvec : ", kvec(1:3), maxnorm
179+ write(*,*) " Not converged at kvec : ", kvec(1:3), ", norm : ", maxnorm
177180 end if
178181 !
179182 evec(1:npw,1:nbnd) = wxp(1:npw,1:nbnd,2)
183+ eval(1:nbnd) = esub(1:nbnd)
180184 !
181185 end subroutine lobpcg_main
182186 !
--- a/src/scf.F90
+++ b/src/scf.F90
@@ -17,41 +17,35 @@ contains
1717 use constant, only : htr2ev
1818 !
1919 integer :: istep, jstep
20- real(8) :: jacob1(g_rh%nr,electron_maxstep), &
21- & jacob2(g_rh%nr,electron_maxstep), &
22- & alpha
23- real(8) :: rhs0(g_rh%nr)
24- real(8) :: dVks(g_rh%nr)
25- real(8) :: res, rhs(g_rh%nr)
26- real(8) :: drhs(g_rh%nr)
20+ real(8) :: alpha, res
21+ real(8),allocatable :: jacob1(:,:), jacob2(:,:), rhs0(:), dVks(:), rhs(:), drhs(:)
22+ !
23+ allocate(rhs0(g_rh%nr),dVks(g_rh%nr),rhs(g_rh%nr), drhs(g_rh%nr), &
24+ & jacob1(g_rh%nr,electron_maxstep), jacob2(g_rh%nr,electron_maxstep))
2725 !
2826 istep = 0
29- write(*,*) "Iteration ", istep
27+ write(*,*) " Iteration ", istep
3028 !
3129 call kohn_sham_eq(.true., rhs)
3230 res = sqrt(dot_product(rhs(1:g_rh%nr), rhs(1:g_rh%nr))) / dble(g_rh%nr)
33- write(*,*) " delta Vks [eV] : ", res * htr2eV
31+ write(*,*) " delta Vks [eV] : ", res * htr2eV
3432 !
35- if(res < conv_thr) GOTO 5
33+ if(res < conv_thr) electron_maxstep = 0
3634 !
3735 dVks(1:g_rh%nr) = - mixing_beta * rhs(1:g_rh%nr)
3836 !
3937 do istep = 1, electron_maxstep
4038 !
41- write(*,*) "Iteration ", istep
39+ write(*,*) " Iteration ", istep
4240 !
4341 Vks(1:g_rh%nr) = Vks(1:g_rh%nr) + dVks(1:g_rh%nr)
4442 !
4543 rhs0(1:g_rh%nr) = rhs(1:g_rh%nr)
4644 call kohn_sham_eq(.false., rhs)
4745 res = sqrt(dot_product(rhs(1:g_rh%nr), rhs(1:g_rh%nr))) / dble(g_rh%nr)
48- write(*,*) " delta Vks [eV] : ", res * htr2eV
46+ write(*,*) " delta Vks [eV] : ", res * htr2eV
4947 !
50- if(res < conv_thr) then
51- !
52- GOTO 5
53- !
54- end if
48+ if(res < conv_thr) exit
5549 !
5650 ! Update Jacobian with drhs
5751 !
@@ -76,12 +70,11 @@ contains
7670 !
7771 end do ! istep
7872 !
79- write(*,'(/,7x,"Not converged ! res = ",e12.5,/)') res
80- RETURN
81- !
82-5 CONTINUE
83- !
84- write(*,'(/,7x,"Converged ! iter = ",i0,/)') istep
73+ if(istep >= electron_maxstep) then
74+ write(*,'(/,7x,"Not converged ! res = ",e12.5,/)') res
75+ else
76+ write(*,'(/,7x,"Converged ! iter = ",i0,/)') istep
77+ end if
8578 !
8679 end subroutine scf_loop
8780 !
@@ -99,17 +92,20 @@ contains
9992 logical,intent(in) :: linit
10093 real(8),intent(out) :: rhs(g_rh%nr)
10194 !
102- integer :: ik
95+ integer :: ik, istep, avestep
10396 !
10497 if(calculation == "direct") then
10598 do ik = 1, nk
10699 call direct(g_wf%npw,kvec(1:3,ik),evec(1:g_wf%npw,1:nbnd,ik),eval(1:nbnd,ik))
107100 end do
108101 else
102+ avestep = 0
109103 do ik = 1, nk
110104 call lobpcg_main(linit, g_wf%npw, kvec(1:3,ik),&
111- & evec(1:g_wf%npw,1:nbnd,ik),eval(1:nbnd,ik))
105+ & evec(1:g_wf%npw,1:nbnd,ik),eval(1:nbnd,ik), istep)
106+ avestep = avestep + istep
112107 end do
108+ write(*,*) " Average LOBPCG steps : ", avestep / nk
113109 end if
114110 !
115111 if(calculation == "scf") then
Show on old repository browser