1#if defined HAVE_CONFIG_H
2#include "config.h"
3#endif
4
5PROGRAM gridXCtest2
6
7  ! Compares potential and the numerical derivative of the energy
8  ! calculated by atomXC. J.M.Soler. Sept.2009
9
10  ! Used module procedures
11  USE gridXC, only: atomXC => gridxc_atomXC
12  USE gridXC, only: setXC  => gridxc_setXC
13  USE gridXC, only: gridxc_init
14
15  ! Used module parameters
16  USE gridXC, only: dp
17
18! Used MPI procedures and types
19! Note that current MPI modules do not include the actual
20! routines
21#ifdef HAVE_MPI
22  USE mpi
23#endif
24
25  implicit none
26
27  ! Tester parameters
28  integer, parameter:: irel  =  0 ! Relativistic? 0=>no, 1=>yes
29  integer, parameter:: nSpin =  2 ! Number of spin components
30  integer, parameter:: nfTot = 19 ! Number of functionals
31  integer, parameter:: nr    = 101 ! Number of radial points
32  integer, parameter:: n1cut =  8 ! Cutoff parameter
33  integer, parameter:: n2cut =  2 ! Cutoff parameter:
34                                  !    fCut(r)=(1-(r/rMax)**n1cut)**n2cut
35  real(dp),parameter:: dWidth = 2._dp ! Width of density distribution, in Bohr
36  real(dp),parameter:: Qtot = 10._dp  ! Integral of density distribution
37  real(dp),parameter:: spinPol= 2._dp ! Integral of densUp - densDown
38  real(dp),parameter:: rMax = 20._dp  ! Cutoff radius, in Bohr
39  real(dp),parameter:: deltaDens = 1.e-8_dp  ! Finite diff. change
40  real(dp),parameter:: densMin  = 1.e-9_dp  ! Min. density to proceed
41
42  ! List of functionals to be tested
43!  integer, parameter:: nf = nfTot-4   ! Number of tested functionals
44!  integer:: indexf(nf) = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,  18/)
45!                         ! Indexes from list below (only one VDW allowed)
46
47  ! Same to test a single functional
48  integer, parameter:: nf = 1        ! Number of tested functionals
49  integer:: indexf(nf) = (/18/)      ! Indexes from list below
50
51  ! All functionals available
52  !                  1,           2,          3,           4,
53  !                  5,           6,          7,           8,
54  !                  9,          10,         11,          12,
55  !                 13,          14,         15,          16,
56  !                 17,          18,         19
57  character(len=3):: &
58    func(nfTot)=(/'LDA',       'LDA',       'GGA',       'GGA', &
59                  'GGA',       'GGA',       'GGA',       'GGA', &
60                  'GGA',       'GGA',       'GGA',       'GGA', &
61                  'GGA',       'GGA',       'VDW',       'VDW', &
62                  'VDW',       'VDW',       'VDW'       /)
63  character(len=10):: &
64    auth(nfTot)=(/'PZ        ','PW92      ','PW91      ','PBE       ', &
65                  'RPBE      ','revPBE    ','LYP       ','WC        ', &
66                  'PBEJsJrLO ','PBEJsJrHEG','PBEGcGxLO ','PBEGcGxHEG', &
67                  'PBESOL    ','AM05      ','DRSLL     ','LMKLL     ', &
68                  'C09       ','BH        ','VV        '/)
69
70  ! Tester variables and arrays
71  integer :: iDelta, ir, irmax, ismax, iSpin, one, two
72  real(dp):: avgDiffVxc, dDensdr, dens(nr,nSpin), dens0(nr,nSpin), &
73             d0tot, d0(nSpin), dEdDens, dDens, diffVxc, &
74             Dc, Dc0, dr, dVol, Dx, Dx0, Ec, Ec0, Ex, Ex0, &
75             kf, kg, maxDiffVxc, pi, r, rMesh(nr), &
76             Vxc(nr,nSpin), Vxc0(nr,nSpin), wc(nfTot), wr, wx(nfTot)
77
78#ifdef HAVE_MPI
79  ! Initialize MPI, even though this test is intended to be run serially
80  integer:: MPIerror, myNode, nNodes
81  call MPI_Init( MPIerror )
82  call MPI_Comm_Rank( MPI_Comm_World, myNode, MPIerror )
83  call MPI_Comm_Size( MPI_Comm_World, nNodes, MPIerror )
84#endif
85
86  ! Initialize cocktail XC functional with all tested functionals
87  wx = 1._dp / nf
88  wc = 1._dp / nf
89  call setXC( nf, func(indexf), auth(indexf), wx(indexf), wc(indexf) )
90
91  ! Find radial mesh points and gaussian density
92  pi = acos(-1._dp)
93  d0tot = Qtot / (2*pi*dWidth**2)**1.5_dp    ! Total density at origin
94  if (nSpin==1) then
95    d0(1) = d0tot
96  else
97    one = 1   ! A silly thing to satisfy the compiler when nSpin=1
98    two = 2
99    d0(one) = d0tot * (Qtot + spinPol) / Qtot / 2 ! Spin up density at origin
100    d0(two) = d0tot * (Qtot - spinPol) / Qtot / 2 ! Spin down density at origin
101  end if
102  dr = rmax / (nr-1)                      ! Interval between radial points
103  do ir = 1,nr
104    rMesh(ir) = dr * (ir-1)               ! Radial point values
105    dens0(ir,:) = DensOfR( d0(:), rMesh(ir) )
106  end do
107
108  ! Find exchange and correlation energy and potential from radial density
109  call atomXC( irel, nr, nr, rMesh, nSpin, dens0, Ex0, Ec0, Dx0, Dc0, Vxc0 )
110
111  ! Print parameters
112  print'(/,a,10(a3,4x))', 'funcs= ', func(indexf)
113  print  '(a,10(a6,1x))', 'auths= ', auth(indexf)
114  print'(a,2f12.6)', 'dr, rMax = ', dr, rMax
115
116  ! Calculate finite-difference derivatives
117  open( unit=44, file='test2.Vxc' )
118  avgDiffVxc = 0
119  maxDiffVxc = 0
120  do iSpin = 1,nSpin
121    do ir = 2,nr
122      if (dens0(ir,iSpin)<densMin) cycle ! Do nothing if dens=0
123      dVol = 4*pi*rMesh(ir)**2 * dr
124      dDens = min( deltaDens, dens0(ir,iSpin)/100 )
125      dEdDens = 0
126      do iDelta = -1,1,2
127        dens = dens0
128        dens(ir,iSpin) = dens0(ir,iSpin) + iDelta*dDens
129        call atomXC( irel, nr, nr, rMesh, nSpin, dens, Ex, Ec, Dx, Dc, Vxc )
130        dEdDens = dEdDens + iDelta * (Ex+Ec) / (2*dDens) / dVol
131      end do ! iDelta
132      diffVxc = Vxc0(ir,iSpin) - dEdDens
133      avgDiffVxc = avgDiffVxc + diffVxc**2
134      if (abs(diffVxc) > maxDiffVxc) then
135        maxDiffVxc = abs(diffVxc)
136        irMax = ir
137        isMax = iSpin
138      end if
139      kf = (3*pi**2*sum(dens0(ir,:)))**(1.0_dp/3)
140      dDensdr = sum(dens0(ir+1,:)-dens0(ir-1,:))/(rMesh(ir+1)-rMesh(ir-1))
141      kg = abs(dDensdr)/sum(dens0(ir,:))
142      if (ir==2) then
143        print'(a5,a9,4a15)','iSpin','r','dens','Vxc','dExc/dDens','diff'
144        print'(i5,f9.3,4f15.9)', &
145          ispin, rMesh(1), dens0(1,iSpin), Vxc0(1,iSpin)
146!        print'(a5,a9,6a15)','iSpin','r','dens','kf','kg','Vxc', &
147!                            'dExc/dDens','diff'
148!        print'(i5,f9.3,6f15.9)', &
149!          ispin, rMesh(1), dens0(1,iSpin), kf, kg, Vxc0(1,iSpin)
150      end if
151      print'(i5,f9.3,4f15.9)', &
152        ispin, rMesh(ir), dens0(ir,iSpin), Vxc0(ir,iSpin), dEdDens, diffVxc
153!      print'(i5,f9.3,6f15.9)', &
154!        ispin, rMesh(ir), dens0(ir,iSpin), kf, kg, &
155!        Vxc0(ir,iSpin), dEdDens, diffVxc
156      write(44,'(f9.3,4f15.9)') &
157        rMesh(ir), dens0(ir,iSpin), Vxc0(ir,iSpin), dEdDens, diffVxc
158    end do ! ir
159  end do ! iSpin
160  avgDiffVxc = sqrt( avgDiffVxc / nSpin / nr )
161  print'(a,2f15.9)', 'avgDiffVxc, maxDiffVxc = ', avgDiffVxc, maxDiffVxc
162!  print'(a,2i6)', 'irMax, iSpinMax = ', irmax, ismax
163  close( unit=44 )
164
165! Finalize MPI
166#ifdef HAVE_MPI
167  call MPI_Finalize( MPIerror )
168#endif
169
170CONTAINS
171
172FUNCTION DensOfR( d0, r )
173
174  ! Returns a radial density distribution
175
176  implicit none
177  real(dp),intent(in):: d0(nSpin)  ! Density at center of charge distribution
178  real(dp),intent(in):: r          ! Distance to center of charge distribution
179  real(dp)           :: DensOfR(nSpin)  ! Electron density
180
181  ! Use a simple gaussian distribution
182  DensOfR = d0 * exp(-r**2/2/dWidth**2)
183
184  ! Impose a smooth radial cutoff
185  DensOfR = DensOfR * ( 1 - (r/rMax)**n1cut )**n2cut
186
187END FUNCTION DensOfR
188
189END PROGRAM gridXCtest2
190
191