• R/O
  • HTTP
  • SSH
  • HTTPS

pwdft: Commit

Source code and sample for Educational-PWDFT


Commit MetaInfo

Revisionb5c65176a1eed97a33574a243b9e14da0367aa4b (tree)
Zeit2018-12-13 12:14:23
Autormitsuaki1987 <kawamitsuaki@gmai...>
Commitermitsuaki1987

Log Message

Broyden's method is unified to the description of slide
FFTW is called with the modern interface

Ändern Zusammenfassung

Diff

--- a/src/fftw_wrapper.F90
+++ b/src/fftw_wrapper.F90
@@ -1,20 +1,22 @@
11 module fftw_wrapper
22 !
3+ use, intrinsic :: iso_c_binding
4+ !
35 implicit none
46 !
5- integer(8),save :: &
7+ include 'fftw3.f03'
8+ !
9+ type(c_ptr),save :: &
610 & plan_g2r, & !< FFTW plan for F(G)e^{iGr} -> f(r) (backward)
711 & plan_r2g !< FFTW plan for f(r)e^{-iGr} -> F(G) (forward)
812 !
913 integer,allocatable,save :: &
1014 & w2r(:) !< (g_wf%npw) g_wf%npw -> g_rh%nr
1115 !
12- complex(8),allocatable,save :: &
16+ complex(c_double_complex),allocatable,save :: &
1317 & fft_in(:), & !< (g_rh%nr) FFT buffer
1418 & fft_out(:) !< (g_rh%nr) FFT buffer
1519 !
16- include "fftw3.f"
17- !
1820 private plan_g2r, plan_r2g, fft_in, fft_out
1921 !
2022 contains
@@ -24,16 +26,17 @@ contains
2426 use gvec, only : g_rh, g_wf
2527 !
2628 integer :: ipw, igv(3)
29+ integer(c_int) :: fftwdim(3)
30+ !
31+ fftwdim(1:3) = (/g_rh%nft(3), g_rh%nft(2), g_rh%nft(1)/)
2732 !
2833 allocate(fft_in(g_rh%nr), fft_out(g_rh%nr))
2934 !
30- call dfftw_plan_dft_3d(plan_G2r, g_rh%nft(1), g_rh%nft(2), g_rh%nft(3), &
31- & fft_in, fft_out, &
32- & fftw_backward, fftw_estimate)
35+ plan_G2r = fftw_plan_dft(3, fftwdim, fft_in, fft_out, &
36+ & fftw_backward, fftw_estimate)
3337 !
34- call dfftw_plan_dft_3d(plan_r2G, g_rh%nft(1), g_rh%nft(2), g_rh%nft(3), &
35- & fft_in, fft_out, &
36- & fftw_forward, fftw_estimate)
38+ plan_r2G = fftw_plan_dft(3, fftwdim, fft_in, fft_out, &
39+ & fftw_forward, fftw_estimate)
3740 !
3841 allocate(w2r(g_wf%npw))
3942 do ipw = 1, g_wf%npw
@@ -54,7 +57,7 @@ contains
5457 complex(8),intent(out) :: VlocG(g_rh%nr)
5558 !
5659 fft_in(1:g_rh%nr) = VlocR(1:g_rh%nr)
57- call dfftw_execute_dft(plan_r2g, fft_in, fft_out)
60+ call fftw_execute_dft(plan_r2g, fft_in, fft_out)
5861 VlocG(1:g_rh%nr) = fft_out(1:g_rh%nr) / dble(g_rh%nr)
5962 !
6063 end subroutine fft_r2g
@@ -69,7 +72,7 @@ contains
6972 real(8),intent(out) :: VlocR(g_rh%nr)
7073 !
7174 fft_in(1:g_rh%nr) = VlocG(1:g_rh%nr)
72- call dfftw_execute_dft(plan_g2r, fft_in, fft_out)
75+ call fftw_execute_dft(plan_g2r, fft_in, fft_out)
7376 VlocR(1:g_rh%nr) = dble(fft_out(1:g_rh%nr))
7477 !
7578 end subroutine fft_g2r
@@ -85,7 +88,7 @@ contains
8588 complex(8),intent(out) :: wfG(g_wf%npw)
8689 !
8790 fft_in(1:g_rh%nr) = wfR(1:g_rh%nr)
88- call dfftw_execute_dft(plan_r2g, fft_in, fft_out)
91+ call fftw_execute_dft(plan_r2g, fft_in, fft_out)
8992 wfG(1:g_wf%npw) = fft_out(w2r(1:g_wf%npw)) * sqrt(Vcell) / dble(g_rh%nr)
9093 !
9194 end subroutine fft_r2g_w
@@ -102,7 +105,7 @@ contains
102105 !
103106 fft_in(1:g_rh%nr) = 0.0d0
104107 fft_in(w2r(1:g_wf%npw)) = wfG(1:g_wf%npw)
105- call dfftw_execute_dft(plan_g2r, fft_in, fft_out)
108+ call fftw_execute_dft(plan_g2r, fft_in, fft_out)
106109 wfR(1:g_rh%nr) = fft_out(1:g_rh%nr) / sqrt(Vcell)
107110 !
108111 end subroutine fft_g2r_w
--- a/src/scf.F90
+++ b/src/scf.F90
@@ -17,56 +17,50 @@ contains
1717 use constant, only : htr2ev
1818 !
1919 integer :: istep, jstep
20- real(8) :: alpha, res
21- real(8),allocatable :: jacob1(:,:), jacob2(:,:), rhs0(:), dVks(:), rhs(:), drhs(:)
20+ logical :: linit
21+ real(8) :: alpha, res, &
22+ & Fvec0(g_rh%nr), dVks(g_rh%nr), Fvec(g_rh%nr), dFvec(g_rh%nr), &
23+ & jacob1(g_rh%nr,2:electron_maxstep), jacob2(g_rh%nr,2:electron_maxstep)
2224 !
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))
25- !
26- istep = 0
27- write(*,*) " Iteration ", istep
28- !
29- call kohn_sham_eq(.true., rhs)
30- res = sqrt(dot_product(rhs(1:g_rh%nr), rhs(1:g_rh%nr))) / dble(g_rh%nr)
31- write(*,*) " delta Vks [eV] : ", res * htr2eV
32- !
33- if(res < conv_thr) electron_maxstep = 0
34- !
35- dVks(1:g_rh%nr) = - mixing_beta * rhs(1:g_rh%nr)
25+ linit = .true.
3626 !
3727 do istep = 1, electron_maxstep
3828 !
3929 write(*,*) " Iteration ", istep
4030 !
41- Vks(1:g_rh%nr) = Vks(1:g_rh%nr) + dVks(1:g_rh%nr)
42- !
43- rhs0(1:g_rh%nr) = rhs(1:g_rh%nr)
44- call kohn_sham_eq(.false., rhs)
45- res = sqrt(dot_product(rhs(1:g_rh%nr), rhs(1:g_rh%nr))) / dble(g_rh%nr)
31+ call kohn_sham_eq(linit, Fvec)
32+ linit = .false.
33+ res = sqrt(dot_product(Fvec(1:g_rh%nr), Fvec(1:g_rh%nr))) / dble(g_rh%nr)
4634 write(*,*) " delta Vks [eV] : ", res * htr2eV
47- !
4835 if(res < conv_thr) exit
4936 !
50- ! Update Jacobian with drhs
51- !
52- drhs(1:g_rh%nr) = rhs(1:g_rh%nr) - rhs0(1:g_rh%nr)
53- !
54- jacob1(1:g_rh%nr,istep) = - mixing_beta * drhs(1:g_rh%nr)
55- do jstep = 1, istep - 1
56- alpha = dot_product(jacob2(1:g_rh%nr,jstep), drhs(1:g_rh%nr))
57- jacob1(1:g_rh%nr,istep) = jacob1(1:g_rh%nr,istep) - jacob1(1:g_rh%nr,jstep) * alpha
58- end do
59- jacob1(1:g_rh%nr,istep) = dVks(1:g_rh%nr) + jacob1(1:g_rh%nr,istep)
60- alpha = dot_product(drhs(1:g_rh%nr), drhs(1:g_rh%nr))
61- jacob2(1:g_rh%nr,istep) = drhs(1:g_rh%nr) / alpha
37+ if(istep > 1) then
38+ !
39+ ! Update Jacobian with dFvec
40+ !
41+ dFvec(1:g_rh%nr) = Fvec(1:g_rh%nr) - Fvec0(1:g_rh%nr)
42+ !
43+ jacob1(1:g_rh%nr,istep) = mixing_beta * dFvec(1:g_rh%nr) + dVks(1:g_rh%nr)
44+ do jstep = 2, istep - 1
45+ alpha = dot_product(jacob2(1:g_rh%nr,jstep), dFvec(1:g_rh%nr))
46+ jacob1(1:g_rh%nr,istep) = jacob1(1:g_rh%nr,istep) - jacob1(1:g_rh%nr,jstep) * alpha
47+ end do
48+ !
49+ alpha = dot_product(dFvec(1:g_rh%nr), dFvec(1:g_rh%nr))
50+ jacob2(1:g_rh%nr,istep) = dFvec(1:g_rh%nr) / alpha
51+ !
52+ end if
6253 !
63- ! Compute dVks with new Jacobian & rhs
54+ ! Compute dVks with new Jacobian & Fvec
6455 !
65- dVks(1:g_rh%nr) = - mixing_beta * rhs(1:g_rh%nr)
66- do jstep = 1, istep
67- alpha = dot_product(jacob2(1:g_rh%nr,jstep), rhs(1:g_rh%nr))
56+ dVks(1:g_rh%nr) = mixing_beta * Fvec(1:g_rh%nr)
57+ do jstep = 2, istep
58+ alpha = dot_product(jacob2(1:g_rh%nr,jstep), Fvec(1:g_rh%nr))
6859 dVks(1:g_rh%nr) = dVks(1:g_rh%nr) - jacob1(1:g_rh%nr,jstep) * alpha
6960 end do
61+ Fvec0(1:g_rh%nr) = Fvec(1:g_rh%nr)
62+ !
63+ Vks(1:g_rh%nr) = Vks(1:g_rh%nr) + dVks(1:g_rh%nr)
7064 !
7165 end do ! istep
7266 !
@@ -80,7 +74,7 @@ contains
8074 !
8175 !
8276 !
83- subroutine kohn_sham_eq(linit,rhs)
77+ subroutine kohn_sham_eq(linit,Fvec)
8478 !
8579 use gvec, only : g_wf, g_rh
8680 use kohn_sham, only : nbnd, eval, evec, calculation
@@ -90,7 +84,7 @@ contains
9084 use lobpcg, only : lobpcg_main
9185 !
9286 logical,intent(in) :: linit
93- real(8),intent(out) :: rhs(g_rh%nr)
87+ real(8),intent(out) :: Fvec(g_rh%nr)
9488 !
9589 integer :: ik, istep, avestep
9690 !
@@ -115,14 +109,14 @@ contains
115109 !
116110 call ksum_rho()
117111 !
118- rhs(1:g_rh%nr) = Vps(1:g_rh%nr)
119- call hartree_pot(rhs)
120- call xc_pot(rhs)
112+ Fvec(1:g_rh%nr) = Vps(1:g_rh%nr)
113+ call hartree_pot(Fvec)
114+ call xc_pot(Fvec)
121115 !
122- rhs(1:g_rh%nr) = Vks(1:g_rh%nr) - rhs(1:g_rh%nr)
116+ Fvec(1:g_rh%nr) = Fvec(1:g_rh%nr) - Vks(1:g_rh%nr)
123117 !
124118 else
125- rhs(1:g_rh%nr) = 0.0d0
119+ Fvec(1:g_rh%nr) = 0.0d0
126120 end if
127121 !
128122 end subroutine kohn_sham_eq
Show on old repository browser