1PROGRAM siestaXCtest5 2 3 ! Finds and writes the integrand of the nonlocal vdW functional 4 ! for fixed kf2 and kg2=grad(n2)/n2 in a mesh of kf1, kg1, and r12. 5 ! J.M.Soler. Jul.2012 6 7 ! Used module procedures 8 USE siestaXC, only: setXC 9 USE m_vdwxc, only: phiofr 10 USE m_vdwxc, only: vdw_get_qmesh 11 USE m_vdwxc, only: vdw_set_kcut 12 USE m_vdwxc, only: vdw_theta 13 USE m_vv_vdwxc, only: vv_vdw_phi_val 14 15 ! Used module parameters 16 USE siestaXC, only: dp 17 USE siestaXC, only: gp => grid_p 18 19 implicit none 20 21 ! Tester parameters 22 integer, parameter:: nspin = 1 ! Number of spin components 23 integer, parameter:: nfunc = 1 ! Number of functionals 24 integer, parameter:: nr = 6 ! Number of radial points 25 real(dp),parameter:: r(nr) = (/0.1_dp, 1.0_dp, 2.0_dp, 5.0_dp, & 26 10._dp, 20._dp/) ! |r1-r2| distances 27 real(dp),parameter:: kcut = 15._dp ! Planewave vector cutoff 28 real(dp),parameter:: kf2 = 1.0_dp ! Fermi wavevector at point 2 29 real(dp),parameter:: kg2 = 2.5_dp ! grad(n)/n at point 2 30 real(dp),parameter:: kfmax = 5._dp ! Max. Fermi wavevector at point 1 31 real(dp),parameter:: kgmax = 5._dp ! Max. grad(n)/n at point 1 32 real(dp),parameter:: dkf = 0.1_dp ! Fermi wavevector mesh interval 33 real(dp),parameter:: dkg = 0.1_dp ! grad(n)/n mesh interval 34 character(len=*),parameter:: func = 'VDW' 35 character(len=*),parameter:: auth = 'VV' ! 'DRSLL'|'LMKLL'|'VV' 36 ! KBM, C09, and BH use DRSLL kernel 37 38 ! Tester variables and arrays 39 integer :: ik, ikf, ikg, ir, nk, nkf, nkg 40 real(dp):: kf1, kg1, gn1(3,nspin), gn2(3,nspin), n1(nspin), n2(nspin), & 41 n1n2phi_exact, n1n2phi_interp, pi, r12, wc(nfunc), wx(nfunc) 42 real(dp),allocatable:: & 43 dtdn1(:,:), dtdn2(:,:), dtdgn1(:,:,:), dtdgn2(:,:,:), & 44 n1phi(:), phi(:,:), theta1(:), theta2(:) 45 46 ! Initialize XC functional 47 wx = 1._dp/nfunc 48 wc = 1._dp/nfunc 49 call setXC( nfunc, func, auth, wx, wc ) 50 51 ! Set planewave cutoff 52 call vdw_set_kcut( kcut ) 53 54 ! Find number of interpolation points and allocate interpolation arrays 55 call vdw_get_qmesh( nk ) 56 allocate( dtdn1(nk,nspin), dtdn2(nk,nspin) ) 57 allocate( dtdgn1(3,nk,nspin), dtdgn2(3,nk,nspin) ) 58 allocate( n1phi(nk), phi(nk,nk), theta1(nk), theta2(nk) ) 59 60 ! Open file for output 61 open(1,file='n1n2phi.table') 62 63 ! Iterate on kf1, kg1, and r12 64 nkf = nint(kfmax/dkf) 65 nkg = nint(kgmax/dkg) 66 do ir = 1,nr 67 r12 = r(ir) 68 call phiofr( r12, phi ) 69! print'(a,/,(6f12.6))','phi(:,13)=',phi(:,13) 70! print'(a,/,(6f12.6))','phi(13,:)=',phi(13,:) 71 do ikf = 0,nkf 72 do ikg = 0,nkg 73 kf1 = ikf*dkf 74 kg1 = ikg*dkg 75 kf1 = max(kf1,1.e-6) 76 pi = acos(-1._dp) 77 n1 = 0 78 n2 = 0 79 gn1 = 0 80 gn2 = 0 81 n1(1) = kf1**3/(3*pi**2) 82 n2(1) = kf2**3/(3*pi**2) 83 gn1(1,1) = n1(1)*kg1 84 gn2(1,1) = n2(1)*kg2 85 call vdw_theta( nspin, n1, gn1, theta1, dtdn1, dtdgn1 ) 86 call vdw_theta( nspin, n2, gn2, theta2, dtdn2, dtdgn2 ) 87! if (ir==1 .and. ikf==0) then 88! print'(a,6f12.6)','r12,kf1,kg1=',r12,kf1,kg1 89! print'(a,/,(6f12.6))','theta1=',theta1 90! endif 91 do ik = 1,nk 92 n1phi(ik) = sum(theta1(:)*phi(:,ik)) 93 end do 94 95 ! Assume phi_val returns phi, with which we want to compare 96! n1n2phi_interp = sum(n1phi*theta2)/n1(1)/n2(1) 97 98 ! Assume phi_val returns sqrt(n1*n2)*phi 99! n1n2phi_interp = sum(n1phi*theta2) /sqrt(n1(1)*n2(1)) 100 101 ! Assume phi_val returns kf1*kf2*phi 102! n1n2phi_interp = sum(n1phi*theta2) *kf1*kf2/n1(1)/n2(1) 103 104 ! Assume phi_val returns (kf1*kf2)**2*phi 105 n1n2phi_interp = sum(n1phi*theta2) *(kf1*kf2)**2/n1(1)/n2(1) 106 107 ! Assume phi_val returns n1*n2*phi 108! n1n2phi_interp = sum(n1phi*theta2) 109 110 n1n2phi_exact = vv_vdw_phi_val( kf1, kf2, kg1, kg2, r12 ) 111 write(1,'(3f9.3,3e15.6)') r12, kf1, kg1, & 112 n1n2phi_exact, n1n2phi_interp, n1n2phi_interp-n1n2phi_exact 113 end do 114 end do 115 end do 116 117 close(1) 118 119END PROGRAM siestaXCtest5 120 121