Skip to content

Commit

Permalink
solved merge confict. Added headed files attrashell and nuclearattra
Browse files Browse the repository at this point in the history
  • Loading branch information
@vtripath65
vtripath65 committed Jun 5, 2024
2 parents 882c427 + ef6ff5b commit ad42ae7
Show file tree
Hide file tree
Showing 6 changed files with 807 additions and 726 deletions.
137 changes: 137 additions & 0 deletions src/modules/include/attrashell.fh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
!-----------------------------------------------------------------------------------------!
! This subroutine computes the V_elec contribution for a shell pair for each grid point !
! It loops over each gridpoint, calls esp_1pdm and stores the value in esp_electronic() !
! This is - \sum_{mu nu} P_{mu nu} * V_{mu nu} !
! See Eqn. A14 of Obara-Saika [J. Chem. Phys. 84, 3963 (1986)] !
! First, calculates 〈 phi_mu | phi_nu 〉 for all mu and nu !
! Then, P_{mu nu} * 〈 phi_mu | 1/|r-C| | phi_nu 〉 !
!-----------------------------------------------------------------------------------------!
#if defined (OEPROP)
subroutine esp_shell_pair(IIsh, JJsh, esp_electronic)
#elif defined (OEI)
subroutine attrashell(IIsh, JJsh)
#endif
use quick_method_module, only: quick_method
use quick_basis_module, only: quick_basis, attraxiao
use quick_molspec_module, only: quick_molspec, xyz
use quick_overlap_module, only: gpt, opf, overlap_core
use quick_constants_module, only : Pi

#if defined (OEI)
use quick_molspec_module, only: natom
#endif

implicit none

integer :: IIsh, JJsh, ips, jps, L, Maxm, NII2, NIJ1, NJJ2
double precision :: a, b, Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz, g, U
double precision :: attra, const, constanttemp, PCsquare, Px, Py, Pz

#if defined (OEI)
double precision :: Z
#endif

double precision, dimension(0:20) :: aux
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),inv_g,g_table(200), valopf

#if defined (OEPROP)
integer :: igridpoint
double precision, dimension(:), intent(inout) :: esp_electronic
#elif defined (OEI)
integer :: iatom
#endif

! Related to positions of "QM" atoms
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))

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

! Calculation of V_elec starts here
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

!Eqn 14 Obara-Saika 86
call gpt(a,b,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,0,g_table)
g = a+b
!Eqn 15 Obara-Saika 86
inv_g = 1.0d0 / dble(g)

! Calculate first two terms of Obara & Saika Eqn A20
constanttemp=dexp(-((a*b*((Ax - Bx)**2.d0 + (Ay - By)**2.d0 + (Az - Bz)**2.d0))*inv_g))
const = overlap_core(a,b,0,0,0,0,0,0,Ax,Ay,Az,Bx,By,Bz,Px,Py,Pz,g_table) * 2.d0 * sqrt(g/Pi)*constanttemp

#if defined (OEPROP)
! Loops over external grid points/MM atoms
do igridpoint=1,quick_molspec%nextpoint
Cx=quick_molspec%extxyz(1,igridpoint)
Cy=quick_molspec%extxyz(2,igridpoint)
Cz=quick_molspec%extxyz(3,igridpoint)
#elif defined (OEI)
!nextatom=number of external MM point charges. set to 0 if none used
do iatom=1,natom+quick_molspec%nextatom
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=-quick_molspec%extchg(iatom-natom)
endif
#endif
! Calculate the last term of Obara--Saika Eqn A21
PCsquare = (Px-Cx)**2 + (Py -Cy)**2 + (Pz -Cz)**2
! Compute Obara--Saika Eqn A21
U = g* PCsquare

! Calculate the last term of Obara--Saika Eqn A20
call FmT(Maxm,U,aux)
! Calculate all the auxilary integrals and store in attraxiao array
do L = 0,maxm
#if defined (OEPROP)
! sign (-1.0d0) is used to ensure the auxilary integrals are negative
aux(L) = -1.0d0*aux(L)*const
#elif defined (OEI)
aux(L) = aux(L)*const*Z
#endif
attraxiao(1,1,L)=aux(L)
enddo

