Skip to content

Commit

Permalink
attrashellopt / efield_shell_pair moved to attrashellopt.fh
Browse files Browse the repository at this point in the history
  • Loading branch information
@etiennepalos
etiennepalos committed Jun 19, 2024
1 parent fa2de18 commit b6ec8ac
Show file tree
Hide file tree
Showing 4 changed files with 326 additions and 176 deletions.
151 changes: 151 additions & 0 deletions src/modules/include/attrashell.fh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,3 +135,154 @@ subroutine attrashell(IIsh, JJsh)
#elif defined (OEI)
end subroutine attrashell
#endif

#if defined (OEPROP)
subroutine efield_shell_pair(IIsh,JJsh,efield_electronic)
#elif defined (OEI)
subroutine attrashellopt(IIsh,JJsh)
#endif
use allmod
use quick_overlap_module, only: opf, overlap
! use xiaoconstants
#ifdef MPIV
use mpi
#endif
implicit double precision(a-h,o-z)
dimension aux(0:20)
double precision AA(3),BB(3),CC(3),PP(3)
common /xiaoattra/attra,aux,AA,BB,CC,PP,g

double precision RA(3),RB(3),RP(3), valopf, g_table(200)

#if defined (OEPROP)
double precision, intent(inout) :: efield_electronic(3,quick_molspec%nextpoint)
#endif

Ax=xyz(1,quick_basis%katom(IIsh))
Ay=xyz(2,quick_basis%katom(IIsh))
Az=xyz(3,quick_basis%katom(IIsh))

Bx=xyz(1,quick_basis%katom(JJsh))
By=xyz(2,quick_basis%katom(JJsh))
Bz=xyz(3,quick_basis%katom(JJsh))

! The purpose of this subroutine is to calculate the nuclear attraction
! of an electron distributed between gtfs with orbital exponents a
! and b on A and B with angular momentums defined by i,j,k (a's x, y
! and z exponents, respectively) and ii,jj,k and kk on B with the core at
! (Cx,Cy,Cz) with charge Z. m is the "order" of the integral which
! arises from the recusion relationship.

! The this is taken from the recursive relation found in Obara and Saika,
! J. Chem. Phys. 84 (7) 1986, 3963.

! The first step is generating all the necessary auxillary integrals.
! These are (0|1/rc|0)^(m) = 2 Sqrt (g/Pi) (0||0) Fm(g(Rpc)^2)
! The values of m range from 0 to i+j+k+ii+jj+kk.

NII2=quick_basis%Qfinal(IIsh)
NJJ2=quick_basis%Qfinal(JJsh)
Maxm=NII2+NJJ2+1+1

do ips=1,quick_basis%kprim(IIsh)
a=quick_basis%gcexpo(ips,quick_basis%ksumtype(IIsh))
do jps=1,quick_basis%kprim(JJsh)
b=quick_basis%gcexpo(jps,quick_basis%ksumtype(JJsh))

valopf = opf(a, b, quick_basis%gccoeff(ips,quick_basis%ksumtype(IIsh)),&
quick_basis%gccoeff(jps,quick_basis%ksumtype(JJsh)), Ax, Ay, Az, Bx, By, Bz)

if(abs(valopf) .gt. quick_method%coreIntegralCutoff) then

g = a+b
Px = (a*Ax + b*Bx)/g
Py = (a*Ay + b*By)/g
Pz = (a*Az + b*Bz)/g
g_table = g**(-1.5)

constant = overlap(a,b,0,0,0,0,0,0,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,g_table) &
* 2.d0 * sqrt(g/Pi)

