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