• R/O
  • HTTP
  • SSH
  • HTTPS

pwdft: Commit

Source code and sample for Educational-PWDFT


Commit MetaInfo

Revision279cbef29e152b4a0e9108e0039e1e5143ed4568 (tree)
Zeit2018-11-30 00:10:09
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Backup

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/sample/Al/band.in
@@ -0,0 +1,22 @@
1+&CONTROL
2+ calculation = 'scf'
3+/
4+&SYSTEM
5+nbnd = 5
6+ nat = 1
7+ ntyp = 1
8+ ecutwfc = 30.000000
9+ ecutrho = 120.000000
10+/
11+&ELECTRONS
12+/
13+CELL_PARAMETERS
14+ 0.000000 2.024700 2.024700
15+ 2.024700 0.000000 2.024700
16+ 2.024700 2.024700 0.000000
17+ATOMIC_SPECIES
18+ Al Al.pz-n-nc.UPF
19+ATOMIC_POSITIONS
20+ Al 0.000000 0.000000 0.000000
21+K_POINTS
22+ 8 8 8
--- /dev/null
+++ b/sample/Al/scf.in
@@ -0,0 +1,26 @@
1+&CONTROL
2+ calculation = 'scf'
3+ pseudo_dir = '/work/i0012/i001200/pseudo/'
4+/
5+&SYSTEM
6+ ibrav = 0
7+ nat = 1
8+ ntyp = 1
9+ ecutwfc = 50.000000
10+ ecutrho = 200.000000
11+ occupations = 'tetrahedra_opt'
12+/
13+&ELECTRONS
14+ mixing_beta = 0.3
15+ conv_thr = 1.000000e-10
16+/
17+CELL_PARAMETERS angstrom
18+ 0.000000 2.024700 2.024700
19+ 2.024700 0.000000 2.024700
20+ 2.024700 2.024700 0.000000
21+ATOMIC_SPECIES
22+ Al 26.981539 Al.pz-n-nc.UPF
23+ATOMIC_POSITIONS crystal
24+ Al 0.000000 0.000000 0.000000
25+K_POINTS automatic
26+ 8 8 8 0 0 0
--- a/sample/MgB2/scf.in
+++ b/sample/MgB2/scf.in
@@ -4,7 +4,8 @@
44 &SYSTEM
55 nat = 3
66 ntyp = 2
7-ecutwfc = 48.000000
7+ecutwfc = 25.000000
8+ecutrho = 100.0
89 nbnd = 10
910 /
1011 CELL_PARAMETERS
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,10 +1,13 @@
11 include ../make.inc
22
3-include make.depend
4-
53 OBJS=\
64 atm_spec.o \
5+diag_direct.o \
6+gvec.o \
7+hamiltonian.o \
8+io_vloc.o \
79 k_point.o \
10+lobpcg.o \
811 pwdft.o \
912 solver.o \
1013 stdin.o
@@ -22,3 +25,5 @@ pwdft.x:$(OBJS)
2225
2326 clean:
2427 rm -rf *.o *.mod *.x
28+
29+include make.depend
--- a/src/atm_spec.F90
+++ b/src/atm_spec.F90
@@ -13,9 +13,9 @@ module atm_spec
1313 character(256) :: ps_file
1414 end type spec_t
1515 !
16- integer :: nat, ntyp
17- real(8) :: avec(3,3)
18- type(atm_t),allocatable :: atm(:)
19- type(spec_t),allocatable :: spec(:)
16+ integer,save :: nat, ntyp
17+ real(8),save :: avec(3,3), bvec(3,3)
18+ type(atm_t),allocatable,save :: atm(:)
19+ type(spec_t),allocatable,save :: spec(:)
2020 !
2121 end module atm_spec
--- /dev/null
+++ b/src/diag_direct.F90
@@ -0,0 +1,62 @@
1+module diag_direct
2+ !
3+ implicit none
4+ !
5+contains
6+ !
7+ subroutine direct(kvec,evec,eval)
8+ !
9+ use solver, only : nbnd
10+ use gvec, only : rfft, wfft
11+ use io_vloc, only : Vloc
12+ !
13+ real(8),intent(in) :: kvec(3)
14+ complex(8),intent(out) :: evec(wfft%npw,nbnd)
15+ complex(8),intent(out) :: eval(nbnd)
16+ !
17+ integer :: info, ipw, jpw, lwork, dmill(3)
18+ real(8) :: rwork(3*wfft%npw-2), eval_full(wfft%npw), kgv(3)
19+ complex(8) :: ham(wfft%npw,wfft%npw), &
20+ & VlocG(rfft%npw3(1), rfft%npw3(2), rfft%npw3(3))
21+ complex(8),allocatable :: work(:)
22+ !
23+ include 'fftw3.f'
24+ integer(8) :: plan
25+ !
26+ ! Local potential term
27+ !
28+ call dfftw_plan_dft_3d(plan, rfft%npw3(1), rfft%npw3(2), rfft%npw3(3), Vloc, VlocG, &
29+ & fftw_forward, fftw_estimate)
30+ call dfftw_execute_dft(plan, Vloc, VlocG)
31+ call dfftw_destroy_plan(plan)
32+ !
33+ do ipw = 1, wfft%npw
34+ do jpw = 1, wfft%npw
35+ dmill(1:3) = wfft%mill(1:3,jpw) - wfft%mill(1:3,ipw)
36+ dmill(1:3) = modulo(dmill(1:3), rfft%npw3(1:3)) + 1
37+ ham(jpw,ipw) = VlocG(dmill(1), dmill(2), dmill(3))
38+ end do
39+ end do
40+ !
41+ ! Kinetic energy term
42+ !
43+ do ipw = 1, wfft%npw
44+ kgv(1:3) = kvec(1:3) + wfft%gv(1:3,ipw)
45+ ham(ipw,ipw) = ham(ipw,ipw) + 0.5d0 * dot_product(kgv,kgv)
46+ end do
47+ !
48+ lwork = -1
49+ allocate(work(1))
50+ call zheev('V', 'U', wfft%npw, ham, wfft%npw, eval_full, work, lwork, rwork, info)
51+ lwork = nint(dble(work(1)))
52+ deallocate(work)
53+ allocate(work(lwork))
54+ call zheev('V', 'U', wfft%npw, ham, wfft%npw, eval_full, work, lwork, rwork, info)
55+ deallocate(work)
56+ !
57+ eval( 1:nbnd) = eval_full( 1:nbnd)
58+ evec(1:wfft%npw,1:nbnd) = ham(1:wfft%npw,1:nbnd)
59+ !
60+ end subroutine direct
61+ !
62+end module diag_direct
--- /dev/null
+++ b/src/gvec.F90
@@ -0,0 +1,100 @@
1+module gvec
2+ !
3+ implicit none
4+ !
5+ type fft
6+ integer :: npw
7+ integer :: npw3(3)
8+ integer,allocatable :: mill(:,:)
9+ real(8),allocatable :: gv(:,:)
10+ end type fft
11+ !
12+ type(fft),save :: rfft
13+ type(fft),save :: wfft
14+ !
15+contains
16+ !
17+ subroutine setup_gvec(xfft,g2cut)
18+ !
19+ use atm_spec, only : avec, bvec
20+ !
21+ real(8),intent(in) :: g2cut
22+ type(fft),intent(out) :: xfft
23+ !
24+ integer :: ii, i1, i2, i3, nmax3(3), nmax
25+ real(8) :: pi = acos(-1.0d0), norm, gv0(3), g2
26+ !
27+ do ii = 1, 3
28+ norm = sqrt(dot_product(avec(1:3,ii), avec(1:3,ii)))
29+ nmax3(ii) = ceiling(sqrt(g2cut)*norm/(2.0d0*pi))
30+ end do
31+ !
32+ nmax = product(nmax3(1:3)*2+1)
33+ allocate(xfft%gv(3,nmax), xfft%mill(3,nmax))
34+ !
35+ xfft%npw = 0
36+ do i1 = -nmax3(1), nmax3(1)
37+ do i2 = -nmax3(2), nmax3(2)
38+ do i3 = -nmax3(3), nmax3(3)
39+ gv0(1:3) = matmul(bvec(1:3,1:3), dble((/i1,i2,i3/)))
40+ g2 = dot_product(gv0,gv0)
41+ if(g2 < g2cut) then
42+ xfft%npw = xfft%npw + 1
43+ xfft%mill(1:3,xfft%npw) = (/i1,i2,i3/)
44+ xfft%gv(1:3,xfft%npw) = gv0(1:3)
45+ end if
46+ end do
47+ end do
48+ end do
49+ !
50+ write(*,*) " Numver of PW : ", xfft%npw
51+ !
52+ ! Find FFT grid
53+ !
54+ do ii = 1, 3
55+ xfft%npw3(ii) = maxval(xfft%mill(ii,1:xfft%npw))*2 + 1
56+ call base2357(xfft%npw3(ii))
57+ end do
58+ !
59+ write(*,*) " FFT grid : ", xfft%npw3(1:3)
60+ !
61+ end subroutine setup_gvec
62+ !
63+ subroutine base2357(nfft)
64+ !
65+ integer,intent(inout) :: nfft
66+ !
67+ integer :: base(4) = (/2,3,5,7/), ibase, iexp, &
68+ & start, last, nexp, nfft1, nfft2
69+ !
70+ start = nfft
71+ last = 2 * nfft
72+ !
73+ do nfft1 = start, last
74+ !
75+ nfft2 = nfft1
76+ !
77+ do ibase = 1, 4
78+ !
79+ nexp = ceiling(log(dble(nfft2))/log(dble(base(ibase))))
80+ !
81+ do iexp = 1, nexp
82+ if(mod(nfft2, base(ibase)) == 0) then
83+ nfft2 = nfft2 / base(ibase)
84+ else
85+ exit
86+ end if
87+ end do
88+ !
89+ if(nfft2 == 1) then
90+ nfft = nfft1
91+ return
92+ end if
93+ !
94+ end do
95+ !
96+ end do
97+ !
98+ end subroutine base2357
99+ !
100+end module gvec
--- /dev/null
+++ b/src/hamiltonian.F90
@@ -0,0 +1,30 @@
1+module hamiltonian
2+ !
3+ implicit none
4+ !
5+contains
6+ !
7+ subroutine h_psi(kvec,psi, hpsi)
8+ !
9+ use gvec, only : wfft
10+ use solver, only : nbnd
11+ !
12+ real(8),intent(in) :: kvec(3)
13+ complex(8),intent(in) :: psi(wfft%npw,nbnd)
14+ complex(8),intent(out) :: hpsi(wfft%npw,nbnd)
15+ !
16+ integer :: ipw
17+ real(8) :: kgv(3)
18+ !
19+ hpsi(1:wfft%npw,1:nbnd) = 0.0d0
20+ !
21+ ! Kinetic energy term
22+ !
23+ do ipw = 1, wfft%npw
24+ kgv(1:3) = kvec(1:3) + wfft%gv(1:3,ipw)
25+ hpsi(ipw,1:nbnd) = 0.5d0 * dot_product(kgv,kgv) * psi(ipw,1:nbnd)
26+ end do
27+ !
28+ end subroutine h_psi
29+ !
30+end module hamiltonian
--- /dev/null
+++ b/src/io_vloc.F90
@@ -0,0 +1,64 @@
1+module io_vloc
2+ !
3+ implicit none
4+ !
5+ complex(8),allocatable :: Vloc(:,:,:)
6+ !
7+contains
8+ !
9+ subroutine read_vloc()
10+ !
11+ use atm_spec, only : nat, atm, bvec, avec
12+ use gvec, only : rfft
13+ !
14+ integer :: itemp(3), fi = 10, iat
15+ character(256) :: ctemp
16+ real(8) :: avec0(3,3), pi = acos(-1.0d0), &
17+ & Vloc0(1:rfft%npw3(1)+1,1:rfft%npw3(2)+1,1:rfft%npw3(3)+1)
18+ !
19+ open(fi, file = "vloc.xsf")
20+ !
21+ read(fi,*) ctemp
22+ read(fi,*) ctemp
23+ read(fi,*) avec0(1:3,1:3)
24+ avec0(1:3,1:3) = avec0(1:3,1:3) / 0.529177249d0
25+ if(any(abs(avec0(1:3,1:3)-avec(1:3,1:3)) > 1.0d-3)) then
26+ write(*,*) "Error in read_vloc"
27+ write(*,*) "Direct lattice vector is different."
28+ stop 'error in read_vloc'
29+ end if
30+ read(fi,*) ctemp
31+ read(fi,*) itemp(1:2)
32+ if(nat /= itemp(1)) then
33+ write(*,*) "Error in read_vloc"
34+ write(*,*) "Number of atoms is different."
35+ stop 'error in read_vloc'
36+ end if
37+ do iat = 1, nat
38+ read(fi,*) ctemp, avec0(1:3,1)
39+ avec0(1:3,1) = avec0(1:3,1) / 0.529177249d0
40+ avec0(1:3,1) = matmul(avec0(1:3,1), bvec(1:3,1:3)) / (2.0d0*pi)
41+ if(any(abs(avec0(1:3,1)-atm(iat)%pos(1:3)) > 1.0d-3)) then
42+ write(*,*) "Error in read_vloc"
43+ write(*,*) "Position of atom ", iat, " is different."
44+ stop 'error in read_vloc'
45+ end if
46+ end do
47+ read(fi,*) ctemp
48+ read(fi,*) ctemp
49+ read(fi,*) ctemp
50+ read(fi,*) itemp(1:3)
51+ if(any(itemp(1:3) /= rfft%npw3(1:3)+1)) then
52+ write(*,*) "Error in read_vloc"
53+ write(*,*) "FFT grid is different."
54+ stop 'error in read_vloc'
55+ end if
56+ read(fi,*) avec0(1:3,1)
57+ read(fi,*) avec0(1:3,1:3)
58+ read(fi,*) Vloc0(1:itemp(1),1:itemp(2),1:itemp(3))
59+ Vloc( 1:rfft%npw3(1),1:rfft%npw3(2),1:rfft%npw3(3)) &
60+ & = Vloc0(1:rfft%npw3(1),1:rfft%npw3(2),1:rfft%npw3(3))
61+ !
62+ end subroutine read_vloc
63+ !
64+end module io_vloc
--- a/src/lobpcg.F90
+++ b/src/lobpcg.F90
@@ -4,17 +4,41 @@ module lobpcg
44 !
55 contains
66 !
7- subroutine diag_ovrp()
7+ subroutine initialize(psi)
8+ !
9+ use solver, only : nbnd
10+ use gvec, only : wfft
11+ !
12+ complex(8),intent(out) :: psi(wfft%npw,nbnd)
13+ !
14+ real(8) :: rpsi(wfft%npw,nbnd), ipsi(wfft%npw,nbnd)
15+ !
16+ call random_number(rpsi)
17+ call random_number(ipsi)
18+ !
19+ psi(1:wfft%npw,1:nbnd) = cmplx(rpsi(1:wfft%npw,1:nbnd), &
20+ & ipsi(1:wfft%npw,1:nbnd), 8)
21+ !
22+ end subroutine initialize
23+ !
24+ subroutine diag_ovrp(nsub,hsub,ovlp,eval)
825 !
926 integer,intent(in) :: nsub
1027 complex(8),intent(inout) :: hsub(nsub,nsub), ovlp(nsub,nsub)
11- complex(8),intent(out) :: eval(nsub)
28+ real(8),intent(out) :: eval(nsub)
1229 !
13- liwork = 5 * nsub + 3
14- lwork = nsub*nsub + 2 * nsub
15- lrwork = 3 * nsub*nsub + (4 + (int)log2(nsub) + 1) * nsub + 1
30+ integer :: lwork, isub, nsub2, info
31+ real(8) :: rwork(3*nsub-2)
32+ complex(8),allocatable :: work(:)
1633 !
34+ lwork = -1
35+ allocate(work(1))
36+ call zheev('V', 'U', nsub, ovlp, nsub, eval, work, lwork, rwork, info)
37+ lwork = nint(dble(work(1)))
38+ deallocate(work)
39+ allocate(work(lwork))
1740 call zheev('V', 'U', nsub, ovlp, nsub, eval, work, lwork, rwork, info)
41+ deallocate(work)
1842 !
1943 nsub2 = 0
2044 do isub = 1, nsub
@@ -28,26 +52,93 @@ contains
2852 hsub(1:nsub2,1:nsub2) = matmul(conjg(transpose(ovlp(1:nsub,1:nsub2))), &
2953 & matmul(hsub(1:nsub,1:nsub), ovlp(1:nsub,1:nsub2)))
3054 !
55+ lwork = -1
56+ allocate(work(1))
57+ call zheev('V', 'U', nsub2, hsub, nsub, eval, work, lwork, rwork, info)
58+ lwork = nint(dble(work(1)))
59+ deallocate(work)
60+ allocate(work(lwork))
3161 call zheev('V', 'U', nsub2, hsub, nsub, eval, work, lwork, rwork, info)
62+ deallocate(work)
3263 !
3364 hsub(1:nsub, 1:nsub2) = matmul(ovlp(1:nsub,1:nsub2), hsub(1:nsub2,1:nsub2))
3465 hsub(1:nsub, nsub2+1:nsub) = 0.0d0
3566 !
3667 end subroutine diag_ovrp
3768 !
38- subroutine lobpcg_main()
69+ subroutine lobpcg_main(linit,npw,nbnd,kvec,evec,eval)
70+ !
71+ use solver, only : electron_maxstep
72+ use hamiltonian, only : h_psi
73+ !
74+ logical,intent(in) :: linit
75+ integer,intent(in) :: npw, nbnd
76+ real(8),intent(in) :: kvec(3)
77+ real(8),intent(out) :: eval(nbnd)
78+ complex(8),intent(out) :: evec(npw,nbnd)
79+ !
80+ integer :: ii, ibnd, iter, nsub
81+ real(8) :: norm
82+ complex(8) :: wxp(npw,nbnd,3), hwxp(npw,nbnd,3), xp(npw,nbnd,2:3), &
83+ & hsub(nbnd,3,3*nbnd), ovlp(3*nbnd,3*nbnd), rotmat(nbnd,3,nbnd,2:3)
3984 !
4085 nsub = 3 * nbnd
86+ xp( 1:npw,1:nbnd,2:3) = 0.0d0
87+ wxp( 1:npw,1:nbnd,1:3) = 0.0d0
88+ hwxp(1:npw,1:nbnd,1:3) = 0.0d0
89+ !
90+ if(linit) call initialize(wxp(1:npw,1:nbnd,2))
4191 !
42- call initialize(wxp(1:npw,1:nbnd,1))
92+ call h_psi(kvec,wxp(1:npw,1:nbnd,2), hwxp(1:npw,1:nbnd,2))
4393 !
44- hwxp(1:npw,1:nbnd,1) = 0.0d0
94+ do ibnd = 1, nbnd
95+ eval(ibnd) = dble(dot_product(wxp(1:npw,ibnd,2), hwxp(1:npw,ibnd,2)))
96+ end do
97+ !
98+ do iter = 1, electron_maxstep
99+ !
100+ do ibnd = 1, nbnd
101+ !
102+ wxp(1:npw,ibnd,1) = hwxp(1:npw,ibnd,2) - eval(ibnd) * wxp(1:npw,ibnd,2)
103+ ! todo precondition
104+ norm = sqrt(dble(dot_product(wxp(1:npw,ibnd,1), wxp(1:npw,ibnd,1))))
105+ wxp(1:npw,ibnd,1) = wxp(1:npw,ibnd,1) / norm
106+ !
107+ end do
108+ !
109+ call h_psi(kvec,wxp(1:npw,1:nbnd,1), hwxp(1:npw,1:nbnd,1))
110+ !
111+ call zgemm("C", "N", 3*nbnd, 3*nbnd, npw, &
112+ & (1.0d0,0.0d0), wxp, npw, hwxp, npw, (0.0d0,0.0d0), hsub, 3*nbnd)
113+ call zgemm("C", "N", 3*nbnd, 3*nbnd, npw, &
114+ & (1.0d0,0.0d0), wxp, npw, wxp, npw, (0.0d0,0.0d0), ovlp, 3*nbnd)
115+ !
116+ call diag_ovrp(nsub,hsub,ovlp,eval)
117+ !
118+ rotmat(1:nbnd,1:3,1:nbnd,2) = hsub(1:nbnd,1:3,1:nbnd)
119+ rotmat(1:nbnd, 1,1:nbnd,3) = hsub(1:nbnd, 1,1:nbnd)
120+ rotmat(1:nbnd, 2,1:nbnd,3) = 0.0d0
121+ rotmat(1:nbnd, 3,1:nbnd,3) = hsub(1:nbnd, 3,1:nbnd)
122+ !
123+ call zgemm("C", "N", npw, 2*nbnd, 3*nbnd, &
124+ & (1.0d0,0.0d0), wxp, npw, rotmat, 3*nbnd, (0.0d0,0.0d0), xp, npw)
125+ wxp(1:npw,1:nbnd,2:3) = xp(1:npw,1:nbnd,2:3)
126+ call zgemm("C", "N", npw, 2*nbnd, 3*nbnd, &
127+ & (1.0d0,0.0d0), hwxp, npw, rotmat, 3*nbnd, (0.0d0,0.0d0), xp, npw)
128+ hwxp(1:npw,1:nbnd,2:3) = xp(1:npw,1:nbnd,2:3)
129+ !
130+ do ii = 2, 3
131+ do ibnd = 1, nbnd
132+ norm = sqrt(dble(dot_product(wxp(1:npw,ibnd,ii), wxp(1:npw,ibnd,ii))))
133+ wxp( 1:npw,ibnd,ii) = wxp(1:npw,ibnd,ii) / norm
134+ hwxp(1:npw,ibnd,ii) = hwxp(1:npw,ibnd,ii) / norm
135+ end do
136+ end do
137+ !
138+ end do
45139 !
46- call h_psi(wxp(1:npw,1:nbnd,1), hwxp(1:npw,1:nbnd,1))
140+ evec(1:npw,1:nbnd) = wxp(1:npw,1:nbnd,2)
47141 !
48-
49-
50-
51142 end subroutine lobpcg_main
52143 !
53144 end module lobpcg
--- a/src/make.depend
+++ b/src/make.depend
@@ -1,5 +1,10 @@
11 atm_spec.o : atm_spec.F90
2+diag_direct.o : diag_direct.F90 gvec.o
3+gvec.o : gvec.F90 atm_spec.o
4+hamiltonian.o : hamiltonian.F90 solver.o gvec.o
5+io_vloc.o : io_vloc.F90 gvec.o atm_spec.o
26 k_point.o : k_point.F90
3-pwdft.o : pwdft.F90 stdin.o
7+lobpcg.o : lobpcg.F90 hamiltonian.o solver.o
8+pwdft.o : pwdft.F90 io_vloc.o gvec.o stdin.o solver.o
49 solver.o : solver.F90
5-stdin.o : stdin.F90 atm_spec.o k_point.o solver.o
10+stdin.o : stdin.F90 gvec.o atm_spec.o k_point.o solver.o
--- a/src/pwdft.F90
+++ b/src/pwdft.F90
@@ -1,9 +1,15 @@
11 program pwdft
22 !
3+ use solver, only : calculation
34 use stdin, only : read_stdin
5+ use gvec, only : rfft
6+ use io_vloc, only : vloc, read_vloc
47 !
58 implicit none
69 !
710 call read_stdin()
811 !
12+ allocate(Vloc(rfft%npw3(1),rfft%npw3(2),rfft%npw3(3)))
13+ if(calculation /= "scf") call read_vloc()
14+ !
915 end program pwdft
--- a/src/solver.F90
+++ b/src/solver.F90
@@ -2,8 +2,9 @@ module solver
22 !
33 implicit none
44 !
5- character(256) :: calculation
6- integer :: nbnd
7- real(8) :: ecutwfc
5+ character(256),save :: calculation
6+ integer,save :: nbnd
7+ integer,save :: electron_maxstep
8+ real(8),save :: conv_thr
89 !
910 end module solver
--- a/src/stdin.F90
+++ b/src/stdin.F90
@@ -6,29 +6,41 @@ contains
66 !
77 subroutine read_stdin()
88 !
9- use solver, only : calculation, ecutwfc, nbnd
9+ use solver, only : calculation, nbnd, electron_maxstep, conv_thr
1010 use k_point, only : nk, kvec, kgrd
11- use atm_spec, only : atm, spec, avec, nat, ntyp
11+ use atm_spec, only : atm, spec, avec, nat, ntyp, bvec
12+ use gvec, only : rfft, wfft, setup_gvec
1213 !
13- integer :: iat, jtyp, ii, ik, jj
14+ integer :: iat, jtyp, ii, ik, jj, lwork = 3, ipiv(3), info
15+ real(8) :: work(3), pi = acos(-1.0d0), ecutrho, ecutwfc
1416 character(256) :: key
1517 !
1618 namelist/control/ calculation
17- namelist/system/ nat, ntyp, ecutwfc, nbnd
19+ namelist/system/ nat, ntyp, ecutwfc, ecutrho, nbnd
20+ namelist/electrons/ electron_maxstep, conv_thr
1821 !
1922 read(*,control)
2023 !
21- write(*,*) " calculation : ", trim(calculation)
24+ write(*,*) " calculation : ", trim(calculation)
2225 !
2326 read(*,system)
2427 !
25- write(*,*) " Number of atoms : ", nat
26- write(*,*) " Number of species : ", ntyp
27- write(*,*) " Plane-wave cutoff [Ry] : ", ecutwfc
28- write(*,*) " Number of bands : ", nbnd
28+ write(*,*) " Number of atoms : ", nat
29+ write(*,*) " Number of species : ", ntyp
30+ write(*,*) " Plane-wave cutoff (wfc) [Ry] : ", ecutwfc
31+ write(*,*) " Plane-wave cutoff (rho,V) [Ry] : ", ecutrho
32+ write(*,*) " Number of bands : ", nbnd
2933 !
3034 allocate(atm(nat), spec(ntyp))
3135 !
36+ electron_maxstep = 100
37+ conv_thr = 1.0d-10
38+ !
39+ read(*,electrons)
40+ !
41+ write(*,*) " Max iteration : ", electron_maxstep
42+ write(*,*) " Convergence threshold : ", conv_thr
43+ !
3244 do jj = 1, 4
3345 !
3446 read(*,*) key
@@ -63,12 +75,21 @@ contains
6375 else if(key == "CELL_PARAMETERS") then
6476 !
6577 read(*,*) avec(1:3,1:3)
66- write(*,*) " Cell parameters:"
78+ avec(1:3,1:3) = avec(1:3,1:3) / 0.529177249d0
79+ write(*,*) " Cell parameters [Bohr] :"
6780 !
6881 do ii = 1, 3
6982 write(*,*) avec(1:3,ii)
7083 end do
7184 !
85+ ! Reciprocal lattice vectors
86+ !
87+ bvec(1:3,1:3) = transpose(avec(1:3,1:3))
88+ call dgetrf(3, 3, bvec, 3, ipiv, info)
89+ call dgetri(3, bvec, 3, ipiv, work, lwork, info)
90+ if(info /= 0) stop 'read_stdin : inverce of avec.'
91+ bvec(1:3,1:3) = 2.0d0 * pi * bvec(1:3,1:3)
92+ !
7293 end if
7394 !
7495 end do
@@ -88,6 +109,11 @@ contains
88109 !
89110 end do
90111 !
112+ write(*,*) " FFT and G-vector for wfc"
113+ call setup_gvec(wfft, ecutwfc)
114+ write(*,*) " FFT and G-vector for rho and V"
115+ call setup_gvec(rfft, ecutrho)
116+ !
91117 end subroutine read_stdin
92118 !
93119 end module stdin
Show on old repository browser