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