1c
2c     == calculate EFGZ4-scalar relativistic (SR) contribution
3      subroutine calc_zora_EFGZ4_SR(ao_bas_han,    ! in: AO basis handle
4     &                              geom,          ! in: geometry handle
5     &                              ipol,          ! in: nr. of polarizations
6     &                              g_dens,        ! in: superposit. atomic densities
7     &                              delchi_ao,     ! in: deriv. of basis functions
8     &                              qxyz,          ! in: grid points
9     &                              qwght,         ! in: weighting coeffs.
10     &                              nbf,           ! in: nr. basis functions
11     &                              npts,          ! in: nr. grid points
12     &                              natoms,        ! in: nr. atoms
13     &                              zora_Qpq,      ! in : type of EFG potential
14     &                              xyz_EFGcoords, ! in : EFG-nuclear coordinates
15     &                              use_modelpotential,
16     &                              gexpo,
17     &                              gcoef,
18     &                              zora_efgz4)    ! out: munu-matrix
19      implicit none
20#include "errquit.fh"
21#include "mafdecls.fh"
22#include "stdio.fh"
23#include "zora.fh"
24#include "global.fh"
25#include "bas.fh"
26
27      integer nbf, npts, ao_bas_han, natoms, geom
28      integer g_dens(2),ipol
29      double precision qwght(npts),pot(npts)
30      double precision qxyz(3,npts)
31      double precision delchi_ao(npts,3,nbf)
32      integer i,j,k
33      double precision amat_coul(npts,ipol)
34      double precision amat_nucl(npts)
35      double precision amat_Qnucl(npts)
36      integer ipt,closegridpts(npts)
37      double precision clight_au2,tol
38      double precision amat_tot
39      double precision denom
40      double precision fac_arr(npts)
41      integer zora_Qpq
42      double precision xyz_EFGcoords(3)
43      double precision ac_efgz4
44      double precision zora_efgz4(nbf,nbf)
45      logical use_modelpotential
46      double precision gexpo(natoms,50)
47      double precision gcoef(natoms,50)
48c
49      external get_ints_zora_efgz4_sr,gridQpqPotential
50
51      clight_au2 = clight_au*clight_au
52c
53c     == preliminaries ==
54      do ipt = 1,npts
55        do i=1,ipol
56         amat_coul(ipt,i) = 0.d0
57        end do
58        amat_nucl(ipt) = 0.d0
59        closegridpts(ipt) = 0
60      end do
61c
62c     == calculate the total hartree potential on the grid ==
63      call gridHartreePotential(use_modelpotential,
64     &    ao_bas_han, geom, natoms, ipol, g_dens, npts, qxyz, qwght,
65     &    closegridpts, gexpo, gcoef, amat_coul)
66c
67c     == calculate the total nuclear potential on the grid ==
68      call gridNuclearPotentialPoint(geom,natoms,npts,qxyz,qwght,
69     &                          closegridpts,amat_nucl)
70      call gridQpqPotential(zora_Qpq,xyz_EFGcoords,
71     &                      npts,qxyz,
72     &                      amat_Qnucl, ! out: EFG potential
73     &                      closegridpts)
74      do k = 1,npts
75        if (k.eq.closegridpts(k)) qwght(k) = 0.d0
76      end do
77c     === define fac1_arr,fac2_arr -- FA
78      do k = 1,npts
79c      == assemble hartree and nuclear contributions ==
80       amat_tot = amat_coul(k,1)+amat_nucl(k)
81       denom = (2.d0*clight_au2 - amat_tot)
82       fac_arr(k)=(clight_au2/denom/denom)*amat_Qnucl(k)
83     &             *qwght(k)
84      end do
85c     == assemble zora correction ==
86c -----main diagonal --- START
87      do i = 1, nbf
88         j=i
89         call get_ints_zora_efgz4_sr(nbf,npts,delchi_ao,i,j,
90     &                               fac_arr,
91     &                               ac_efgz4)  ! out
92         zora_efgz4(i,j) = zora_efgz4(i,j) + ac_efgz4
93      enddo ! end-loop-i
94c -----main diagonal --- END
95c ----- off diagonal --- START
96      do i = 1, nbf
97        do j = i+1, nbf
98         call get_ints_zora_efgz4_sr(nbf,npts,delchi_ao,i,j,
99     &                               fac_arr,
100     &                               ac_efgz4)  ! out
101         zora_efgz4(i,j) = zora_efgz4(i,j)  + 2.0d0*ac_efgz4
102        enddo ! end-loop-j
103      enddo ! end-loop-i
104c ----- off diagonal --- END
105      return
106      end
107      subroutine get_ints_zora_efgz4_sr(nbf,            ! in: nr. basis functions
108     &                                  npts,           ! in: grid points
109     &                                  delchi_ao,      ! in: deriv. of basis functions
110     &                                  i,j,            ! in: (i,j) indices for delchi_ao
111     &                                  fac_arr,        ! in
112     &                                  ac_efgz4)       ! out
113      implicit none
114#include "errquit.fh"
115#include "stdio.fh"
116#include "global.fh"
117      integer nbf,npts,i,j,k
118      double precision delchi_ao(npts,3,nbf)
119      double precision fac_arr(npts)
120      double precision ac_efgz4
121      double precision prod0
122          ac_efgz4 = 0.0d0
123          do k = 1, npts
124           prod0 = delchi_ao(k,1,i)*delchi_ao(k,1,j)
125     &            +delchi_ao(k,2,i)*delchi_ao(k,2,j)
126     &            +delchi_ao(k,3,i)*delchi_ao(k,3,j)
127           ac_efgz4 = ac_efgz4 + fac_arr(k)*prod0
128          enddo ! end-loo-k
129      return
130      end
131c $Id$
132