mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-09 07:33:53 +01:00
626 lines
17 KiB
Fortran
626 lines
17 KiB
Fortran
!*****************************************************************************
|
|
subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result)
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
! Compute the overlap integral between a Gaussian function
|
|
! with arbitrary angular momemtum and a s-type Slater function
|
|
END_DOC
|
|
|
|
! Input variables
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
double precision,intent(out) :: result
|
|
|
|
! Final value of the integrals
|
|
double precision :: ss,ps,ds
|
|
double precision :: pxs,pys,pzs
|
|
double precision :: dxxs,dyys,dzzs,dxys,dxzs,dyzs
|
|
|
|
double precision :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k
|
|
|
|
pi = 4d0*atan(1d0)
|
|
|
|
! calculate the length AB between the two centers and other usful quantities
|
|
|
|
AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
|
|
AB = dsqrt(AB)
|
|
|
|
AxBx = (cGau(1)-cSla(1))/2d0
|
|
AyBy = (cGau(2)-cSla(2))/2d0
|
|
AzBz = (cGau(3)-cSla(3))/2d0
|
|
ds = 0.d0
|
|
|
|
! intermediate variables
|
|
|
|
t = expSla*dsqrt(0.25d0/expGau)
|
|
u = dsqrt(expGau)*AB
|
|
|
|
double precision :: d, et2
|
|
if(AB > 0d0) then
|
|
|
|
! (s|s)
|
|
ss = 0.d0
|
|
|
|
d = derfc(t+u)
|
|
if (dabs(d) > 1.d-30) then
|
|
ss = (t+u)*d*dexp(2d0*t*(t+u))
|
|
endif
|
|
|
|
d = derfc(t-u)
|
|
if (dabs(d) > 1.d-30) then
|
|
ss -= (t-u)*d*dexp(2d0*t*(t-u))
|
|
endif
|
|
|
|
! (p|s)
|
|
ps = 0.d0
|
|
if (t*t-u*u > 300.d0) then
|
|
et2 = huge(1.0)
|
|
else
|
|
et2 = dexp(t*t-u*u)
|
|
endif
|
|
if (et2 /= 0.d0) then
|
|
d = derfc(t-u)
|
|
if (d /= 0.d0) then
|
|
ps += dexp((t-u)**2)*(1d0+2d0*t*(t-u))*d
|
|
endif
|
|
d = derfc(t+u)
|
|
if (d /= 0.d0) then
|
|
ps += dexp((t+u)**2)*(1d0+2d0*t*(t+u))*d
|
|
endif
|
|
ps *= dsqrt(pi)
|
|
ps -= 4d0*t
|
|
ps *= et2/dsqrt(pi)
|
|
endif
|
|
|
|
! (d|s)
|
|
! ds = 4d0*dexp(2d0*t*(t-u))*t*(-((1d0+t**2-t*u)*derfc(t-u))+dexp(4d0*t*u)*(1d0+t*(t+u))*derfc(t+u))
|
|
ds = 0.d0
|
|
d = derfc(t+u)
|
|
if (d /= 0.d0) then
|
|
ds = dexp(4d0*t*u)*(1d0+t*(t+u))*d
|
|
endif
|
|
d = derfc(t-u)
|
|
if (d /= 0.d0) then
|
|
ds -= (1d0+t*t-t*u)*d
|
|
endif
|
|
|
|
if ( dabs(ds) > 1.d-100) then
|
|
ds *= 4d0*dexp(2d0*t*(t-u))*t
|
|
endif
|
|
|
|
! backward scaling
|
|
ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0
|
|
ps = ps/u**2-ss/u**3d0
|
|
ss = ss/u
|
|
|
|
else
|
|
|
|
! concentric case
|
|
d = derfc(t)
|
|
if (d /= 0.d0) then
|
|
et2 = dexp(t*t)
|
|
ss = 2d0*et2*((-2d0*t)/dsqrt(pi)+et2*(1d0+2d0*t*t)*d)
|
|
ps = (8d0*et2*t*(-2d0*(1d0+t*t)+et2*dsqrt(pi)*t*(3d0+2d0*t*t)*d))/(3d0*dsqrt(pi))
|
|
else
|
|
ss = 0.d0
|
|
ps = 0.d0
|
|
endif
|
|
|
|
endif
|
|
|
|
k = t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0)
|
|
|
|
! (s|s)
|
|
ss = k*ss
|
|
|
|
! (p|s)
|
|
ps = k*ps
|
|
|
|
pxs = AxBx*ps
|
|
pys = AyBy*ps
|
|
pzs = AzBz*ps
|
|
|
|
! (d|s)
|
|
ds = k*ds
|
|
|
|
dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds
|
|
dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds
|
|
dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds
|
|
|
|
dxys = AxBx*AyBy*ds
|
|
dxzs = AxBx*AzBz*ds
|
|
dyzs = AyBy*AzBz*ds
|
|
|
|
select case (sum(aGau))
|
|
case (0)
|
|
result = ss
|
|
|
|
case (1)
|
|
if (aGau(1) == 1) then
|
|
result = pxs
|
|
else if (aGau(2) == 1) then
|
|
result = pys
|
|
else if (aGau(3) == 1) then
|
|
result = pzs
|
|
endif
|
|
|
|
case (2)
|
|
if (aGau(1) == 2) then
|
|
result = dxxs
|
|
else if (aGau(2) == 2) then
|
|
result = dyys
|
|
else if (aGau(3) == 2) then
|
|
result = dzzs
|
|
else if (aGau(1)+aGau(2) == 2) then
|
|
result = dxys
|
|
else if (aGau(1)+aGau(3) == 2) then
|
|
result = dxzs
|
|
else if (aGau(2)+aGau(3) == 2) then
|
|
result = dyzs
|
|
endif
|
|
|
|
case default
|
|
stop 'GauSlaOverlap not implemented'
|
|
|
|
end select
|
|
|
|
end
|
|
!*****************************************************************************
|
|
|
|
!*****************************************************************************
|
|
subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla,result)
|
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
! Compute the kinetic energy integral between a Gaussian function
|
|
! with arbitrary angular momemtum and a s-type Slater function
|
|
END_DOC
|
|
|
|
! Input variables
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
double precision,intent(out) :: result
|
|
|
|
! Final value of the integrals
|
|
double precision :: ss,ps,ds
|
|
double precision :: pxs,pys,pzs
|
|
double precision :: dxxs,dyys,dzzs,dxys,dxzs,dyzs
|
|
|
|
double precision :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k
|
|
|
|
pi = 4d0*atan(1d0)
|
|
|
|
! calculate the length AB between the two centers
|
|
|
|
AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
|
|
AB = dsqrt(AB)
|
|
|
|
AxBx = (cGau(1)-cSla(1))/2d0
|
|
AyBy = (cGau(2)-cSla(2))/2d0
|
|
AzBz = (cGau(3)-cSla(3))/2d0
|
|
|
|
! intermediate variables
|
|
|
|
t = expSla*dsqrt(0.25d0/expGau)
|
|
u = dsqrt(expGau)*AB
|
|
|
|
if(AB > 0d0) then
|
|
|
|
! (s|s)
|
|
ss = (1d0+t*(t-u))*derfc(t-u)*dexp(2d0*t*(t-u)) - (1d0+t*(t+u))*derfc(t+u)*dexp(2d0*t*(t+u))
|
|
|
|
! (p|s)
|
|
ps = (dexp(t**2-2d0*t*u-u**2)*(4d0*dexp(2d0*t*u)*(1d0+t**2) &
|
|
+ dsqrt(pi)*t*(-(dexp(t**2+u**2)*(3d0+2d0*t*(t-u))*derfc(t-u)) &
|
|
- dexp(2d0*t*u+(t+u)**2)*(3d0+2d0*t*(t+u))*derfc(t+u))))/dsqrt(pi)
|
|
|
|
! (d|s)
|
|
ds = (-8d0*dexp(t**2-u**2)*u+4d0*dexp(2d0*t*(t-u))*dsqrt(pi)*t**2*((2d0+t**2-t*u)*derfc(t-u) &
|
|
- dexp(4d0*t*u)*(2d0+t*(t+u))*derfc(t+u)))/dsqrt(pi)
|
|
|
|
! backward scaling
|
|
ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0
|
|
ps = ps/u**2-ss/u**3d0
|
|
ss = ss/u
|
|
|
|
else
|
|
|
|
! concentric case
|
|
ss = (4d0*dexp(t**2)*(1d0+t**2))/dsqrt(pi)-2d0*dexp(2d0*t**2)*t*(3d0+2d0*t**2)*derfc(t)
|
|
ps = (8d0*dexp(t**2)*(-1d0+4d0*t**2+2d0*t**4d0-dexp(t**2)*dsqrt(pi)*t**3d0*(5d0+2d0*t**2)*derfc(t)))/(3d0*dsqrt(pi))
|
|
|
|
endif
|
|
|
|
k = expSla*dsqrt(expGau)*t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0)
|
|
|
|
! (s|s)
|
|
ss = k*ss
|
|
|
|
! (p|s)
|
|
ps = k*ps
|
|
|
|
pxs = AxBx*ps
|
|
pys = AyBy*ps
|
|
pzs = AzBz*ps
|
|
|
|
! (d|s)
|
|
ds = k*ds
|
|
|
|
dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds
|
|
dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds
|
|
dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds
|
|
|
|
dxys = AxBx*AyBy*ds
|
|
dxzs = AxBx*AzBz*ds
|
|
dyzs = AyBy*AzBz*ds
|
|
|
|
select case (sum(aGau))
|
|
case (0)
|
|
result = ss
|
|
|
|
case (1)
|
|
if (aGau(1) == 1) then
|
|
result = pxs
|
|
else if (aGau(2) == 1) then
|
|
result = pys
|
|
else if (aGau(3) == 1) then
|
|
result = pzs
|
|
endif
|
|
|
|
case (2)
|
|
if (aGau(1) == 2) then
|
|
result = dxxs
|
|
else if (aGau(2) == 2) then
|
|
result = dyys
|
|
else if (aGau(3) == 2) then
|
|
result = dzzs
|
|
else if (aGau(1)+aGau(2) == 2) then
|
|
result = dxys
|
|
else if (aGau(1)+aGau(3) == 2) then
|
|
result = dxzs
|
|
else if (aGau(2)+aGau(3) == 2) then
|
|
result = dyzs
|
|
endif
|
|
|
|
case default
|
|
stop 'GauSlaOverlap not implemented'
|
|
|
|
end select
|
|
|
|
end
|
|
!*****************************************************************************
|
|
|
|
|
|
|
|
!*****************************************************************************
|
|
subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result)
|
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
! Compute the nuclear attraction integral between a Gaussian function
|
|
! with arbitrary angular momemtum and a s-type Slater function
|
|
END_DOC
|
|
|
|
! Input variables
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
double precision,intent(in) :: cNuc(3)
|
|
double precision,intent(in) :: ZNuc
|
|
double precision,intent(out) :: result
|
|
|
|
! Final value of the overlap integral
|
|
double precision :: ss,ps,ds,fs
|
|
double precision :: pxs,pys,pzs
|
|
|
|
double precision :: pi,E,AB,x,y,k
|
|
|
|
pi = 4d0*atan(1d0)
|
|
E = exp(1d0)
|
|
|
|
! calculate the length AB between the two centers
|
|
|
|
AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
|
|
AB = dsqrt(AB)
|
|
|
|
! intermediate variables
|
|
|
|
x = dsqrt(expSla**2/(4d0*expGau))
|
|
y = dsqrt(expGau)*AB
|
|
|
|
if(AB > 0d0) then
|
|
ss = (1d0+x*(x+y))*derfc(x+y)*dexp(2d0*x*(x+y)) - (1d0+x*(x-y))*derfc(x-y)*dexp(2d0*x*(x-y))
|
|
ss = ss/y
|
|
else
|
|
ss = (4d0*E**x**2*(1d0+x**2))/dsqrt(Pi)-2d0*E**(2d0*x**2)*x*(3d0+2d0*x**2)*dErfc(x)
|
|
endif
|
|
|
|
k = expSla*dsqrt(expGau)*x**3d0*dexp(-x*x)*4d0*pi/expSla**(3d0/2d0)
|
|
ss = k*ss
|
|
|
|
! Print result
|
|
! write(*,*) ss
|
|
result = 0.d0
|
|
|
|
end
|
|
!*****************************************************************************
|
|
|
|
double precision function BoysF0(t)
|
|
implicit none
|
|
double precision, intent(in) :: t
|
|
double precision :: pi
|
|
|
|
pi = 4d0*atan(1d0)
|
|
|
|
if(t > 0d0) then
|
|
BoysF0 = 0.5d0*dsqrt(pi/t)*derf(dsqrt(t))
|
|
else
|
|
BoysF0 = 1d0
|
|
endif
|
|
|
|
end
|
|
!*****************************************************************************
|
|
|
|
!TODO
|
|
subroutine GauSlaOverlap_write(expGau,cGau,aGau,expSla,cSla,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
write(iunit, *) &
|
|
'SDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,', {',cSla(1),',',cSla(2),',',cSla(3),'} } ],'
|
|
result = 0.d0
|
|
end
|
|
|
|
subroutine GauSlaOverlap_read(expGau,cGau,aGau,expSla,cSla,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
read(iunit, *) result
|
|
end
|
|
|
|
subroutine GauSlaKinetic_write(expGau,cGau,aGau,expSla,cSla,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
write(iunit, *) &
|
|
'TDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} } ],'
|
|
result = 0.d0
|
|
end
|
|
|
|
subroutine GauSlaKinetic_read(expGau,cGau,aGau,expSla,cSla,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
read(iunit, *) result
|
|
end
|
|
|
|
subroutine GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
double precision,intent(in) :: cNuc(3)
|
|
double precision,intent(in) :: ZNuc
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
write(iunit, *) &
|
|
'VDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{ ', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} }, {', ZNuc, ',{', cNuc(1),',', cNuc(2),',', cNuc(3), '} } ],'
|
|
result = 0.d0
|
|
end
|
|
|
|
subroutine GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit)
|
|
implicit none
|
|
double precision,intent(in) :: expGau,expSla
|
|
double precision,intent(in) :: cGau(3),cSla(3)
|
|
integer,intent(in) :: aGau(3)
|
|
double precision,intent(in) :: cNuc(3)
|
|
double precision,intent(in) :: ZNuc
|
|
integer,intent(in) :: iunit
|
|
double precision,intent(out) :: result
|
|
read(iunit, *) result
|
|
end
|
|
!TODO
|
|
|
|
BEGIN_TEMPLATE
|
|
|
|
BEGIN_PROVIDER [ double precision, GauSla$X_matrix, (ao_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <Gaussian | Slater> overlap matrix
|
|
END_DOC
|
|
integer :: i,j,k
|
|
double precision :: cGau(3)
|
|
double precision :: cSla(3)
|
|
double precision :: expSla, res, expGau
|
|
integer :: aGau(3)
|
|
|
|
!TODO
|
|
! logical :: read
|
|
! integer :: iunit
|
|
! integer :: getunitandopen
|
|
!
|
|
! inquire(FILE=trim(ezfio_filename)//'/work/GauSla$X.dat',EXIST=read)
|
|
! if (read) then
|
|
! print *, 'READ $X'
|
|
! iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.dat','r')
|
|
! else
|
|
! print *, 'WRITE $X'
|
|
! iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.inp','w')
|
|
! write(iunit,*) '{'
|
|
! endif
|
|
!TODO
|
|
|
|
do k=1,nucl_num
|
|
cSla(1:3) = nucl_coord_transp(1:3,k)
|
|
expSla = slater_expo(k)
|
|
|
|
do i=1,ao_num
|
|
cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i))
|
|
aGau(1:3) = ao_power(i,1:3)
|
|
GauSla$X_matrix(i,k) = 0.d0
|
|
|
|
do j=1,ao_prim_num(i)
|
|
expGau = ao_expo_ordered_transp(j,i)
|
|
call GauSla$X(expGau,cGau,aGau,expSla,cSla,res)
|
|
! if (read) then
|
|
! call GauSla$X_read(expGau,cGau,aGau,expSla,cSla,res,iunit)
|
|
! else
|
|
! call GauSla$X_write(expGau,cGau,aGau,expSla,cSla,res,iunit)
|
|
! endif
|
|
GauSla$X_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
|
|
enddo
|
|
|
|
enddo
|
|
|
|
enddo
|
|
! if (.not.read) then
|
|
! write(iunit,*) '0.}'
|
|
! endif
|
|
! close(iunit)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <MO | Slater>
|
|
END_DOC
|
|
call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0, &
|
|
mo_coef, size(mo_coef,1), &
|
|
GauSla$X_matrix, size(GauSla$X_matrix,1), &
|
|
0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1))
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <AO_ortho | Slater>
|
|
END_DOC
|
|
call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, &
|
|
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
|
|
GauSla$X_matrix, size(GauSla$X_matrix,1), &
|
|
0.d0, AO_orthoSla$X_matrix, size(AO_orthoSla$X_matrix,1))
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
SUBST [ X ]
|
|
|
|
Overlap ;;
|
|
Kinetic ;;
|
|
|
|
END_TEMPLATE
|
|
|
|
BEGIN_PROVIDER [ double precision, GauSlaNuclear_matrix, (ao_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <Gaussian | Slater> overlap matrix
|
|
END_DOC
|
|
integer :: i,j,k,A
|
|
double precision :: cGau(3)
|
|
double precision :: cSla(3)
|
|
double precision :: expSla, res, expGau, Znuc, cNuc(3)
|
|
integer :: aGau(3)
|
|
|
|
!TODO
|
|
logical :: read
|
|
integer :: iunit
|
|
integer :: getunitandopen
|
|
|
|
inquire(FILE=trim(ezfio_filename)//'/work/GauSlaNuclear.dat',EXIST=read)
|
|
if (read) then
|
|
print *, 'READ Nuclear'
|
|
iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.dat','r')
|
|
else
|
|
print *, 'WRITE Nuclear'
|
|
iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.inp','w')
|
|
write(iunit,*)'{'
|
|
endif
|
|
!TODO
|
|
|
|
do k=1,nucl_num
|
|
cSla(1:3) = nucl_coord_transp(1:3,k)
|
|
expSla = slater_expo(k)
|
|
|
|
do i=1,ao_num
|
|
cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i))
|
|
aGau(1:3) = ao_power(i,1:3)
|
|
GauSlaNuclear_matrix(i,k) = 0.d0
|
|
|
|
do j=1,ao_prim_num(i)
|
|
expGau = ao_expo_ordered_transp(j,i)
|
|
do A=1,nucl_num
|
|
cNuc(1:3) = nucl_coord_transp(1:3,A)
|
|
Znuc = nucl_charge(A)
|
|
! call GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res)
|
|
if (read) then
|
|
call GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit)
|
|
else
|
|
call GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit)
|
|
endif
|
|
GauSlaNuclear_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
|
|
enddo
|
|
enddo
|
|
|
|
enddo
|
|
|
|
enddo
|
|
if (.not.read) then
|
|
write(iunit,*) '0.}'
|
|
endif
|
|
close(iunit)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, GauSlaH_matrix, (ao_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Core hamiltonian in AO basis
|
|
END_DOC
|
|
GauSlaH_matrix(1:ao_num,1:nucl_num) = &
|
|
GauSlaKinetic_matrix(1:ao_num,1:nucl_num) + &
|
|
GauSlaNuclear_matrix(1:ao_num,1:nucl_num)
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, MOSlaNuclear_matrix, (mo_tot_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <MO | Slater>
|
|
END_DOC
|
|
call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, &
|
|
mo_coef_transp, size(mo_coef_transp,1), &
|
|
GauSlaNuclear_matrix, size(GauSlaNuclear_matrix,1), &
|
|
0.d0, MOSlaNuclear_matrix, size(MOSlaNuclear_matrix,1))
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! <AO ortho | Slater>
|
|
END_DOC
|
|
call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, &
|
|
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
|
|
GauSlaH_matrix, size(GauSlaH_matrix,1), &
|
|
0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1))
|
|
END_PROVIDER
|
|
|