1czora... 2czora... Calculate the zora correction on the grid 3czora... 4 subroutine grid_zoracorr(nqpts,qxyz,qwght,natoms,g_dens,amat) 5c 6 implicit none 7c 8#include "errquit.fh" 9#include "mafdecls.fh" 10#include "cdft.fh" 11#include "stdio.fh" 12#include "zora.fh" 13#include "geom.fh" 14#include "global.fh" ! FA 15c 16 external gridQpqPotential 17 integer nqpts 18 integer g_dens(2),igrid,natoms,npol 19 double precision qxyz(3,nqpts),qwght(nqpts) ![in] 20 double precision amat(nqpts,ipol) 21 double precision amat_coul(nqpts,ipol) 22 double precision amat_nucl(nqpts) 23 double precision amat_Qnucl(nqpts) ! Added by FA 24 double precision tol 25 double precision zoraCorr, totPot 26 integer iatom 27 28 double precision nucCharge, nucCoords(3) 29 character*16 tags(natoms) 30 logical lSuccess 31 double precision rx,ry,rz,dist 32 double precision denomFac 33 integer closegridpts(nqpts) 34c 35 double precision clight_au2 36 clight_au2 = clight_au*clight_au 37c 38c == preliminaries == 39 do igrid = 1,nqpts 40 amat(igrid,1) = 0.d0 41 amat_coul(igrid,1) = 0.d0 42 if (ipol.gt.1) then 43 amat(igrid,2) = 0.d0 44 amat_coul(igrid,2) = 0.d0 45 end if 46 amat_nucl(igrid) = 0.d0 47 amat_Qnucl(igrid) = 0.d0 ! Added by FA 48 closegridpts(igrid) = 0 49 end do 50c 51c == calculate the hartree potential on a supplied list of points == 52 tol = 1d-8 53 call potential_list(ao_bas_han, g_dens(1), nqpts, 54 & qxyz, amat_coul(1,1), tol) 55 56 if (ipol.gt.1) then 57 call potential_list(ao_bas_han, g_dens(2), nqpts, 58 & qxyz, amat_coul(1,2), tol) 59 end if 60c 61c == calculate the total nuclear potential on the grid == 62 call gridNuclearPotential(geom,natoms,nqpts,qxyz,qwght, 63 & closegridpts,amat_nucl) 64 if ((zora_calc_type.eq.3).or.(zora_calc_type.eq.4)) then 65c == calculate Quadrupole potential on the grid == 66 call gridQpqPotential(nqpts,qxyz,amat_Qnucl, 67 & closegridpts) 68 end if 69c 70c == assemble zora correction == 71 zoraCorr = 0.d0 72 totPot = 0.d0 73 do igrid = 1,nqpts 74 if (ipol.gt.1) then 75 totPot = -amat_coul(igrid,1)-amat_coul(igrid,2) 76 & + amat_nucl(igrid) 77 else 78 totPot = -amat_coul(igrid,1)+amat_nucl(igrid) 79 end if 80c 81c == assemble the appropriate correction == 82 if (zora_calc_type.eq.0) then ! pure kinetic test 83 zoraCorr = 0.5d0 84 else if (zora_calc_type.eq.1) then ! zora correction 85 zoraCorr = totPot/(4.d0*clight_au2 - 2.d0*totPot) 86 else if (zora_calc_type.eq.2) then ! zora energy scaling 87 denomFac = (2.d0*clight_au2 - totPot) 88 zoraCorr = clight_au2/denomFac/denomFac 89 else if (zora_calc_type.eq.3) then ! zora EFG - Added by FA 90 denomFac = (2.d0*clight_au2 - totPot) 91 zoraCorr = clight_au2/denomFac/denomFac 92 zoraCorr = zoraCorr*amat_Qnucl(igrid) 93 else if (zora_calc_type.eq.4) then ! Numerical EFG - Added by FA 94 zoraCorr = amat_Qnucl(igrid) 95c zoraCorr = 1.0 ! For Atomic Nr FA-04-05-10 96 else 97 zoraCorr = 0.d0 98 end if 99c 100c == multiply by the integration weight == 101 if (igrid.eq.closegridpts(igrid)) then 102 zoraCorr = 0.d0 103 else 104 zoraCorr = zoraCorr*qwght(igrid) 105 end if 106c 107c == fill array == 108 amat(igrid,1) = zoraCorr 109 if (ipol.gt.1) amat(igrid,2) = zoraCorr 110 end do 111c 112 return 113 end 114c $Id$ 115