NIJ1=10*NII2+NJJ2
#if defined (OEPROP)
! Call and get P_{mu nu} V_{mu nu} into esp_electronic( )
call esp_1pdm(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &
Cx,Cy,Cz,Px,Py,Pz, esp_electronic(igridpoint))
#elif defined (OEI)
call nuclearattra(ips,jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz, &
Cx,Cy,Cz,Px,Py,Pz)
#endif
enddo

endif
enddo
enddo
#if defined (OEPROP)
end subroutine esp_shell_pair
#elif defined (OEI)
end subroutine attrashell
#endif
258 changes: 258 additions & 0 deletions src/modules/include/nuclearattra.fh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,258 @@
#if defined (OEPROP)
subroutine esp_1pdm(Ips,Jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz,Px,Py,Pz,esp)
#elif defined (OEI)
subroutine nuclearattra(Ips,Jps,IIsh,JJsh,NIJ1,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz,Px,Py,Pz)
#endif
use quick_params_module, only: trans
use quick_calculated_module, only : quick_qm_struct
use quick_basis_module, only: attraxiao,quick_basis
use quick_files_module, only: ioutfile

#if defined (OEPROP)
use quick_method_module, only: quick_method
#endif

implicit none

double precision attra,aux(0:20)
integer a(3),b(3)
integer Ips, Jps, IIsh, JJsh, NIJ1
integer Iang, III, III1, III2, itemp, itemp1, itemp2, itempt
integer Jang, JJJ, JJJ1, JJJ2, NBI1, NBI2, NBJ1, NBJ2
double precision X1temp, Xconstant, dense_sym_factor
double precision Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz,Px,Py,Pz,g
double precision AA(3),BB(3),CC(3),PP(3)

#if defined (OEPROP)
double precision, intent(inout) :: esp
#endif

common /xiaoattra/attra,aux,AA,BB,CC,PP,g

AA(1)=Ax
AA(2)=Ay
AA(3)=Az
BB(1)=Bx
BB(2)=By
BB(3)=Bz
CC(1)=Cx
CC(2)=Cy
CC(3)=Cz
PP(1)=Px
PP(2)=Py
PP(3)=Pz

select case (NIJ1)

case(0)
case(10)
call PSattra(0)
case(1)
call SPattra(0)
case(11)
call SPattra(0)
call PSattra(0)
call PSattra(1)
call PPattra(0)
case(20)
call PSattra(0)
call PSattra(1)
call DSattra(0)
case(2)
call SPattra(0)
call SPattra(1)
call SDattra(0)
case(21)
call PSattra(0)
call PSattra(1)
call PSattra(2)
call DSattra(0)
call DSattra(1)
call DPattra(0)
case(12)
call SPattra(0)
call SPattra(1)
call SPattra(2)
call SDattra(0)
call SDattra(1)
call PDattra(0)
case(22)
do itempt=0,3
call PSattra(itempt)
enddo
do itempt=0,1
call PPattra(itempt)
enddo
do itempt=0,2
call DSattra(itempt)
enddo
do itempt=0,1
call DPattra(itempt)
enddo
call DDattra(0)
case(30)
do itemp=0,2
call PSattra(itemp)
enddo
do itemp=0,1
call DSattra(itemp)
enddo
call FSattra(0)
case(3)
do itemp=0,2
call SPattra(itemp)
enddo
do itemp=0,1
call SDattra(itemp)
enddo
call SFattra(0)
case(31)
do itemp=0,3
call PSattra(itemp)
enddo
do itemp=0,2
call PPattra(itemp)
enddo
do itemp=0,2
call DSattra(itemp)
enddo
do itemp=0,1
call DPattra(itemp)
enddo
do itemp=0,1
call FSattra(itemp)
enddo
call FPattra(0)
case(13)
do itemp=0,3
call SPattra(itemp)
call PSattra(itemp)
enddo
do itemp=0,2
call PPattra(itemp)
enddo
do itemp=0,2
call SDattra(itemp)
enddo
do itemp=0,1
call PDattra(itemp)
enddo
do itemp=0,1
call SFattra(itemp)
enddo
call PFattra(0)
case(32)
do itemp=0,4
call PSattra(itemp)
enddo
do itemp=0,3
call PPattra(itemp)
enddo
do itemp=0,3
call DSattra(itemp)
enddo
do itemp=0,2
call DPattra(itemp)
enddo
do itemp=0,2
call FSattra(itemp)
enddo
do itemp=0,1
call FPattra(itemp)
enddo
call FDattra(0)
case(23)
do itemp=0,4
call SPattra(itemp)
call PSattra(itemp)
enddo
do itemp=0,3
call PPattra(itemp)
enddo
do itemp=0,3
call SDattra(itemp)
enddo
do itemp=0,2
call PDattra(itemp)
enddo
do itemp=0,2
call SFattra(itemp)
enddo
do itemp=0,1
call PFattra(itemp)
enddo
call DFattra(0)
case(33)
do itemp=0,5
call PSattra(itemp)
enddo
do itemp=0,4
call PPattra(itemp)
enddo
do itemp=0,4
call DSattra(itemp)
enddo
do itemp=0,3
call DPattra(itemp)
enddo
do itemp=0,2
call DDattra(itemp)
enddo
do itemp=0,3
call FSattra(itemp)
enddo
do itemp=0,2
call FPattra(itemp)
enddo
do itemp=0,1
call FDattra(itemp)
enddo
call FFattra(0)
end select

