• R/O
  • HTTP
  • SSH
  • HTTPS

pwdft: Commit

Source code and sample for Educational-PWDFT


Commit MetaInfo

Revision969c826e5bb50ae8afcf0a7ab89ce98727ca394b (tree)
Zeit2018-12-20 14:53:12
Autormitsuaki1987 <kawamitsuaki@gmai...>
Commitermitsuaki1987

Log Message

Computing total energy

Ändern Zusammenfassung

Diff

--- a/src/Makefile
+++ b/src/Makefile
@@ -4,6 +4,7 @@ OBJS=\
44 atm_spec.o \
55 constant.o \
66 diag_direct.o \
7+energy.o \
78 fftw_wrapper.o \
89 griddata.o \
910 gvec.o \
--- /dev/null
+++ b/src/energy.F90
@@ -0,0 +1,257 @@
1+module energy
2+ !
3+ implicit none
4+ !
5+contains
6+ !
7+ subroutine total_e()
8+ !
9+ use constant, only : htr2ev
10+ !
11+ real(8) :: Etot
12+ !
13+ Etot = 0.0d0
14+ !
15+ write(*,*) " Energy per u.c. [eV]"
16+ !
17+ call kinetic(Etot)
18+ call hartree(Etot)
19+ call atomic(Etot)
20+ call ewald(Etot)
21+ call xc(Etot)
22+ !
23+ write(*,*) " Total energy : ", Etot*htr2ev
24+ !
25+ end subroutine total_e
26+ !
27+ subroutine kinetic(Etot)
28+ !
29+ use constant, only : htr2ev
30+ use atm_spec, only : bvec
31+ use k_point, only : nk, kvec, kgrd
32+ use kohn_sham, only : nbnd, ef, eval, evec
33+ use gvec, only : g_wf
34+ use libtetrabz, only : libtetrabz_occ
35+ !
36+ real(8),intent(inout) :: Etot
37+ !
38+ integer :: ik, ibnd, ipw
39+ real(8) :: occ(nbnd,nk), kgv(3), gg05(g_wf%npw), Ekin
40+ !
41+ eval(1:nbnd,1:nk) = eval(1:nbnd,1:nk) - ef
42+ call libtetrabz_occ(2,bvec,nbnd,kgrd,eval,kgrd,occ)
43+ eval(1:nbnd,1:nk) = eval(1:nbnd,1:nk) + ef
44+ !
45+ Ekin = 0.0d0
46+ do ik = 1, nk
47+ !
48+ do ipw = 1, g_wf%npw
49+ kgv(1:3) = kvec(1:3,ik) + matmul(bvec(1:3,1:3), dble(g_wf%mill(1:3,g_wf%map(ipw))))
50+ gg05(ipw) = 0.5d0 * dot_product(kgv,kgv)
51+ end do
52+ !
53+ do ibnd = 1, nbnd
54+ Ekin = Ekin + occ(ibnd,ik) &
55+ & * dble(sum(conjg(evec(1:g_wf%npw,ibnd,ik)) &
56+ & * evec(1:g_wf%npw,ibnd,ik) &
57+ & * gg05(1:g_wf%npw ) ))
58+ end do
59+ end do
60+ !
61+ ! Spin
62+ !
63+ Ekin = Ekin * 2.0d0
64+ !
65+ write(*,*) " Kinetic energy : ", Ekin*htr2ev
66+ Etot = Etot + Ekin
67+ !
68+ end subroutine kinetic
69+ !
70+ subroutine hartree(Etot)
71+ !
72+ use constant, only : pi, htr2ev
73+ use gvec, only : g_rh
74+ use atm_spec, only : bvec, Vcell
75+ use fftw_wrapper, only : fft_r2g
76+ use rho_v, only : rho
77+ !
78+ real(8),intent(inout) :: Etot
79+ !
80+ integer :: ir
81+ real(8) :: g3(3), EH
82+ complex(8) :: rhog(g_rh%nr)
83+ !
84+ call fft_r2G(rho, rhoG)
85+ !
86+ ! G = 0 : Compensation to the ionic potential
87+ !
88+ EH = 0.0d0
89+ do ir = 2, g_rh%nr
90+ g3(1:3) = matmul(bvec(1:3,1:3), dble(g_rh%mill(1:3,ir)))
91+ EH = EH + 4.0d0 * pi * dble(conjg(rhog(ir))*rhog(ir)) &
92+ & / dot_product(g3(1:3),g3(1:3))
93+ end do
94+ EH = EH * 0.5d0 * Vcell
95+ !
96+ write(*,*) " Hartree energy : ", EH*htr2ev
97+ Etot = Etot + EH
98+ !
99+ end subroutine hartree
100+ !
101+ subroutine atomic(Etot)
102+ !
103+ use constant, only : htr2ev
104+ use rho_v, only : rho, Vps
105+ use gvec, only : g_rh
106+ use atm_spec, only : Vcell
107+ !
108+ real(8),intent(inout) :: Etot
109+ !
110+ real(8) :: Eps
111+ !
112+ Eps = dot_product(rho(1:g_rh%nr), Vps(1:g_rh%nr)) * Vcell / dble(g_rh%nr)
113+ !
114+ write(*,*) " rho*V : ", Eps*htr2ev
115+ Etot = Etot + Eps
116+ !
117+ end subroutine atomic
118+ !
119+ subroutine ewald(Etot)
120+ !
121+ use constant, only : htr2ev, pi
122+ use atm_spec, only : avec, bvec, Vcell, atm, spec, nat, nelec
123+ !
124+ real(8),intent(inout) :: Etot
125+ !
126+ integer :: nmax3(3), iat, jat, i1, i2, i3, ii
127+ real(8) :: Eew, eta, cutoff, dtau(3), ZZ, gv(3), g2, phase, norm, rv(3)
128+ !
129+ Eew = 0.0d0
130+ !
131+ eta = 0.0d0
132+ do ii = 1, 3
133+ norm = sqrt(dot_product(bvec(1:3,ii), bvec(1:3,ii)))
134+ eta = eta + norm
135+ end do
136+ eta = eta / 3.0d0
137+ write(*,*) " Ewald eta [Bohr^-1] : ", eta
138+ !
139+ ! Reciprocal space
140+ !
141+ cutoff = 10.0d0 * eta * 2.0d0
142+ do ii = 1, 3
143+ norm = sqrt(dot_product(avec(1:3,ii), avec(1:3,ii)))
144+ nmax3(ii) = ceiling(cutoff*norm/(2.0d0*pi))
145+ end do
146+ !
147+ write(*,*) " Ewald grid (G) : ", nmax3(1:3)
148+ !
149+ do iat = 1, nat
150+ do jat = 1, nat
151+ !
152+ dtau(1:3) = atm(iat)%pos(1:3) - atm(jat)%pos(1:3)
153+ ZZ = spec(atm(iat)%ityp)%Zion * spec(atm(jat)%ityp)%Zion
154+ !
155+ do i3 = -nmax3(3), nmax3(3)
156+ do i2 = -nmax3(2), nmax3(2)
157+ do i1 = -nmax3(1), nmax3(1)
158+ !
159+ if(all((/i1,i2,i3/)==0)) cycle
160+ !
161+ gv(1:3) = matmul(bvec(1:3,1:3), dble((/i1,i2,i3/)))
162+ g2 = dot_product(gv,gv)
163+ phase = 2.0d0*pi*dot_product(dtau(1:3), dble((/i1,i2,i3/)))
164+ !
165+ Eew = Eew &
166+ & + 0.5d0 * ZZ * cos(phase) * 4.0d0 * pi * exp(-g2/(4.0d0*eta**2)) / (Vcell*g2)
167+ !
168+ end do
169+ end do
170+ end do
171+ end do
172+ end do
173+ !
174+ Eew = Eew - 0.5d0 * pi * nelec**2 / (Vcell*eta**2)
175+ !
176+ ! Real space
177+ !
178+ cutoff = 10.0d0 / eta
179+ do ii = 1, 3
180+ norm = sqrt(dot_product(bvec(1:3,ii), bvec(1:3,ii)))
181+ nmax3(ii) = ceiling(cutoff*norm/(2.0d0*pi))
182+ end do
183+ !
184+ write(*,*) " Ewald grid (R) : ", nmax3(1:3)
185+ !
186+ do iat = 1, nat
187+ do jat = 1, nat
188+ !
189+ dtau(1:3) = atm(iat)%pos(1:3) - atm(jat)%pos(1:3)
190+ ZZ = spec(atm(iat)%ityp)%Zion * spec(atm(jat)%ityp)%Zion
191+ !
192+ do i3 = -nmax3(3), nmax3(3)
193+ do i2 = -nmax3(2), nmax3(2)
194+ do i1 = -nmax3(1), nmax3(1)
195+ !
196+ if(all((/i1,i2,i3/)==0) .and. iat == jat) cycle
197+ !
198+ rv(1:3) = matmul(avec(1:3,1:3), dtau(1:3) + dble((/i1,i2,i3/)))
199+ norm = sqrt(dot_product(rv(1:3), rv(1:3)))
200+ !
201+ Eew = Eew &
202+ & + 0.5d0 * ZZ * erfc(norm*eta) / norm
203+ !
204+ end do
205+ end do
206+ end do
207+ end do
208+ !
209+ Eew = Eew - eta * spec(atm(iat)%ityp)%Zion**2 / sqrt(pi)
210+ !
211+ end do
212+ !
213+ write(*,*) " Eward energy : ", Eew*htr2ev
214+ Etot = Etot + Eew
215+ !
216+ end subroutine ewald
217+ !
218+ subroutine xc(Etot)
219+ !
220+ use constant, only : htr2ev, pi
221+ use gvec, only : g_rh
222+ use rho_v, only : rho
223+ use atm_spec, only : Vcell
224+ !
225+ real(8),intent(inout) :: Etot
226+ !
227+ integer :: ir
228+ real(8) :: rs, exc0, Exc
229+ !
230+ Exc = 0.0d0
231+ do ir = 1, g_rh%nr
232+ !
233+ if(rho(ir) > 1.0d-10) then
234+ !
235+ rs = (3.0d0 / (4.0d0 * pi * rho(ir)))**(1.0d0/3.0d0)
236+ !
237+ if(rs < 1.0d0) then
238+ exc0 = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
239+ & - 0.0480d0 + 0.031d0*log(rs) - 0.0116d0*rs + 0.0020d0*rs*log(rs)
240+ else
241+ exc0 = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
242+ & - 0.1423d0 / (1.0d0 + 1.0529d0*sqrt(rs) + 0.3334d0*rs)
243+ end if
244+ !
245+ Exc = Exc + exc0 * rho(ir)
246+ !
247+ end if
248+ !
249+ end do
250+ Exc = Exc * Vcell / dble(g_rh%nr)
251+ !
252+ write(*,*) " rho*V : ", Exc*htr2ev
253+ Etot = Etot + Exc
254+ !
255+ end subroutine xc
256+ !
257+end module energy
--- a/src/k_point.F90
+++ b/src/k_point.F90
@@ -22,7 +22,7 @@ contains
2222 real(8) :: occ(nbnd,nk)
2323 complex(8) :: psir(g_rh%nr)
2424 !
25- CALL libtetrabz_fermieng(2,bvec,nbnd,kgrd,eval,kgrd,occ,ef,nelec*0.5d0)
25+ call libtetrabz_fermieng(2,bvec,nbnd,kgrd,eval,kgrd,occ,ef,nelec*0.5d0)
2626 !
2727 rho(1:g_rh%nr) = 0.0d0
2828 do ik = 1, nk
--- a/src/make.depend
+++ b/src/make.depend
@@ -1,6 +1,7 @@
11 atm_spec.o : atm_spec.F90
22 constant.o : constant.F90
33 diag_direct.o : diag_direct.F90 fftw_wrapper.o rho_v.o gvec.o kohn_sham.o atm_spec.o
4+energy.o : energy.F90 rho_v.o fftw_wrapper.o gvec.o kohn_sham.o k_point.o atm_spec.o constant.o
45 fftw_wrapper.o : fftw_wrapper.F90 atm_spec.o gvec.o
56 griddata.o : griddata.F90 gvec.o atm_spec.o constant.o
67 gvec.o : gvec.F90 atm_spec.o constant.o
@@ -10,7 +11,7 @@ kohn_sham.o : kohn_sham.F90
1011 lobpcg.o : lobpcg.F90 gvec.o atm_spec.o hamiltonian.o kohn_sham.o
1112 plot.o : plot.F90 gvec.o atm_spec.o kohn_sham.o k_point.o constant.o
1213 pp.o : pp.F90 kohn_sham.o atm_spec.o
13-pwdft.o : pwdft.F90 plot.o griddata.o fftw_wrapper.o scf.o pp.o k_point.o rho_v.o gvec.o constant.o stdin.o kohn_sham.o
14+pwdft.o : pwdft.F90 energy.o plot.o griddata.o fftw_wrapper.o scf.o pp.o k_point.o rho_v.o gvec.o constant.o stdin.o kohn_sham.o
1415 rho_v.o : rho_v.F90 fftw_wrapper.o constant.o atm_spec.o gvec.o
1516 scf.o : scf.F90 lobpcg.o diag_direct.o k_point.o kohn_sham.o constant.o rho_v.o gvec.o
1617 stdin.o : stdin.F90 gvec.o atm_spec.o k_point.o scf.o kohn_sham.o constant.o
--- a/src/pwdft.F90
+++ b/src/pwdft.F90
@@ -18,6 +18,7 @@ program pwdft
1818 use fftw_wrapper, only : init_fft
1919 use griddata, only : read_griddata, write_griddata
2020 use plot, only : fermi_plot, band_plot
21+ use energy, only : total_e
2122 !
2223 implicit none
2324 !
@@ -42,6 +43,8 @@ program pwdft
4243 call cpu_time(t4)
4344 write(*,*) " SCF time : ", t4 - t3, " sec"
4445 !
46+ call total_e()
47+ !
4548 Vks(1:g_rh%nr) = (Vks(1:g_rh%nr) - ef)*htr2ev
4649 call write_griddata("vks.xsf", Vks)
4750 Vks(1:g_rh%nr) = Vks(1:g_rh%nr)/htr2ev + ef
Show on old repository browser