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