1PROGRAM siestaXCtest4
2
3  ! Compares the energy and potential calculated by atomXC and cellXC.
4  ! J.M.Soler. Sept.2009
5
6  ! Used module procedures
7  USE siestaXC, only: atomXC
8  USE siestaXC, only: cellXC
9  USE siestaXC, only: setXC
10
11  ! Used module parameters
12  USE siestaXC, only: dp
13  USE siestaXC, only: gp => grid_p
14
15! Used MPI types
16#ifdef MPI
17  USE mpi_siesta, only: MPI_Double_Precision
18  USE mpi_siesta, only: MPI_Max
19  USE mpi_siesta, only: MPI_Sum
20  USE mpi_siesta, only: MPI_Comm_World
21#endif
22
23  implicit none
24
25  ! Tester parameters
26  integer, parameter:: irel  =  0 ! Relativistic? 0=>no, 1=>yes
27  integer, parameter:: nSpin =  2 ! Number of spin components
28  integer, parameter:: nfTot = 19 ! Number of functionals
29  integer, parameter:: nr = 501   ! Number of radial points
30  integer, parameter:: nx = 60    ! Number of grid points per lattice vector
31  integer, parameter:: n1cut = 8  ! Cutoff parameter
32  integer, parameter:: n2cut = 2  ! Cutoff parameter:
33                                  !    fCut(r)=(1-(r/rMax)**n1cut)**n2cut
34  real(dp),parameter:: dWidth = 2._dp ! Width of density distribution, in Bohr
35  real(dp),parameter:: Qtot = 10._dp  ! Integral of density distribution
36  real(dp),parameter:: spinPol= 2._dp ! Integral of densUp - densDown
37  real(dp),parameter:: rMax = 12._dp  ! Cutoff radius, in Bohr
38  real(dp),parameter:: rBuff = 3._dp  ! Radial buffer of zero density, in Bohr
39  real(dp),parameter:: densMin  = 1.e-9_dp  ! Min. density to proceed
40
41  ! List of functionals to be tested
42!  integer, parameter:: nf = nfTot-4   ! Number of tested functionals
43!  integer:: indexf(nf) = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,  18/)
44!                         ! Indexes from list below (only one VDW allowed)
45
46  ! Same to test a single functional
47  integer, parameter:: nf = 1        ! Number of tested functionals
48  integer:: indexf(nf) = (/18/)      ! Indexes from list below
49
50  ! All functionals available
51  !                  1,           2,          3,           4,
52  !                  5,           6,          7,           8,
53  !                  9,          10,         11,          12,
54  !                 13,          14,         15,          16,
55  !                 17,          18,         19
56  character(len=3):: &
57    func(nfTot)=(/'LDA',       'LDA',       'GGA',       'GGA', &
58                  'GGA',       'GGA',       'GGA',       'GGA', &
59                  'GGA',       'GGA',       'GGA',       'GGA', &
60                  'GGA',       'GGA',       'VDW',       'VDW', &
61                  'VDW',       'VDW',       'VDW'       /)
62  character(len=10):: &
63    auth(nfTot)=(/'PZ        ','PW92      ','PW91      ','PBE       ', &
64                  'RPBE      ','revPBE    ','LYP       ','WC        ', &
65                  'PBEJsJrLO ','PBEJsJrHEG','PBEGcGxLO ','PBEGcGxHEG', &
66                  'PBESOL    ','AM05      ','DRSLL     ','LMKLL     ', &
67                  'C09       ','BH        ','VV        '/)
68
69  ! Tester variables and arrays
70  integer :: cellMesh(3) = (/nx,nx,nx/)
71  integer :: i1, i1max, i2, i2max, i3, i3max, ir, irmax, &
72             lb1, lb2, lb3, mSpin, myNode, nNodes, one, two, ub1, ub2, ub3
73  real(dp):: atomDens(nr,nSpin), atomEc, atomEx, atomDc, atomDx, &
74             atomVxc(nr,nSpin), avgDiffVxc, &
75             cell(3,3), cellEc, cellEx, cellDc, cellDx, &
76             d0, d0s(nSpin), diffVxc(nSpin), dr, dx, Ecut, kCut, &
77             latConst, maxDiffVxc, pi, r, recCell(3,3), rMesh(nr), &
78             stress(3,3), sumDiffVxc, tmp, Vxc(nSpin), &
79             wc(nfTot), wr, wx(nfTot), x(3), x0(3)
80  real(gp),allocatable:: cellDens(:,:,:,:), cellVxc(:,:,:,:)
81
82#ifdef MPI
83  ! MPI-related variables
84  integer :: MPIerror
85  integer :: nLarger, nxNode
86#endif
87
88#ifdef MPI
89  ! Initialize MPI and get myNode and nNodes
90  call MPI_Init( MPIerror )
91  call MPI_Comm_Rank( MPI_Comm_World, myNode, MPIerror )
92  call MPI_Comm_Size( MPI_Comm_World, nNodes, MPIerror )
93#else
94  myNode = 0
95  nNodes = 1
96#endif
97
98  ! Initialize hybrid XC functional with all tested functionals
99  mSpin = min( nSpin, 2 )
100  wx = 1._dp / nf
101  wc = 1._dp / nf
102  call setXC( nf, func(indexf), auth(indexf), wx(indexf), wc(indexf) )
103
104  ! Find radial mesh points and gaussian density
105  pi = acos(-1._dp)
106  d0 = Qtot / (2*pi*dWidth**2)**1.5_dp    ! Total density at origin
107  if (nSpin==1) then
108    d0s(1) = d0
109  else
110    one = 1   ! A silly thing to satisfy the compiler when nSpin=1
111    two = 2
112    d0s(one) = d0 * (Qtot + spinPol) / Qtot / 2 ! Spin up density at origin
113    d0s(two) = d0 * (Qtot - spinPol) / Qtot / 2 ! Spin down density at origin
114  end if
115  dr = rmax / (nr-1)                      ! Interval between radial points
116  do ir = 1,nr
117    rMesh(ir) = dr * (ir-1)               ! Radial point values
118    atomDens(ir,:) = DensOfR( d0s(:), rMesh(ir) )
119  end do
120
121  ! Find exchange and correlation energy and potential from radial density
122  call atomXC( irel, nr, nr, rMesh, nSpin, atomDens, &
123               atomEx, atomEc, atomDx, atomDc, atomVxc )
124
125  ! Define fcc unit cell, such that a sphere of radius rMax+rBuff fits in it
126  latConst = (rMax+rBuff) * 2*sqrt(2._dp)
127  cell(:,1) = (/ 0.0_dp, 0.5_dp, 0.5_dp /)
128  cell(:,2) = (/ 0.5_dp, 0.0_dp, 0.5_dp /)
129  cell(:,3) = (/ 0.5_dp, 0.5_dp, 0.0_dp /)
130  cell(:,:) = cell(:,:) * latConst
131
132  ! Define reciprocal unit cell
133  recCell(:,1) = (/-1.0_dp, 1.0_dp, 1.0_dp /)
134  recCell(:,2) = (/ 1.0_dp,-1.0_dp, 1.0_dp /)
135  recCell(:,3) = (/ 1.0_dp, 1.0_dp,-1.0_dp /)
136  recCell(:,:) = recCell(:,:) * 2*pi/latConst
137  kCut = cellMesh(1) * sqrt(sum(recCell(:,1)**2)) / 2  ! Max. wave vector
138  Ecut = kCut**2                                       ! Mesh cutoff, in Ry
139  dx = pi / kCut                                 ! Dist. between mesh planes
140
141  ! Find the box of mesh points own by my processor
142#ifdef MPI
143  ! Do simplest thing: divide only along first axis
144  nxNode = Nx / nNodes          ! Points per node along first vector
145  nLarger = nx - nxNode*nNodes  ! Number of nodes with one more point
146  if (myNode<nLarger) then      ! My node has nx+1 points
147    lb1 = (nxNode+1)*(myNode-1)
148    ub1 = (nxNode+1)*(myNode-1) - 1
149  else                          ! My node has nx points
150    lb1 = (nxNode+1)*nLarger + nxNode*(myNode-nLarger)
151    ub1 = (nxNode+1)*nLarger + nxNode*(myNode-nLarger+1) - 1
152  end if
153#else
154  ! All points belong to the only processor
155  lb1 = 0
156  ub1 = nx-1
157#endif
158  lb2 = 0
159  lb3 = 0
160  ub2 = nx-1
161  ub3 = nx-1
162
163  ! Allocate arrays for density and potential
164  allocate( cellDens(lb1:ub1,lb2:ub2,lb3:ub3,nSpin), &
165             cellVxc(lb1:ub1,lb2:ub2,lb3:ub3,nSpin) )
166
167  ! Find density at mesh points
168  x0(:) = sum(cell,2) / 2     ! Center of cell
169  do i3 = lb3,ub3
170  do i2 = lb2,ub2
171  do i1 = lb1,ub1
172    x(:) = i1*cell(:,1)/cellMesh(1) &   ! Mesh point position
173         + i2*cell(:,2)/cellMesh(2) &
174         + i3*cell(:,3)/cellMesh(3)
175    r = sqrt( sum((x-x0)**2) )          ! Distance to center of cell
176    cellDens(i1,i2,i3,:) = DensOfR( d0s(:), r )
177  end do ! i1
178  end do ! i2
179  end do ! i3
180
181  ! Find exchange and correlation energy and potential from density in cell
182  call cellXC( irel, cell, cellMesh, lb1, ub1, lb2, ub2, lb3, ub3, nSpin, &
183               cellDens, cellEx, cellEc, cellDx, cellDc, stress, cellVxc )
184
185  ! Print parameters
186  if (myNode==0) then
187    print'(/,a,10(a3,4x))', 'funcs= ', func(indexf)
188    print  '(a,10(a6,1x))', 'auths= ', auth(indexf)
189    print'(a,3f12.6)', 'dr, dx, Ecut = ', dr, dx, Ecut
190    print'(a,2f12.6)', 'rMax, rBuff = ', rMax, rBuff
191  end if
192
193  ! Compare energies
194  if (myNode==0) then
195    print'(a,3f15.9)', 'atomEx, cellEx, diff =', atomEx, cellEx, atomEx-cellEx
196    print'(a,3f15.9)', 'atomEc, cellEc, diff =', atomEc, cellEc, atomEc-cellEc
197    print'(a,3f15.9)', 'atomDx, cellDx, diff =', atomDx, cellDx, atomDx-cellDx
198    print'(a,3f15.9)', 'atomDc, cellDc, diff =', atomDc, cellDc, atomDc-cellDc
199  end if
200
201  ! Write potentials
202  if (nNodes==1) then
203    open( unit=44, file='atomVxc.out' )
204    do ir = 1,nr
205      if (sum(atomDens(ir,1:mSpin)) < densMin) cycle
206      write(44,'(3f15.9)') rMesh(ir), atomVxc(ir,1), atomVxc(ir,mSpin)
207    end do
208    close( unit=44 )
209    open( unit=44, file='cellVxc.out' )
210    do i1 = cellMesh(1)/2,0,-1
211      i2 = cellMesh(2) / 2
212      i3 = cellMesh(3) / 2
213      if (sum(cellDens(i1,i2,i3,1:mSpin)) < densMin) cycle
214      x(:) = i1*cell(:,1)/cellMesh(1) &   ! Mesh point position
215           + i2*cell(:,2)/cellMesh(2) &
216           + i3*cell(:,3)/cellMesh(3)
217      write(44,'(3f15.9)') &
218        sqrt(sum((x-x0)**2)), cellVxc(i1,i2,i3,1), cellVxc(i1,i2,i3,mSpin)
219    end do
220    close( unit=44 )
221  end if
222
223  ! Compare potentials
224  sumDiffVxc = 0
225  maxDiffVxc = 0
226  do i3 = lb3,ub3
227  do i2 = lb2,ub2
228  do i1 = lb1,ub1
229    if (sum(cellDens(i1,i2,i3,1:mSpin)) < densMin) cycle
230    x(:) = cell(:,1)*i1/cellMesh(1) &
231         + cell(:,2)*i2/cellMesh(2) &
232         + cell(:,3)*i3/cellMesh(3)
233    r = sqrt(sum((x-x0)**2))
234    if (r>=rMax) cycle
235    ! Simplest thing: a linear interpolation of atomVxc (requires large nr)
236    ir = r/dr + 1
237    wr = (r - ir*dr) / dr
238    Vxc(:) = atomVxc(ir,:)*(1-wr) + atomVxc(ir+1,:)*wr
239    diffVxc(:) = abs(cellVxc(i1,i2,i3,:)-Vxc(:))
240    if (maxval(diffVxc(:)) > maxDiffVxc) then
241      i1max = i1
242      i2max = i2
243      i3max = i3
244      irmax = ir
245    end if
246    sumDiffVxc = sumDiffVxc + sum(diffVxc(:)**2)
247    maxDiffVxc = max( maxDiffVxc, maxval(diffVxc(:)) )
248  end do ! i1
249  end do ! i2
250  end do ! i3
251#ifdef MPI
252  ! Find sumDiffVxc and maxDiffVxc accross all processor nodes
253  tmp = sumDiffVxc
254  call MPI_AllReduce( tmp, sumDiffVxc, 1, MPI_double_precision, &
255                      MPI_Sum, MPI_Comm_World, MPIerror )
256  tmp = maxDiffVxc
257  call MPI_AllReduce( tmp, maxDiffVxc, 1, MPI_double_precision, &
258                      MPI_Max, MPI_Comm_World, MPIerror )
259#endif
260  avgDiffVxc = sumDiffVxc / nSpin / nx**3
261  if (myNode==0) then
262    print'(a,2f15.9)', 'avgDiffVxc, maxDiffVxc = ', avgDiffVxc, maxDiffVxc
263!    print'(a,4i6)', 'i1max,i2max,i3max, irmax = ', i1max, i2max, i3max, irmax
264  end if
265
266! Finalize MPI
267#ifdef MPI
268  call MPI_Finalize( MPIerror )
269#endif
270
271CONTAINS
272
273FUNCTION DensOfR( d0, r )
274
275  ! Returns a radial density distribution
276
277  implicit none
278  real(dp),intent(in):: d0(nSpin)  ! Density at center of charge distribution
279  real(dp),intent(in):: r          ! Distance to center of charge distribution
280  real(dp)           :: DensOfR(nSpin)  ! Electron density
281
282  ! Use a simple gaussian distribution
283  DensOfR = d0 * exp(-r**2/2/dWidth**2)
284
285  ! Impose a smooth radial cutoff
286  DensOfR = DensOfR * ( 1 - (r/rMax)**n1cut )**n2cut
287
288END FUNCTION DensOfR
289
290END PROGRAM siestaXCtest4
291
292