File: atom.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (189 lines) | stat: -rw-r--r-- 5,653 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: atom
! !INTERFACE:
subroutine atom(sol,ptnucl,zn,nst,n,l,k,occ,xctype,xcgrad,nr,r,eval,rho,vr,rwf)
! !USES:
use modxcifc
! !INPUT/OUTPUT PARAMETERS:
!   sol    : speed of light in atomic units (in,real)
!   ptnucl : .true. if the nucleus is a point particle (in,logical)
!   zn     : nuclear charge (in,real)
!   nst    : number of states to solve for (in,integer)
!   n      : priciple quantum number of each state (in,integer(nst))
!   l      : quantum number l of each state (in,integer(nst))
!   k      : quantum number k (l or l+1) of each state (in,integer(nst))
!   occ    : occupancy of each state (inout,real(nst))
!   xctype : exchange-correlation type (in,integer(3))
!   xcgrad : 1 for GGA functional, 0 otherwise (in,integer)
!   nr     : number of radial mesh points (in,integer)
!   r      : radial mesh (in,real(nr))
!   eval   : eigenvalue without rest-mass energy for each state (out,real(nst))
!   rho    : charge density (out,real(nr))
!   vr     : self-constistent potential (out,real(nr))
!   rwf    : major and minor components of radial wavefunctions for each state
!            (out,real(nr,2,nst))
! !DESCRIPTION:
!   Solves the Dirac-Kohn-Sham equations for an atom using the
!   exchange-correlation functional {\tt xctype} and returns the self-consistent
!   radial wavefunctions, eigenvalues, charge densities and potentials. Requires
!   the exchange-correlation interface routine {\tt xcifc}.
!
! !REVISION HISTORY:
!   Created September 2002 (JKD)
!   Fixed s.c. convergence problem, October 2003 (JKD)
!   Added support for GGA functionals, June 2006 (JKD)
!
!EOP
!BOC
implicit none
! arguments
real(8), intent(in) :: sol
logical, intent(in) :: ptnucl
real(8), intent(in) :: zn
integer, intent(in) :: nst
integer, intent(in) :: n(nst),l(nst),k(nst)
real(8), intent(inout) :: occ(nst)
integer, intent(in) :: xctype(3),xcgrad
integer, intent(in) :: nr
real(8), intent(in) :: r(nr)
real(8), intent(out) :: eval(nst)
real(8), intent(out) :: rho(nr),vr(nr)
real(8), intent(out) :: rwf(nr,2,nst)
! local variables
integer, parameter :: maxscl=200
integer ir,ist,iscl
real(8), parameter :: fourpi=12.566370614359172954d0
! potential convergence tolerance
real(8), parameter :: eps=1.d-6
real(8) sum,dv,dvp,ze,beta,t1
! allocatable arrays
real(8), allocatable :: vn(:),vh(:),ex(:),ec(:),vx(:),vc(:),vrp(:)
real(8), allocatable :: ri(:),fr1(:),fr2(:),gr1(:),gr2(:)
real(8), allocatable :: grho(:),g2rho(:),g3rho(:)
if (nst.le.0) then
  write(*,*)
  write(*,'("Error(atom): invalid nst : ",I8)') nst
  write(*,*)
  stop
end if
! allocate local arrays
allocate(vn(nr),vh(nr),ex(nr),ec(nr),vx(nr),vc(nr),vrp(nr))
allocate(ri(nr),fr1(nr),fr2(nr),gr1(nr),gr2(nr))
if (xcgrad.eq.1) then
  allocate(grho(nr),g2rho(nr),g3rho(nr))
end if
! find total electronic charge
ze=0.d0
do ist=1,nst
  ze=ze+occ(ist)
end do
! set up nuclear potential
call potnucl(ptnucl,nr,r,zn,vn)
do ir=1,nr
  ri(ir)=1.d0/r(ir)
! initialise the Kohn-Sham potential to the nuclear potential
  vr(ir)=vn(ir)
end do
dvp=0.d0
vrp(:)=0.d0
! initialise mixing parameter
beta=0.5d0
! initialise eigenvalues to relativistic values (minus the rest mass energy)
do ist=1,nst
  t1=sqrt(dble(k(ist)**2)-(zn/sol)**2)
  t1=(dble(n(ist)-abs(k(ist)))+t1)**2
  t1=1.d0+((zn/sol)**2)/t1
  eval(ist)=sol**2/sqrt(t1)-sol**2
end do
! start of self-consistent loop
do iscl=1,maxscl
! solve the Dirac equation for each state
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO
  do ist=1,nst
    call rdirac(sol,n(ist),l(ist),k(ist),nr,r,vr,eval(ist),rwf(:,1,ist), &
     rwf(:,2,ist))
  end do
!$OMP END DO
!$OMP END PARALLEL
! compute the charge density
  do ir=1,nr
    sum=0.d0
    do ist=1,nst
      sum=sum+occ(ist)*(rwf(ir,1,ist)**2+rwf(ir,2,ist)**2)
    end do
    fr1(ir)=sum
    fr2(ir)=sum*ri(ir)
    rho(ir)=(1.d0/fourpi)*sum*ri(ir)**2
  end do
  call splintp(nr,r,fr1,gr1)
  call splintp(nr,r,fr2,gr2)
! find the Hartree potential
  t1=gr2(nr)
  do ir=1,nr
    vh(ir)=gr1(ir)*ri(ir)+t1-gr2(ir)
  end do
! normalise charge density and potential
  t1=ze/gr1(nr)
  rho(:)=t1*rho(:)
  vh(:)=t1*vh(:)
! compute the exchange-correlation energy and potential
  if (xcgrad.eq.1) then
! GGA functional
! |grad rho|
    call fderiv(1,nr,r,rho,grho)
! grad^2 rho
    call fderiv(2,nr,r,rho,g2rho)
    do ir=1,nr
      g2rho(ir)=g2rho(ir)+2.d0*ri(ir)*grho(ir)
    end do
! approximate (grad rho).(grad |grad rho|)
    do ir=1,nr
      g3rho(ir)=grho(ir)*g2rho(ir)
    end do
    call xcifc(xctype,n=nr,rho=rho,grho=grho,g2rho=g2rho,g3rho=g3rho,ex=ex, &
     ec=ec,vx=vx,vc=vc)
  else
! LDA functional
    call xcifc(xctype,n=nr,rho=rho,ex=ex,ec=ec,vx=vx,vc=vc)
  end if
! self-consistent potential
  vr(:)=vh(:)+vx(:)+vc(:)
! determine change in potential
  sum=0.d0
  do ir=1,nr
    sum=sum+(vr(ir)-vrp(ir))**2
  end do
  dv=sqrt(sum)/dble(nr)
  if (iscl.gt.2) then
! reduce beta if change in potential is diverging
    if (dv.gt.dvp) beta=beta*0.8d0
    beta=max(beta,0.01d0)
  end if
  dvp=dv
  do ir=1,nr
! mix old and new potentials
    vr(ir)=(1.d0-beta)*vrp(ir)+beta*vr(ir)
    vrp(ir)=vr(ir)
! add nuclear potential
    vr(ir)=vr(ir)+vn(ir)
  end do
! check for convergence
  if ((iscl.gt.2).and.(dv.lt.eps)) goto 10
! end self-consistent loop
end do
write(*,*)
write(*,'("Warning(atom): maximum iterations exceeded")')
10 continue
deallocate(vn,vh,ex,ec,vx,vc,vrp)
deallocate(ri,fr1,fr2,gr1,gr2)
if (xcgrad.eq.1) deallocate(grho,g2rho,g3rho)
return
end subroutine
!EOC