Loading [MathJax]/extensions/tex2jax.js
pwdft  0.1
PW-DFT code for education
All Classes Namespaces Files Functions Variables
rho_v Module Reference

Functions/Subroutines

subroutine init_rho_v ()
 Initialize rho, Vps, Vks. More...
 
subroutine generate_vps ()
 Pseudopotential Vps is computed. More...
 
subroutine hartree_pot (Vloc)
 Add Hartree potential. More...
 
subroutine xc_pot (Vloc)
 Add XC potential (LDA) More...
 

Variables

real(8), dimension(:), allocatable vks
 (g_rhnr) Kohn-Sham potential [Htr] More...
 
real(8), dimension(:), allocatable vps
 (g_rhnr) Pseudopotential [Htr] More...
 
real(8), dimension(:), allocatable rho
 (g_rhnr) Charge density More...
 

Function/Subroutine Documentation

◆ init_rho_v()

subroutine rho_v::init_rho_v

Initialize rho, Vps, Vks.

Definition at line 37 of file rho_v.F90.

37  !
38  use gvec, only : g_rh
39  use atm_spec, only : vcell, nelec
40  !
41  allocate(vps(g_rh%nr), vks(g_rh%nr), rho(g_rh%nr))
42  !
43  rho(1:g_rh%nr) = nelec / vcell
44  !
45  call generate_vps()
46  vks(1:g_rh%nr) = vps(1:g_rh%nr)
47  call hartree_pot(vks)
48  call xc_pot(vks)
49  !

References gvec::g_rh, generate_vps(), hartree_pot(), atm_spec::nelec, rho, atm_spec::vcell, vks, vps, and xc_pot().

Referenced by pwdft().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ generate_vps()

subroutine rho_v::generate_vps

Pseudopotential Vps is computed.

Definition at line 55 of file rho_v.F90.

55  !
56  use constant, only : pi
57  use gvec, only : g_rh
58  use atm_spec, only : bvec, spec, atm, nat, ntyp, vcell
59  use fftw_wrapper, only : fft_g2r
60  !
61  integer :: jtyp, iat, nat2, ipw
62  real(8) :: prod(nat), pos2(3,nat), gv(3), ctail, glen, dr, &
63  & frmfac(g_rh%nr)
64  complex(8) :: strfac(g_rh%nr), VpsG(g_rh%nr)
65  real(8),allocatable :: sinc(:)
66  !
67  vpsg(1:g_rh%nr) = 0.0d0
68  !
69  do jtyp = 1, ntyp
70  !
71  ! Structure factor
72  !
73  nat2 = 0
74  do iat = 1, nat
75  if(atm(iat)%ityp == jtyp) then
76  nat2 = nat2 + 1
77  pos2(1:3,nat2) = atm(iat)%pos(1:3)
78  end if
79  end do
80  !
81  do ipw = 1, g_rh%nr
82  prod(1:nat2) = matmul(dble(g_rh%mill(1:3,ipw)),pos2(1:3,1:nat2))
83  strfac(ipw) = sum(exp(cmplx(0.0d0, - 2.0d0 * pi * prod(1:nat2), 8)))
84  end do
85  strfac(1:g_rh%nr) = strfac(1:g_rh%nr) / vcell
86  !
87  ! Form factor
88  !
89  allocate(sinc(spec(jtyp)%mmax))
90  dr = spec(jtyp)%psr(spec(jtyp)%mmax) / dble(spec(jtyp)%mmax - 1)
91  !
92  ! G = 0 : Compensation to the Hartree potential
93  !
94  ctail = 0.5d0 * spec(jtyp)%Zion * spec(jtyp)%psr(spec(jtyp)%mmax)**2
95  !
96  sinc(1:spec(jtyp)%mmax) = dr
97  sinc(1 ) = dr * 0.5d0
98  sinc(spec(jtyp)%mmax) = dr * 0.5d0
99  !
100  frmfac(1) = sum( sinc(1:spec(jtyp)%mmax) &
101  & * spec(jtyp)%psr(1:spec(jtyp)%mmax)**2 &
102  & * spec(jtyp)%psV(1:spec(jtyp)%mmax)) &
103  & + ctail
104  frmfac(1) = 4.0d0 * pi * frmfac(1)
105  !
106  ! G /= 0
107  !
108  do ipw = 2, g_rh%nr
109  !
110  gv(1:3) = matmul(bvec(1:3,1:3), dble(g_rh%mill(1:3,ipw)))
111  glen = sqrt(dot_product(gv(1:3),gv(1:3)))
112  sinc(1:spec(jtyp)%mmax) = glen*spec(jtyp)%psr(1:spec(jtyp)%mmax)
113  !
114  ctail = - spec(jtyp)%Zion * cos(sinc(spec(jtyp)%mmax)) / glen**2
115  !
116  sinc(2:spec(jtyp)%mmax) = sin(sinc(2:spec(jtyp)%mmax)) &
117  & / sinc( 2:spec(jtyp)%mmax) &
118  & * dr
119  sinc(1 ) = dr * 0.5d0
120  sinc(spec(jtyp)%mmax) = sinc(spec(jtyp)%mmax) * 0.5d0
121  !
122  frmfac(ipw) = sum( sinc(1:spec(jtyp)%mmax) &
123  & * spec(jtyp)%psr(1:spec(jtyp)%mmax)**2 &
124  & * spec(jtyp)%psV(1:spec(jtyp)%mmax)) &
125  & + ctail
126  frmfac(ipw) = 4.0d0 * pi * frmfac(ipw)
127  !
128  end do
129  !
130  deallocate(sinc)
131  !
132  vpsg(1:g_rh%nr) = vpsg(1:g_rh%nr) + strfac(1:g_rh%nr) * frmfac(1:g_rh%nr)
133  !
134  end do
135  !
136  call fft_g2r(vpsg, vps)
137  !

