1c =============== Fredy Aquino's routines ======= START 2c Another version of grid_zoracorr 3 subroutine grid_zoracorr_FA(nqpts,qxyz,qwght, 4 & 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" 14c 15 integer nqpts,iatom 16 integer g_dens(2),igrid,natoms,npol 17 integer closegridpts(nqpts) 18 double precision qxyz(*),qwght(*) ! [ in ] 19 double precision amat(nqpts,ipol) ! [ out ] 20 double precision amat_coul(nqpts,ipol) 21 double precision amat_nucl(nqpts) 22 double precision amat_Qnucl(nqpts) ! Added by FA 23 double precision tol,rx,ry,rz,dist 24 double precision nucCharge, nucCoords(3) 25 character*16 tags(natoms) 26 logical lSuccess 27 external fn0,fn1,fn2,fn3,fn4 28 external gridQpqPotential1 29c 30c == preliminaries == 31 do igrid = 1,nqpts 32 amat(igrid,1) = 0.d0 33 amat_coul(igrid,1) = 0.d0 34 if (ipol.gt.1) then 35 amat(igrid,2) = 0.d0 36 amat_coul(igrid,2) = 0.d0 37 end if 38 amat_nucl(igrid) = 0.d0 39 amat_Qnucl(igrid) = 0.d0 ! Added by FA 40 closegridpts(igrid) = 0 41 end do 42c 43c == calculate the hartree potential on a supplied list of points == 44 tol = 1d-8 45 call potential_list(ao_bas_han, g_dens(1), nqpts, 46 & qxyz, amat_coul(1,1), tol) 47 48 if (ipol.gt.1) then 49 call potential_list(ao_bas_han, g_dens(2), nqpts, 50 & qxyz, amat_coul(1,2), tol) 51 end if 52c 53c == calculate the total nuclear potential on the grid == 54 call gridNuclearPotential(geom,natoms,nqpts,qxyz,qwght, 55 & closegridpts,amat_nucl) 56 if ((zora_calc_type.eq.3).or.(zora_calc_type.eq.4)) then 57c == calculate Quadrupole potential on the grid == 58 call gridQpqPotential1(nqpts,qxyz,amat_Qnucl, 59 & amat_Qnucl) 60 end if 61c 62c == assemble zora correction == 63 if (zora_calc_type.eq.0) then ! pure kinetic test 64 call get_amat(fn0,nqpts,qwght, 65 & amat_coul,amat_nucl,amat_Qnucl, 66 & amat,closegridpts) 67 else if (zora_calc_type.eq.1) then ! zora correction 68 call get_amat(fn1,nqpts,qwght, 69 & amat_coul,amat_nucl,amat_Qnucl, 70 & amat,closegridpts) 71 else if (zora_calc_type.eq.2) then ! zora energy scaling 72 call get_amat(fn2,nqpts,qwght, 73 & amat_coul,amat_nucl,amat_Qnucl, 74 & amat,closegridpts) 75 else if (zora_calc_type.eq.3) then ! zora EFG 76 call get_amat(fn3,nqpts,qwght, 77 & amat_coul,amat_nucl,amat_Qnucl, 78 & amat,closegridpts) 79 else if (zora_calc_type.eq.4) then ! num EFG 80 call get_amat(fn4,nqpts,qwght, 81 & amat_coul,amat_nucl,amat_Qnucl, 82 & amat,closegridpts) 83 endif 84 return 85 end 86 87 subroutine get_amat(fcn,nqpts,qwght, 88 & amat_coul,amat_nucl,amat_Qnucl, 89 & amat,closegridpts) 90#include "errquit.fh" 91#include "mafdecls.fh" 92#include "cdft.fh" 93#include "stdio.fh" 94#include "zora.fh" 95#include "geom.fh" 96 integer nqpts,igrid 97 integer closegridpts(*) 98 double precision qwght(*) ![in] 99 double precision amat(nqpts,ipol) 100 double precision amat_coul(nqpts,ipol) 101 double precision amat_nucl(nqpts),amat_Qnucl(nqpts) 102 double precision valfn,totPot 103 external fcn 104 105 do igrid = 1,nqpts 106 if (ipol.gt.1) then 107 totPot = -amat_coul(igrid,1)-amat_coul(igrid,2) 108 & + amat_nucl(igrid) 109 else 110 totPot = -amat_coul(igrid,1)+amat_nucl(igrid) 111 end if 112 valfn=0.0 113 if (igrid.ne.closegridpts(igrid)) then 114 valfn=fcn(totPot,amat_Qnucl(igrid)) 115 endif 116 amat(igrid,1)=valfn*qwght(igrid) 117 if (ipol.gt.1) amat(igrid,2) = valfn*qwght(igrid) 118 end do 119 return 120 end 121 122 double precision function fn0(totPot,Qnucl) 123 double precision toPot,Qnucl 124 fn0=0.5d0 125 return 126 end 127 128 double precision function fn1(totPot,Qnucl) 129#include "errquit.fh" 130#include "mafdecls.fh" 131#include "cdft.fh" 132#include "stdio.fh" 133#include "zora.fh" 134#include "geom.fh" 135 double precision totPot,Qnucl 136 double precision clight_au2 137 clight_au2 = clight_au*clight_au 138 fn1=totPot/(4.d0*clight_au2-2.d0*totPot) 139 return 140 end 141 142 double precision function fn2(totPot,Qnucl) 143#include "errquit.fh" 144#include "mafdecls.fh" 145#include "cdft.fh" 146#include "stdio.fh" 147#include "zora.fh" 148#include "geom.fh" 149 double precision totPot,Qnucl 150 double precision clight_au2 151 double precision denomFac 152 clight_au2 = clight_au*clight_au 153 denomFac = (2.d0*clight_au2-totPot) 154 fn2=clight_au2/denomFac/denomFac 155 return 156 end 157 158 double precision function fn3(totPot,Qnucl) 159#include "errquit.fh" 160#include "mafdecls.fh" 161#include "cdft.fh" 162#include "stdio.fh" 163#include "zora.fh" 164#include "geom.fh" 165 double precision totPot,Qnucl 166 double precision clight_au2 167 double precision denomFac 168 clight_au2 = clight_au*clight_au 169 denomFac = (2.d0*clight_au2-totPot) 170 fn3=clight_au2/denomFac/denomFac*Qnucl 171 return 172 end 173 174 double precision function fn4(totPot,Qnucl) 175#include "errquit.fh" 176#include "mafdecls.fh" 177#include "cdft.fh" 178#include "stdio.fh" 179#include "zora.fh" 180#include "geom.fh" 181 double precision totPot,Qnucl 182 fn4=Qnucl 183 return 184 end 185czora... 186czora...Calculate the Quadrupole potential on the grid 187czora... 188 subroutine gridQpqPotential1(nqpts,qxyz, 189 & amat_Qnucl, 190 & closegridpts) 191c Purpose: Evaluates Q_pq(\vec{r},\vec{r}') 192c Input: 193c zora_Qpq = 1 -> Evaluates Qxx 194c = 2 -> Evaluates Qyy 195c = 3 -> Evaluates Qzz 196c = 4 -> Evaluates Qxy 197c = 5 -> Evaluates Qxz 198c = 6 -> Evaluates Qyz 199c ===> zora_Qpq, defined in zora.fh 200c nqpts , number of grid points 201c qxyz , grid points 202c xyz_EFGcoords, Quadrupole potential is evaluated 203c in this point (\vec{r}) 204c ===> xyz_EFGcoords, defined in zora.fh 205c Output: 206c amat_Qnucl, Quadrupole potential ev. in the grid 207c for integration purpose (\vec{r}') 208c 209 implicit none 210c 211#include "global.fh" 212#include "stdio.fh" 213#include "zora.fh" 214 215c 216 integer igrid,nqpts 217 double precision qxyz(3,nqpts) 218 double precision rx,ry,rz,dist 219 double precision amat_Qnucl(nqpts),dist5 220 integer closegridpts(*) 221 external fxx,fyy,fzz,fxy,fxz,fyz 222 223c == loop over the grid points == 224 if (zora_Qpq.eq.1) then ! Qxx 225 call ev_Qpot(fxx,nqpts,qxyz,closegridpts,amat_Qnucl) 226 else if (zora_Qpq.eq.2) then ! Qyy 227 call ev_Qpot(fyy,nqpts,qxyz,closegridpts,amat_Qnucl) 228 else if (zora_Qpq.eq.3) then ! Qzz 229 call ev_Qpot(fzz,nqpts,qxyz,closegridpts,amat_Qnucl) 230 else if (zora_Qpq.eq.4) then ! Qxy 231 call ev_Qpot(fxy,nqpts,qxyz,closegridpts,amat_Qnucl) 232 else if (zora_Qpq.eq.5) then ! Qxz 233 call ev_Qpot(fxz,nqpts,qxyz,closegridpts,amat_Qnucl) 234 else if (zora_Qpq.eq.6) then ! Qyz 235 call ev_Qpot(fyz,nqpts,qxyz,closegridpts,amat_Qnucl) 236 endif 237 return 238 end 239 240 subroutine ev_Qpot(fcn,nqpts,qxyz, 241 & closegridpts,amat_Qnucl) 242#include "stdio.fh" 243#include "zora.fh" 244c 245 integer igrid,nqpts 246 double precision qxyz(3,nqpts) 247 double precision rx,ry,rz,dist 248 double precision amat_Qnucl(nqpts),dist5 249 integer closegridpts(*) 250 external fcn 251 252 do igrid = 1,nqpts 253 amat_Qnucl(igrid) = 0.d0 254c == distance from the grid points to given xyz_EFGcoords() == 255 rx = xyz_EFGcoords(1) - qxyz(1,igrid) 256 ry = xyz_EFGcoords(2) - qxyz(2,igrid) 257 rz = xyz_EFGcoords(3) - qxyz(3,igrid) 258 dist = dsqrt(rx*rx + ry*ry + rz*rz) 259 if (dist.gt.zoracutoff_EFG) then 260 dist5=dist*dist*dist*dist*dist 261 amat_Qnucl(igrid)=fcn(rx,ry,rz,dist5) 262 else 263 closegridpts(igrid) = igrid 264 end if 265 end do ! end-grid 266 return 267 end 268 269 double precision function fxx(rx,ry,rz,dist5) 270 double precision rx,ry,rz,dist5 271 fxx=(2.d0*rx*rx-(ry*ry+rz*rz))/dist5 272 return 273 end 274 275 double precision function fyy(rx,ry,rz,dist5) 276 double precision rx,ry,rz,dist5 277 fyy=(2.d0*ry*ry-(rx*rx+rz*rz))/dist5 278 return 279 end 280 281 double precision function fzz(rx,ry,rz,dist5) 282 double precision rx,ry,rz,dist5 283 fzz=(2.d0*rz*rz-(rx*rx+ry*ry))/dist5 284 return 285 end 286 287 double precision function fxy(rx,ry,rz,dist5) 288 double precision rx,ry,rz,dist5 289 fxy=(3.d0*rx*ry)/dist5 290 return 291 end 292 293 double precision function fxz(rx,ry,rz,dist5) 294 double precision rx,ry,rz,dist5 295 fxz=(3.d0*rx*rz)/dist5 296 return 297 end 298 299 double precision function fyz(rx,ry,rz,dist5) 300 double precision rx,ry,rz,dist5 301 fyz=(3.d0*ry*rz)/dist5 302 return 303 end 304c =============== Fredy Aquino's routines ======= END 305c 306c $Id$ 307