1 subroutine dimqm_screening(rtdb, geom, basis, dimxyz, screen) 2c 3c Subroutine to calculate the electric field due to the QM electrons 4c 5c Called from dimqm_main.F 6 7c use constants 8 implicit none 9#include "nwc_const.fh" 10#include "errquit.fh" 11#include "basP.fh" 12#include "basdeclsP.fh" 13#include "geomP.fh" 14#include "geobasmapP.fh" 15#include "mafdecls.fh" 16#include "bas_exndcf_dec.fh" 17#include "bas_ibs_dec.fh" 18#include "int_nbf.fh" 19#include "stdio.fh" 20#include "apiP.fh" 21#include "util.fh" 22#include "bas.fh" 23#include "dimqm.fh" 24 25c Input Variables 26 integer rtdb 27 integer geom 28 integer basis 29 double precision dimxyz(3, nDIM) 30 double precision screen(nDIM) 31c 32 integer nshell, nbf_max 33 integer ishell, ilo, ihi, jatom 34 integer ibas, ucont, itype, inp, igen, iexp, icf, iatom, igeom 35 integer lens 36 double precision S(isz_1e) 37 double precision scr(isz_1e) 38 double precision overlap(nDIM) 39 integer i_nbf_x, i_nbf_s, j_nbf_x, j_nbf_s 40#include "bas_exndcf_sfn.fh" 41#include "bas_ibs_sfn.fh" 42 43 if (.not. bas_geom(basis, geom)) call errquit 44 $ ('screening: bad basis', 555, BASIS_ERR) 45 if (.not. bas_numcont(basis, nshell)) call errquit 46 $ ('screening: bas_numcont failed for basis', basis, BASIS_ERR) 47 if (.not. bas_nbf_cn_max(basis,nbf_max)) call errquit 48 & ('screening: bas_nbf_cn_max failed',555, BASIS_ERR) 49 lens = sqrt(real(isz_1e)) 50 S = 0.0d0 51c write(luout,*) "len?", isz_1e 52 overlap = 0.0d0 53 do ishell = 1, nshell 54 if (.not. bas_cn2bfr(basis, ishell, ilo, ihi)) call errquit 55 & ('screening: bas_cn2bfr failed for i',basis,BASIS_ERR) 56 call int_nogencont_check(basis,'screening:basis') 57 ibas = basis + basis_handle_offset 58 ucont = (sf_ibs_cn2ucn(ishell,ibas)) 59 itype = infbs_cont(CONT_TYPE ,ucont,ibas) 60 inp = infbs_cont(CONT_NPRIM,ucont,ibas) 61 igen = infbs_cont(CONT_NGEN ,ucont,ibas) 62 iexp = infbs_cont(CONT_IEXP ,ucont,ibas) 63 icf = infbs_cont(CONT_ICFP ,ucont,ibas) 64 iatom = (sf_ibs_cn2ce(ishell,ibas)) 65 igeom = ibs_geom(ibas) 66 do jatom = 1, nDIM 67c write(luout,*) "DIM: ", jatom 68c write(luout,*) "Shell: ", ishell 69 call dimqm_screen2(coords(1,iatom,igeom), 70 $ dbl_mb(mb_exndcf(iexp,ibas)), 71 $ dbl_mb(mb_exndcf(icf, ibas)), 72 $ inp, igen, itype, dimxyz(:, jatom), 73 $ S, lens) 74c i_nbf_x = int_nbf_x(itype) 75c i_nbf_s = int_nbf_s(itype) 76c j_nbf_x = int_nbf_x(0) 77c j_nbf_s = j_nbf_x 78c call spcart_tran1e(S, scr, 79c $ j_nbf_x, i_nbf_x, 0, 1, 80c $ j_nbf_s, i_nbf_s, itype, igen, .false.) 81c write(luout,*) "Transform" 82c write(luout,*) S(1:lens) 83 overlap(jatom) = overlap(jatom) + sum(abs(S)) 84 end do 85 86 end do 87c 88 do jatom = 1, nDIM 89 screen(jatom) = erfc(overlap(jatom)) 90 end do 91c write(luout,*) "Total Overlap:" 92c write(luout,*) overlap 93c write(luout,*) "Screen factors:" 94c write(luout,*) erfc(overlap) 95 end subroutine dimqm_screening 96 97 subroutine dimqm_screen2(xyzi, expi, coefi, i_nprim, i_ngen, Li, 98 $ dimxyz, s, lens) 99 implicit none 100#include "nwc_const.fh" 101#include "hnd_rys.fh" 102#include "stdio.fh" 103#include "hnd_tol.fh" 104#include "dimqm.fh" 105#include "dimqm_constants.fh" 106 107 double precision xyzi(3) 108 double precision expi(i_nprim) 109 double precision coefi(i_nprim, i_ngen) 110 integer i_nprim 111 integer i_ngen 112 integer Li 113 double precision dimxyz(3) 114 double precision s(lens) 115 integer lens 116c 117 double precision xi, yi, zi 118 integer lit, maxi 119 double precision xj, yj, zj 120 integer ljt, maxj, ljtmod 121 double precision rr 122 double precision tol 123c 124 integer ig 125 double precision ai, arri, axi, ayi, azi, csi, csj 126 double precision aj, aa, aa1, dum, ax, ay, az, fac 127c 128 double precision dij, t 129 integer i, j, ij 130 double precision xint, yint, zint 131 double precision xs(Li+1, 3), ys(Li+1, 3), zs(Li+1, 3) ! Lj = 0 since all DIM atoms are spherical gaussians 132 integer Nxyz(3), ix, iy, iz, jx, jy, jz 133 134 tol = 2.30258d+00*itol 135c 136c ---- i shell ---- 137c 138 xi = xyzi(1) 139 yi = xyzi(2) 140 zi = xyzi(3) 141 lit = Li + 1 142 maxi = lit*(lit+1)/2 143c write(luout,*) "ngen:", i_ngen 144c 145c ---- DIM atom ---- 146c 147 xj = dimxyz(1) 148 yj = dimxyz(2) 149 zj = dimxyz(3) 150 ljt = 1 ! All DIM atoms are spherical (Lj = 0, so Lj+1 = 1) 151 ljtmod = ljt + 2 152 maxj = ljt*(ljt+1)/2 ! Always 1 153 rr = (xi - xj)**2 + (yi - yj)**2 + (zi - zj)**2 154 nroots = (lit+ljt - 2)/2 + 1 155 s = zero 156c 157c --- i primitive --- 158c 159 do 7000 ig = 1, i_nprim 160 ai = expi(ig) 161 arri = ai*rr 162 axi = ai*xi 163 ayi = ai*yi 164 azi = ai*zi 165 csi = coefi(ig, i_ngen) 166c write(luout,*) csi 167c 168c --- DIM primitive --- 169c 170 ! No loop cause DIM atoms only have 1 gaussian 171 aj = 1.0d0 ! Width of DIM gaussian 172 aa = ai+aj 173 aa1 = one/aa 174 dum = aj*arri*aa1 175 if(dum .gt. tol) then 176c write(luout,*) "value", dum 177c write(luout,*) "greater than tolerance", tol 178c write(luout,*) "Overlap is zero" 179 go to 7000 180 end if 181 fac = exp(-dum) 182 csj = one ! Contraction coefficent to DIM gaussian 183 ax = (axi + aj*xj)*aa1 184 ay = (ayi + aj*yj)*aa1 185 az = (azi + aj*zj)*aa1 186c 187c ---- density factor ---- 188c 189 dij = fac * csi * csj 190c write(luout,*) "dij:", dij 191c 192c --- overlap --- 193c 194 t = sqrt(aa1) 195 do j = 1, ljtmod 196 do i = 1, lit 197 call dimqm_overlap(ax, ay, az, xi, yi, zi, 198 $ xj, yj, zj, i, j, t, 199 $ xint, yint, zint) 200 xs(i,j) = xint*t 201 ys(i,j) = yint*t 202 zs(i,j) = zint*t 203c write(luout,*) xint, yint, zint 204 end do ! i = 1, lit 205 end do ! j = 1, ljtmod 206c 207 ij = 0 208 do i = 1, maxi 209 call getNxyz(Li, i, Nxyz) 210 ix = Nxyz(1) + 1 211 iy = Nxyz(2) + 1 212 iz = Nxyz(3) + 1 213 do j = 1, maxj 214 call getNxyz(0, j, Nxyz) 215 jx = Nxyz(1) + 1 ! Should always be 1 216 jy = Nxyz(2) + 1 ! Should always be 1 217 jz = Nxyz(3) + 1 ! Should always be 1 218 ij = ij + 1 219 dum = xs(ix,jx)*ys(iy,jy)*zs(iz,jz) 220 s(ij) = s(ij) + dum*dij 221 end do ! j = 1, maxj 222 end do ! i = 1, maxi 2237000 end do ! ig = 1, i_nprim 224c write(luout,*) s 225 226 end subroutine dimqm_screen2 227 228 subroutine dimqm_overlap(x0, y0, z0, xi, yi, zi, xj, yj, zj, 229 $ ni, nj, t, xint, yint, zint) 230 implicit none 231#include "hnd_whermt.fh" 232#include "dimqm_constants.fh" 233#include "stdio.fh" 234 double precision x0, y0, z0 235 double precision xi, yi, zi 236 double precision xj, yj, zj 237 integer ni, nj 238 double precision t, xint, yint, zint 239c 240 integer i, imin, imax, ii, jj, npts 241 double precision dum 242 double precision px, py, pz 243 double precision ax, ay, az 244 double precision bx, by, bz 245 double precision ptx, pty, ptz 246 247 xint = zero 248 yint = zero 249 zint = zero 250 npts = (ni + nj - 2)/2 + 1 251c write(luout,*) "npts:", npts 252 imin = hermin(npts) 253 imax = hermax(npts) 254c write(luout,*) "min, max", imin, imax 255 do i = imin, imax 256 dum = w(i) 257c write(luout,*) "w", dum 258 px = dum 259 py = dum 260 pz = dum 261 dum = h(i) * t 262 ptx = dum+x0 263 pty = dum+y0 264 ptz = dum+z0 265 ax = ptx-xi 266 ay = pty-yi 267 az = ptz-zi 268 bx = ptx-xj 269 by = pty-yj 270 bz = ptz-zj 271 do ii = 1, ni-1 272 px = px*ax 273 py = py*ay 274 pz = pz*az 275 end do 276 do jj = 1, nj-1 277 px = px*bx 278 py = py*by 279 pz = pz*bz 280 end do 281 xint = xint + px 282 yint = yint + py 283 zint = zint + pz 284 end do 285 return 286 287 end subroutine dimqm_overlap 288