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