#if defined (OEPROP)
do igridpoint=1,quick_molspec%nextpoint
Cx=quick_molspec%extxyz(1,igridpoint)
Cy=quick_molspec%extxyz(2,igridpoint)
Cz=quick_molspec%extxyz(3,igridpoint)
PCsquare = (Px-Cx)**2 + (Py -Cy)**2 + (Pz -Cz)**2
U = g* PCsquare
! Maxm = i+j+k+ii+jj+kk
call FmT(Maxm,U,aux)
do L = 0,maxm
aux(L) = -1.0d0*aux(L)*constant
attraxiao(1,1,L)=aux(L)
enddo
do L = 0,maxm-1
attraxiaoopt(1,1,1,L)=2.0d0*g*(Px-Cx)*aux(L+1)
attraxiaoopt(2,1,1,L)=2.0d0*g*(Py-Cy)*aux(L+1)
attraxiaoopt(3,1,1,L)=2.0d0*g*(Pz-Cz)*aux(L+1)
enddo
! At this point all the auxillary integrals have been calculated.
! It is now time to decompase the attraction integral to it's
! auxillary integrals through the recursion scheme. To do this we use
! a recursive function.
! attraction = attrecurse(i,j,k,ii,jj,kk,0,aux,Ax,Ay,Az,Bx,By,Bz, &
! Cx,Cy,Cz,Px,Py,Pz,g)
NIJ1=10*NII2+NJJ2
call efield_1pdm(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &
Cx,Cy,Cz,Px,Py,Pz,efield_electronic(1,igridpoint))
!endif
#elif defined (OEI)
do iatom=1,natom+quick_molspec%nextatom
if(quick_basis%katom(IIsh).eq.iatom.and.quick_basis%katom(JJsh).eq.iatom)then
continue
else
if(iatom<=natom)then
Cx=xyz(1,iatom)
Cy=xyz(2,iatom)
Cz=xyz(3,iatom)
Z=-1.0d0*quick_molspec%chg(iatom)
else
Cx=quick_molspec%extxyz(1,iatom-natom)
Cy=quick_molspec%extxyz(2,iatom-natom)
Cz=quick_molspec%extxyz(3,iatom-natom)
Z=-1.0d0*quick_molspec%extchg(iatom-natom)
endif
PCsquare = (Px-Cx)**2 + (Py -Cy)**2 + (Pz -Cz)**2
U = g* PCsquare
! Maxm = i+j+k+ii+jj+kk
call FmT(Maxm,U,aux)
do L = 0,maxm
aux(L) = aux(L)*constant*Z
attraxiao(1,1,L)=aux(L)
enddo
do L = 0,maxm-1
attraxiaoopt(1,1,1,L)=2.0d0*g*(Px-Cx)*aux(L+1)
attraxiaoopt(2,1,1,L)=2.0d0*g*(Py-Cy)*aux(L+1)
attraxiaoopt(3,1,1,L)=2.0d0*g*(Pz-Cz)*aux(L+1)
enddo
! At this point all the auxillary integrals have been calculated.
! It is now time to decompase the attraction integral to it's
! auxillary integrals through the recursion scheme. To do this we use
! a recursive function.
! attraction = attrecurse(i,j,k,ii,jj,kk,0,aux,Ax,Ay,Az,Bx,By,Bz, &
! Cx,Cy,Cz,Px,Py,Pz,g)
NIJ1=10*NII2+NJJ2
call nuclearattraopt(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &
Cx,Cy,Cz,Px,Py,Pz,iatom)
endif
#endif
enddo
endif
enddo
enddo

! Xiao HE remember to multiply Z 01/12/2008
! attraction = attraction*(-1.d0)* Z