References atm_spec::atm, atm_spec::bvec, fftw_wrapper::fft_g2r(), gvec::g_rh, atm_spec::nat, atm_spec::ntyp, constant::pi, atm_spec::spec, atm_spec::vcell, and vps.

Referenced by init_rho_v().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ hartree_pot()

subroutine rho_v::hartree_pot ( real(8), dimension(g_rh%nr), intent(inout)  Vloc)

Add Hartree potential.

Definition at line 143 of file rho_v.F90.

143  !
144  use constant, only : pi
145  use gvec, only : g_rh
146  use atm_spec, only : bvec
147  use fftw_wrapper, only : fft_g2r, fft_r2g
148  !
149  real(8),intent(inout) :: Vloc(g_rh%nr)
150  !
151  integer :: ir
152  real(8) :: g3(3), VH(g_rh%nr)
153  complex(8) :: rhog(g_rh%nr)
154  !
155  call fft_r2g(rho, rhog)
156  !
157  ! G = 0 : Compensation to the ionic potential
158  !
159  rhog(1) = 0.0d0
160  !
161  do ir = 2, g_rh%nr
162  g3(1:3) = matmul(bvec(1:3,1:3), dble(g_rh%mill(1:3,ir)))
163  rhog(ir) = 4.0d0 * pi * rhog(ir) / dot_product(g3(1:3),g3(1:3))
164  end do
165  !
166  call fft_g2r(rhog, vh)
167  !
168  vloc(1:g_rh%nr) = vloc(1:g_rh%nr) + vh(1:g_rh%nr)
169  !

References atm_spec::bvec, fftw_wrapper::fft_g2r(), fftw_wrapper::fft_r2g(), gvec::g_rh, constant::pi, and rho.

Referenced by init_rho_v(), and scf::kohn_sham_eq().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ xc_pot()

subroutine rho_v::xc_pot ( real(8), dimension(g_rh%nr), intent(inout)  Vloc)

Add XC potential (LDA)

Definition at line 175 of file rho_v.F90.

