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