#if defined (OEPROP)
return
end subroutine efield_shell_pair
#elif defined (OEI)
return
end subroutine attrashellopt
#endif
3 changes: 1 addition & 2 deletions src/modules/include/nuclearattra.fh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1172,7 +1172,6 @@ subroutine nuclearattraopt(Ips,Jps,IIsh,JJsh,NIJ1, &
efield(1) = efield(1)- CGrad1
efield(2) = efield(2)- CGrad2
efield(3) = efield(3)- CGrad3

#elif defined (OEI)
itemp1new=trans(quick_basis%KLMN(1,III)+1,quick_basis%KLMN(2,III),quick_basis%KLMN(3,III))
Agrad1=Agrad1+2.0d0*Xconstant2* &
Expand Down Expand Up @@ -1259,4 +1258,4 @@ endif
End subroutine efield_1pdm
#elif defined (OEI)
End subroutine nuclearattraopt
#endif
#endif
184 changes: 92 additions & 92 deletions src/modules/quick_oei_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -445,122 +445,122 @@ end subroutine kineticO



subroutine attrashellopt(IIsh,JJsh)
use allmod
use quick_overlap_module, only: opf, overlap
! use xiaoconstants
#ifdef MPIV
use mpi
#endif
implicit double precision(a-h,o-z)
dimension aux(0:20)
double precision AA(3),BB(3),CC(3),PP(3)
common /xiaoattra/attra,aux,AA,BB,CC,PP,g
! subroutine attrashellopt(IIsh,JJsh)
! use allmod
! use quick_overlap_module, only: opf, overlap
! ! use xiaoconstants
! #ifdef MPIV
! use mpi
! #endif
! implicit double precision(a-h,o-z)
! dimension aux(0:20)
! double precision AA(3),BB(3),CC(3),PP(3)
! common /xiaoattra/attra,aux,AA,BB,CC,PP,g

double precision RA(3),RB(3),RP(3), valopf, g_table(200)
! double precision RA(3),RB(3),RP(3), valopf, g_table(200)

Ax=xyz(1,quick_basis%katom(IIsh))
Ay=xyz(2,quick_basis%katom(IIsh))
Az=xyz(3,quick_basis%katom(IIsh))
! Ax=xyz(1,quick_basis%katom(IIsh))
! Ay=xyz(2,quick_basis%katom(IIsh))
! Az=xyz(3,quick_basis%katom(IIsh))

Bx=xyz(1,quick_basis%katom(JJsh))
By=xyz(2,quick_basis%katom(JJsh))
Bz=xyz(3,quick_basis%katom(JJsh))
! Bx=xyz(1,quick_basis%katom(JJsh))
! By=xyz(2,quick_basis%katom(JJsh))
! Bz=xyz(3,quick_basis%katom(JJsh))

! The purpose of this subroutine is to calculate the nuclear attraction
! of an electron distributed between gtfs with orbital exponents a
! and b on A and B with angular momentums defined by i,j,k (a's x, y
! and z exponents, respectively) and ii,jj,k and kk on B with the core at
! (Cx,Cy,Cz) with charge Z. m is the "order" of the integral which
! arises from the recusion relationship.
! ! The purpose of this subroutine is to calculate the nuclear attraction
! ! of an electron distributed between gtfs with orbital exponents a
! ! and b on A and B with angular momentums defined by i,j,k (a's x, y
! ! and z exponents, respectively) and ii,jj,k and kk on B with the core at
! ! (Cx,Cy,Cz) with charge Z. m is the "order" of the integral which
! ! arises from the recusion relationship.

! The this is taken from the recursive relation found in Obara and Saika,
! J. Chem. Phys. 84 (7) 1986, 3963.
! ! The this is taken from the recursive relation found in Obara and Saika,
! ! J. Chem. Phys. 84 (7) 1986, 3963.

! The first step is generating all the necessary auxillary integrals.
! These are (0|1/rc|0)^(m) = 2 Sqrt (g/Pi) (0||0) Fm(g(Rpc)^2)
! The values of m range from 0 to i+j+k+ii+jj+kk.
! ! The first step is generating all the necessary auxillary integrals.
! ! These are (0|1/rc|0)^(m) = 2 Sqrt (g/Pi) (0||0) Fm(g(Rpc)^2)
! ! The values of m range from 0 to i+j+k+ii+jj+kk.

NII2=quick_basis%Qfinal(IIsh)
NJJ2=quick_basis%Qfinal(JJsh)
Maxm=NII2+NJJ2+1+1
! NII2=quick_basis%Qfinal(IIsh)
! NJJ2=quick_basis%Qfinal(JJsh)
! Maxm=NII2+NJJ2+1+1

do ips=1,quick_basis%kprim(IIsh)
a=quick_basis%gcexpo(ips,quick_basis%ksumtype(IIsh))
do jps=1,quick_basis%kprim(JJsh)
b=quick_basis%gcexpo(jps,quick_basis%ksumtype(JJsh))
! do ips=1,quick_basis%kprim(IIsh)
! a=quick_basis%gcexpo(ips,quick_basis%ksumtype(IIsh))
! do jps=1,quick_basis%kprim(JJsh)
! b=quick_basis%gcexpo(jps,quick_basis%ksumtype(JJsh))

valopf = opf(a, b, quick_basis%gccoeff(ips,quick_basis%ksumtype(IIsh)),&
quick_basis%gccoeff(jps,quick_basis%ksumtype(JJsh)), Ax, Ay, Az, Bx, By, Bz)
! valopf = opf(a, b, quick_basis%gccoeff(ips,quick_basis%ksumtype(IIsh)),&
! quick_basis%gccoeff(jps,quick_basis%ksumtype(JJsh)), Ax, Ay, Az, Bx, By, Bz)

if(abs(valopf) .gt. quick_method%coreIntegralCutoff) then
! if(abs(valopf) .gt. quick_method%coreIntegralCutoff) then

g = a+b
Px = (a*Ax + b*Bx)/g
Py = (a*Ay + b*By)/g
Pz = (a*Az + b*Bz)/g
g_table = g**(-1.5)
! g = a+b
! Px = (a*Ax + b*Bx)/g
! Py = (a*Ay + b*By)/g
! Pz = (a*Az + b*Bz)/g
! g_table = g**(-1.5)

constant = overlap(a,b,0,0,0,0,0,0,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,g_table) &
* 2.d0 * sqrt(g/Pi)
! constant = overlap(a,b,0,0,0,0,0,0,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,g_table) &
! * 2.d0 * sqrt(g/Pi)

do iatom=1,natom+quick_molspec%nextatom
if(quick_basis%katom(IIsh).eq.iatom.and.quick_basis%katom(JJsh).eq.iatom)then
continue
else
if(iatom<=natom)then
Cx=xyz(1,iatom)
Cy=xyz(2,iatom)
Cz=xyz(3,iatom)
Z=-1.0d0*quick_molspec%chg(iatom)
else
Cx=quick_molspec%extxyz(1,iatom-natom)
Cy=quick_molspec%extxyz(2,iatom-natom)
Cz=quick_molspec%extxyz(3,iatom-natom)
Z=-1.0d0*quick_molspec%extchg(iatom-natom)
endif
! do iatom=1,natom+quick_molspec%nextatom
! if(quick_basis%katom(IIsh).eq.iatom.and.quick_basis%katom(JJsh).eq.iatom)then
! continue
! else
! if(iatom<=natom)then
! Cx=xyz(1,iatom)
! Cy=xyz(2,iatom)
! Cz=xyz(3,iatom)
! Z=-1.0d0*quick_molspec%chg(iatom)
! else
! Cx=quick_molspec%extxyz(1,iatom-natom)
! Cy=quick_molspec%extxyz(2,iatom-natom)
! Cz=quick_molspec%extxyz(3,iatom-natom)
! Z=-1.0d0*quick_molspec%extchg(iatom-natom)
! endif

PCsquare = (Px-Cx)**2 + (Py -Cy)**2 + (Pz -Cz)**2
! PCsquare = (Px-Cx)**2 + (Py -Cy)**2 + (Pz -Cz)**2

U = g* PCsquare
! Maxm = i+j+k+ii+jj+kk
call FmT(Maxm,U,aux)
do L = 0,maxm
aux(L) = aux(L)*constant*Z
attraxiao(1,1,L)=aux(L)
enddo
! U = g* PCsquare
! ! Maxm = i+j+k+ii+jj+kk
! call FmT(Maxm,U,aux)
! do L = 0,maxm
! aux(L) = aux(L)*constant*Z
! attraxiao(1,1,L)=aux(L)
! enddo

do L = 0,maxm-1
attraxiaoopt(1,1,1,L)=2.0d0*g*(Px-Cx)*aux(L+1)
attraxiaoopt(2,1,1,L)=2.0d0*g*(Py-Cy)*aux(L+1)
attraxiaoopt(3,1,1,L)=2.0d0*g*(Pz-Cz)*aux(L+1)
enddo
! do L = 0,maxm-1
! attraxiaoopt(1,1,1,L)=2.0d0*g*(Px-Cx)*aux(L+1)
! attraxiaoopt(2,1,1,L)=2.0d0*g*(Py-Cy)*aux(L+1)
! attraxiaoopt(3,1,1,L)=2.0d0*g*(Pz-Cz)*aux(L+1)
! enddo

! At this point all the auxillary integrals have been calculated.
! It is now time to decompase the attraction integral to it's
! auxillary integrals through the recursion scheme. To do this we use
! a recursive function.
! ! At this point all the auxillary integrals have been calculated.
! ! It is now time to decompase the attraction integral to it's
! ! auxillary integrals through the recursion scheme. To do this we use
! ! a recursive function.

! attraction = attrecurse(i,j,k,ii,jj,kk,0,aux,Ax,Ay,Az,Bx,By,Bz, &
! ! attraction = attrecurse(i,j,k,ii,jj,kk,0,aux,Ax,Ay,Az,Bx,By,Bz, &

! Cx,Cy,Cz,Px,Py,Pz,g)
NIJ1=10*NII2+NJJ2
! ! Cx,Cy,Cz,Px,Py,Pz,g)
! NIJ1=10*NII2+NJJ2

call nuclearattraopt(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &
! call nuclearattraopt(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &

Cx,Cy,Cz,Px,Py,Pz,iatom)
! Cx,Cy,Cz,Px,Py,Pz,iatom)

endif
! endif

enddo
endif
enddo
enddo
! enddo
! endif
! enddo
! enddo

! Xiao HE remember to multiply Z 01/12/2008
! attraction = attraction*(-1.d0)* Z
return
end subroutine attrashellopt
! ! Xiao HE remember to multiply Z 01/12/2008
! ! attraction = attraction*(-1.d0)* Z
! return
! end subroutine attrashellopt


double precision function ekinetic(a,b,i,j,k,ii,jj,kk,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,g_table)
Expand Down
Loading

0 comments on commit b6ec8ac

Please sign in to comment.