175  !
176  use constant, only : pi
177  use gvec, only : g_rh
178  !
179  real(8),intent(inout) :: Vloc(g_rh%nr)
180  !
181  integer :: ir
182  real(8) :: rs, exc, Vxc
183  !
184  do ir = 1, g_rh%nr
185  !
186  if(rho(ir) > 1.0d-10) then
187  !
188  rs = (3.0d0 / (4.0d0 * pi * rho(ir)))**(1.0d0/3.0d0)
189  !
190  if(rs < 1.0d0) then
191  exc = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
192  & - 0.0480d0 + 0.031d0*log(rs) - 0.0116d0*rs + 0.0020d0*rs*log(rs)
193  !
194  vxc = exc -1.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs**2 &
195  & -0.0096d0*rs + 0.031d0 + 0.002d0*rs*log(rs)
196  else
197  exc = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
198  & - 0.1423d0 / (1.0d0 + 1.0529d0*sqrt(rs) + 0.3334d0*rs)
199  !
200  vxc = exc -1.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs**2 &
201  & - 0.0474333d0*(0.3334d0*rs + 0.52645*sqrt(rs)) &
202  & /(1.0d0 + 1.0529d0*sqrt(rs) + 0.3334d0*rs)**2
203  end if
204  !
205  vloc(ir) = vloc(ir) + vxc
206  !
207  end if
208  !
209  end do
210  !

References gvec::g_rh, constant::pi, and rho.

Referenced by init_rho_v(), and scf::kohn_sham_eq().

Here is the caller graph for this function:

Variable Documentation

◆ vks

real(8), dimension(:), allocatable rho_v::vks

(g_rhnr) Kohn-Sham potential [Htr]

Definition at line 27 of file rho_v.F90.

27  real(8),allocatable :: &
28  & Vks(:), & !< (g_rh%nr) Kohn-Sham potential [Htr]
29  & Vps(:), & !< (g_rh%nr) Pseudopotential [Htr]
30  & rho(:) !< (g_rh%nr) Charge density

Referenced by diag_direct::direct(), hamiltonian::h_psi(), init_rho_v(), scf::kohn_sham_eq(), pwdft(), and scf::scf_loop().

◆ vps

real(8), dimension(:), allocatable rho_v::vps

(g_rhnr) Pseudopotential [Htr]

Definition at line 27 of file rho_v.F90.

Referenced by energy::atomic(), generate_vps(), init_rho_v(), and scf::kohn_sham_eq().

◆ rho

real(8), dimension(:), allocatable rho_v::rho

(g_rhnr) Charge density

Definition at line 27 of file rho_v.F90.

Referenced by energy::atomic(), energy::hartree(), hartree_pot(), init_rho_v(), k_point::ksum_rho(), energy::xc(), and xc_pot().

atm_spec::spec
type(spec_t), dimension(:), allocatable, save spec
(ntyp) Species
Definition: atm_spec.F90:57
atm_spec::vcell
real(8), save vcell
Unit cell volume [Bohr^3].
Definition: atm_spec.F90:50
fftw_wrapper::fft_r2g
subroutine fft_r2g(VlocR, VlocG)
v(r) e^{-iGr} -> V(G)
Definition: fftw_wrapper.F90:75
constant
Definition: constant.F90:23
rho_v::rho
real(8), dimension(:), allocatable rho
(g_rhnr) Charge density
Definition: rho_v.F90:27
fftw_wrapper
Definition: fftw_wrapper.F90:23
atm_spec::ntyp
integer, save ntyp
Number of species (elements)
Definition: atm_spec.F90:47
atm_spec::atm
type(atm_t), dimension(:), allocatable, save atm
(nat) Atom
Definition: atm_spec.F90:55
gvec
Definition: gvec.F90:23
atm_spec::nat
integer, save nat
Number of atoms.
Definition: atm_spec.F90:47
fftw_wrapper::fft_g2r
subroutine fft_g2r(VlocG, VlocR)
V(G) e^{iGr} -> v(r)
Definition: fftw_wrapper.F90:90
atm_spec
Definition: atm_spec.F90:23