do Iang=quick_basis%Qstart(IIsh),quick_basis%Qfinal(IIsh)
X1temp=quick_basis%gccoeff(ips,quick_basis%ksumtype(IIsh)+Iang)
do Jang=quick_basis%Qstart(JJsh),quick_basis%Qfinal(JJsh)
NBI1=quick_basis%Qsbasis(IIsh,Iang)
NBI2=quick_basis%Qfbasis(IIsh,Iang)
NBJ1=quick_basis%Qsbasis(JJsh,Jang)
NBJ2=quick_basis%Qfbasis(JJsh,Jang)

III1=quick_basis%ksumtype(IIsh)+NBI1
III2=quick_basis%ksumtype(IIsh)+NBI2
JJJ1=quick_basis%ksumtype(JJsh)+NBJ1
JJJ2=quick_basis%ksumtype(JJsh)+NBJ2

Xconstant=X1temp*quick_basis%gccoeff(jps,quick_basis%ksumtype(JJsh)+Jang)
do III=III1,III2

itemp1=trans(quick_basis%KLMN(1,III),quick_basis%KLMN(2,III),quick_basis%KLMN(3,III))
do JJJ=max(III,JJJ1),JJJ2
itemp2=trans(quick_basis%KLMN(1,JJJ),quick_basis%KLMN(2,JJJ),quick_basis%KLMN(3,JJJ))

#if defined (OEPROP)
dense_sym_factor = 1.0d0
if (III /= JJJ) dense_sym_factor = 2.0d0

if (quick_method%UNRST) quick_qm_struct%denseSave(JJJ,III) = (quick_qm_struct%dense(JJJ,III)&
+quick_qm_struct%denseb(JJJ,III))

esp = esp + dense_sym_factor*quick_qm_struct%denseSave(JJJ,III)*Xconstant &
*quick_basis%cons(III)*quick_basis%cons(JJJ)*attraxiao(itemp1,itemp2,0)
#elif defined (OEI)
quick_qm_struct%o(JJJ,III)=quick_qm_struct%o(JJJ,III)+ &
Xconstant*quick_basis%cons(III)*quick_basis%cons(JJJ)*attraxiao(itemp1,itemp2,0)
#endif

enddo
enddo

enddo
enddo
201 return

#if defined (OEPROP)
End subroutine esp_1pdm
#elif defined (OEI)
End subroutine nuclearattra
#endif
Loading

0 comments on commit ad42ae7

Please sign in to comment.