4
1
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00

AOtoMO OOOOOO

This commit is contained in:
Pierre-Francois Loos 2019-03-13 15:02:16 +01:00
parent 6f694db086
commit 355e5026d9
13 changed files with 157 additions and 101 deletions

View File

@ -1,9 +1,20 @@
1 3 1 5
S 3 1.00 S 6 1.00
38.3600000 0.0238090 1264.5857000 0.0019448
5.7700000 0.1548910 189.9368100 0.0148351
1.2400000 0.4699870 43.1590890 0.0720906
S 1 1.00 12.0986630 0.2371542
0.2976000 1.0000000 3.8063232 0.4691987
P 1 1.00 1.2728903 0.3565202
1.2750000 1.0000000 S 3 1.00
3.1964631 -0.1126487
0.7478133 -0.2295064
0.2199663 1.1869167
S 1 1.00
0.0823099 1.0000000
P 3 1.00
3.1964631 0.0559802
0.7478133 0.2615506
0.2199663 0.7939723
P 1 1.00
0.0823099 1.0000000

View File

@ -5,10 +5,10 @@
# Exposant of the Slater geminal # Exposant of the Slater geminal
1.0 1.0
# One-electron integrals: Ov Kin Nuc # One-electron integrals: Ov Kin Nuc
T T T F F F
# Two-electron integrals: ERI F12 Yuk Erf # Two-electron integrals: ERI F12 Yuk Erf
T F F F T F F F
# Three-electron integrals: Type1 Type2 Type3 # Three-electron integrals: Type1 Type2 Type3
F F F T F F
# Four-electron integrals: Type1 Type2 Type3 # Four-electron integrals: Type1 Type2 Type3
F F F F F F

View File

@ -1,7 +1,7 @@
# HF MOM # HF MOM
T F T F
# MP2 MP3 # MP2 MP3 MP2-F12
F F F F T
# CIS TDHF ADC # CIS TDHF ADC
F F F F F F
# GF2 GF3 # GF2 GF3

View File

@ -1,4 +1,4 @@
# nAt nEl nEla nElb # nAt nEl nCore nRyd
1 2 1 1 1 4 0 0
# Znuc x y z # Znuc x y z
He 0.0 0.0 0.0 Be 0.0 0.0 0.0

View File

@ -1,7 +1,7 @@
# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type # RHF: maxSCF thresh DIIS n_diis guess_type ortho_type
32 0.0000001 T 5 1 1 32 0.0000001 T 5 1 1
# MPn: # MP: MP2-F12
T
# CIS/TDHF: singlet triplet # CIS/TDHF: singlet triplet
T F T F
# GF: maxSCF thresh DIIS n_diis renormalization # GF: maxSCF thresh DIIS n_diis renormalization

View File

@ -1,9 +1,20 @@
1 3 1 5
S 3 1.00 S 6 1.00
38.3600000 0.0238090 1264.5857000 0.0019448
5.7700000 0.1548910 189.9368100 0.0148351
1.2400000 0.4699870 43.1590890 0.0720906
S 1 1.00 12.0986630 0.2371542
0.2976000 1.0000000 3.8063232 0.4691987
P 1 1.00 1.2728903 0.3565202
1.2750000 1.0000000 S 3 1.00
3.1964631 -0.1126487
0.7478133 -0.2295064
0.2199663 1.1869167
S 1 1.00
0.0823099 1.0000000
P 3 1.00
3.1964631 0.0559802
0.7478133 0.2615506
0.2199663 0.7939723
P 1 1.00
0.0823099 1.0000000

View File

@ -1,4 +1,4 @@
subroutine ComputeOv(debug,NBasis,nShell, & subroutine ComputeOv(debug,nBas,nShell, &
CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
npOv,nSigpOv,ncOv,nSigcOv,S) npOv,nSigpOv,ncOv,nSigcOv,S)
@ -11,7 +11,7 @@ subroutine ComputeOv(debug,NBasis,nShell, &
! Input variables ! Input variables
logical,intent(in) :: debug logical,intent(in) :: debug
integer,intent(in) :: NBasis,nShell integer,intent(in) :: nBas,nShell
double precision,intent(in) :: CenterShell(maxShell,3) double precision,intent(in) :: CenterShell(maxShell,3)
integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
@ -39,7 +39,7 @@ subroutine ComputeOv(debug,NBasis,nShell, &
! Output variables ! Output variables
integer,intent(out) :: npOv,nSigpOv,ncOv,nSigcOv integer,intent(out) :: npOv,nSigpOv,ncOv,nSigcOv
double precision,intent(out) :: S(NBasis,NBasis) double precision,intent(out) :: S(nBas,nBas)
! Compute one-electron integrals ! Compute one-electron integrals

View File

@ -19,7 +19,7 @@ program IntPak
logical :: do4eInt(n4eInt) logical :: do4eInt(n4eInt)
integer :: nNuc,nBas,iType integer :: nNuc,nBas,iType
integer :: nEl,nO,nV integer :: nEl,nO,nV,nCore,nRyd
double precision :: ExpS double precision :: ExpS
double precision :: ENuc double precision :: ENuc
integer :: KG integer :: KG
@ -69,7 +69,7 @@ program IntPak
! nBas = number of basis functions (see below) ! nBas = number of basis functions (see below)
! = nO + nV ! = nO + nV
call read_molecule(nNuc,nEl,nO) call read_molecule(nNuc,nEl,nO,nCore,nRyd)
allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3)) allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3))
@ -346,13 +346,13 @@ program IntPak
if(do3eInt(1)) then if(do3eInt(1)) then
iType = 1 iType = 1
KG = 1 ! KG = 1
! KG = 6 KG = 6
allocate(DG(1:KG),ExpG(1:KG)) allocate(DG(1:KG),ExpG(1:KG))
DG = (/ 1d0 /) ! DG = (/ 1d0 /)
ExpG = (/ 1d0 /) ! ExpG = (/ 1d0 /)
! DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /) DG = (/ 0.3144d0, 0.3037d0, 0.1681d0, 0.09811d0, 0.06024d0, 0.03726d0 /)
! ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /) ExpG = (/ 0.2209d0, 1.004d0, 3.622d0, 12.16d0, 45.87d0, 254.4d0 /)
call cpu_time(start_3eInt(iType)) call cpu_time(start_3eInt(iType))
call Compute3eInt(debug,iType,nShell, & call Compute3eInt(debug,iType,nShell, &

View File

@ -34,17 +34,35 @@ subroutine S3eInt(debug,iType,np3eInt,nSigp3eInt, &
double precision :: p3eInt double precision :: p3eInt
p3eInt = 0d0 p3eInt = 0d0
do k=1,KG
do l=1,KG if(iType == 1) then
do k=1,KG
ExpSG13 = ExpG(k)*ExpS**2 ExpSG13 = ExpG(k)*ExpS**2
ExpSG23 = ExpG(l)*ExpS**2 p3eInt = p3eInt &
p3eInt = p3eInt & + DG(k)*G3eInt(debug,iType, &
+ DG(k)*DG(l)*G3eInt(debug,iType, & ExpSG13,ExpSG23, &
ExpSG13,ExpSG23, & ExpBra,CenterBra,AngMomBra, &
ExpBra,CenterBra,AngMomBra, & ExpKet,CenterKet,AngMomKet)
ExpKet,CenterKet,AngMomKet) end do
enddo
enddo end if
if(iType == 2 .or. iType == 3) then
do k=1,KG
do l=1,KG
ExpSG13 = ExpG(k)*ExpS**2
ExpSG23 = ExpG(l)*ExpS**2
p3eInt = p3eInt &
+ DG(k)*DG(l)*G3eInt(debug,iType, &
ExpSG13,ExpSG23, &
ExpBra,CenterBra,AngMomBra, &
ExpKet,CenterKet,AngMomKet)
end do
end do
end if
! Print result ! Print result

View File

@ -10,11 +10,11 @@ program MCQC
logical :: doG0W0,doevGW,doqsGW logical :: doG0W0,doevGW,doqsGW
logical :: doMCMP2,doMinMCMP2 logical :: doMCMP2,doMinMCMP2
logical :: doeNcusp logical :: doeNcusp
integer :: nAt,nBas,nBasCABS,nEl,nC,nO,nV,nR,nS integer :: nNuc,nBas,nBasCABS,nEl,nC,nO,nV,nR,nS
double precision :: ENuc,ERHF,Norm double precision :: ENuc,ERHF,Norm
double precision :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3) double precision :: EcMP2(3),EcMP3,EcMP2F12(3),EcMCMP2(3),Err_EcMCMP2(3),Var_EcMCMP2(3)
double precision,allocatable :: ZNuc(:),rAt(:,:),cHF(:,:),eHF(:),eG0W0(:),PHF(:,:) double precision,allocatable :: ZNuc(:),rNuc(:,:),cHF(:,:),eHF(:),eG0W0(:),PHF(:,:)
integer :: nShell integer :: nShell
integer,allocatable :: TotAngMomShell(:),KShell(:) integer,allocatable :: TotAngMomShell(:),KShell(:)
@ -25,7 +25,7 @@ program MCQC
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:) double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:)
double precision,allocatable :: ERI_AO_basis(:,:,:,:),ERI_MO_basis(:,:,:,:) double precision,allocatable :: ERI_AO_basis(:,:,:,:),ERI_MO_basis(:,:,:,:)
double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
double precision :: start_HF ,end_HF ,t_HF double precision :: start_HF ,end_HF ,t_HF
double precision :: start_MOM ,end_MOM ,t_MOM double precision :: start_MOM ,end_MOM ,t_MOM
@ -77,7 +77,7 @@ program MCQC
! Which calculations do you want to do? ! Which calculations do you want to do?
call read_methods(doHF,doMOM, & call read_methods(doHF,doMOM, &
doMP2,doMP3, & doMP2,doMP3,doMP2F12, &
doCIS,doTDHF,doADC, & doCIS,doTDHF,doADC, &
doGF2,doGF3, & doGF2,doGF3, &
doG0W0,doevGW,doqsGW, & doG0W0,doevGW,doqsGW, &
@ -95,7 +95,6 @@ program MCQC
doeNCusp = .false. doeNCusp = .false.
doMinMCMP2 = .false. doMinMCMP2 = .false.
doMP2F12 = .false.
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Read input information ! Read input information
@ -109,12 +108,12 @@ program MCQC
! nBas = number of basis functions (see below) ! nBas = number of basis functions (see below)
! = nO + nV ! = nO + nV
call read_molecule(nAt,nEl,nC,nO,nR) call read_molecule(nNuc,nEl,nC,nO,nR)
allocate(ZNuc(nAt),rAt(nAt,3)) allocate(ZNuc(nNuc),rNuc(nNuc,3))
! Read geometry ! Read geometry
call read_geometry(nAt,ZNuc,rAt,ENuc) call read_geometry(nNuc,ZNuc,rNuc,ENuc)
allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), & allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), &
DShell(maxShell,maxK),ExpShell(maxShell,maxK)) DShell(maxShell,maxK),ExpShell(maxShell,maxK))
@ -123,14 +122,14 @@ program MCQC
! Read basis set information ! Read basis set information
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call read_basis(nAt,rAt,nBas,nC,nO,nV,nR,nS, & call read_basis(nNuc,rNuc,nBas,nC,nO,nV,nR,nS, &
nShell,TotAngMomShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell) nShell,TotAngMomShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Read auxiliary basis set information ! Read auxiliary basis set information
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! call ReadAuxBasis(nAt,rAt,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell) ! call ReadAuxBasis(nNuc,rNuc,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell)
! Compute the number of basis functions ! Compute the number of basis functions
@ -241,10 +240,10 @@ program MCQC
call cpu_time(start_MP2F12) call cpu_time(start_MP2F12)
! Memory allocation for one- and two-electron integrals ! Memory allocation for one- and two-electron integrals
! allocate(F12(nBas,nBas,nBas,nBas),Yuk(nBas,nBas,nBas,nBas)) allocate(F12(nBas,nBas,nBas,nBas),Yuk(nBas,nBas,nBas,nBas),FC(nBas,nBas,nBas,nBas,nBas,nBas))
! Read integrals ! Read integrals
! call ReadF12Ints(nBas,S,G,F12,Yuk) call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC)
! call MP2F12(nBas,nC,nO,nV,nA,G,F12,Yuk,ERHF,EcMP2(1),c,c,EcMP2F12) call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,EcMP2(1),cHF,cHF,EcMP2F12)
call cpu_time(end_MP2F12) call cpu_time(end_MP2F12)
t_MP2F12 = end_MP2F12 - start_MP2F12 t_MP2F12 = end_MP2F12 - start_MP2F12

View File

@ -1,24 +1,30 @@
subroutine MP2F12(nBas,nC,nO,nV,nA,ERI,F12,Yuk,EHF,EcMP2,c,cA,EcMP2F12) subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12)
! Perform restricted Hartree-Fock calculation ! Perform MP2-F12 calculation
implicit none implicit none
! Input variables ! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nA integer,intent(in) :: nBas,nC,nO,nV
double precision,intent(in) :: EHF,EcMP2 double precision,intent(in) :: EHF,EcMP2
double precision,intent(in) :: c(nBas,nBas),cA(nBas,nBas) double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),F12(nBas,nBas,nBas,nBas),Yuk(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: F12(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Yuk(nBas,nBas,nBas,nBas)
double precision,intent(in) :: FC(nBas,nBas,nBas,nBas,nBas,nBas)
! Local variables ! Local variables
double precision,allocatable :: ooCoo(:,:,:,:),ooFoo(:,:,:,:),ooYoo(:,:,:,:) double precision,allocatable :: ooCoo(:,:,:,:)
double precision,allocatable :: ooCoa(:,:,:,:),ooFoa(:,:,:,:) double precision,allocatable :: ooFoo(:,:,:,:)
double precision,allocatable :: ooCvv(:,:,:,:),ooFvv(:,:,:,:) double precision,allocatable :: ooYoo(:,:,:,:)
double precision,allocatable :: ooCvv(:,:,:,:)
double precision,allocatable :: ooFvv(:,:,:,:)
double precision,allocatable :: oooFCooo(:,:,:,:,:,:)
double precision,allocatable :: cO(:,:),cV(:,:) double precision,allocatable :: cO(:,:),cV(:,:)
double precision :: E2a,E2b,E3a,E3b,E4a,E4b,E4c,E4d double precision :: E2a,E2b,E3a,E3b,E4a,E4b,E4c,E4d
integer :: i,j,k,l,a,b,x integer :: i,j,k,l,a,b
! Output variables ! Output variables
@ -48,24 +54,20 @@ subroutine MP2F12(nBas,nC,nO,nV,nA,ERI,F12,Yuk,EHF,EcMP2,c,cA,EcMP2F12)
! Compute the three-electron part of the MP2-F12 energy ! Compute the three-electron part of the MP2-F12 energy
allocate(ooCoa(nO,nO,nO,nA),ooFoa(nO,nO,nO,nA)) allocate(oooFCooo(nO,nO,nO,nO,nO,nO))
call AOtoMO_oooa(nBas,nO,nA,cO,cA,ERI,ooCoa)
call AOtoMO_oooa(nBas,nO,nA,cO,cA,F12,ooFoa)
E3a = 0d0 E3a = 0d0
E3b = 0d0 E3b = 0d0
do i=1,nO do i=1,nO
do j=1,nO do j=1,nO
do k=1,nO do k=1,nO
do x=1,nA E3a = E3a + oooFCooo(i,j,k,k,j,i)
E3a = E3a + ooCoa(i,j,k,x)*ooFoa(j,i,k,x) E3b = E3b + oooFCooo(i,j,k,k,i,j)
E3b = E3b + ooCoa(i,j,k,x)*ooFoa(i,j,k,x)
enddo
enddo enddo
enddo enddo
enddo enddo
deallocate(ooCoa,ooFoa) deallocate(oooFCooo)
! Compute the four-electron part of the MP2-F12 energy ! Compute the four-electron part of the MP2-F12 energy

View File

@ -1,4 +1,4 @@
subroutine read_F12_integrals(nBas,S,C,F,Y) subroutine read_F12_integrals(nBas,S,C,F,Y,FC)
! Read one- and two-electron integrals from files ! Read one- and two-electron integrals from files
@ -12,12 +12,15 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
! Local variables ! Local variables
logical :: debug logical :: debug
integer :: mu,nu,la,si integer :: mu,nu,la,si,ka,ta
double precision :: ERI,F12,Yuk,ExpS double precision :: ERI,F12,Yuk,ExpS
! Output variables ! Output variables
double precision,intent(out) :: C(nBas,nBas,nBas,nBas),F(nBas,nBas,nBas,nBas),Y(nBas,nBas,nBas,nBas) double precision,intent(out) :: C(nBas,nBas,nBas,nBas)
double precision,intent(out) :: F(nBas,nBas,nBas,nBas)
double precision,intent(out) :: Y(nBas,nBas,nBas,nBas)
double precision,intent(out) :: FC(nBas,nBas,nBas,nBas,nBas,nBas)
debug = .false. debug = .false.
@ -26,8 +29,9 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
open(unit=21,file='int/ERI.dat') open(unit=21,file='int/ERI.dat')
open(unit=22,file='int/F12.dat') open(unit=22,file='int/F12.dat')
open(unit=23,file='int/Yuk.dat') open(unit=23,file='int/Yuk.dat')
open(unit=31,file='int/3eInt_Type1.dat')
! Read electron repulsion integrals ! Read 1/r12 integrals
C = 0d0 C = 0d0
do do
@ -51,7 +55,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
enddo enddo
21 close(unit=21) 21 close(unit=21)
! Read F12 integrals ! Read f12 integrals
F = 0d0 F = 0d0
do do
@ -74,7 +78,8 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
F(si,la,nu,mu) = F12 F(si,la,nu,mu) = F12
enddo enddo
22 close(unit=22) 22 close(unit=22)
! Read electron repulsion integrals
! Read f12/r12 integrals
Y = 0d0 Y = 0d0
do do
@ -98,6 +103,14 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
enddo enddo
23 close(unit=23) 23 close(unit=23)
! Read f13/r12 integrals
FC = 0d0
do
read(31,*,end=31) mu,nu,la,si,ka,ta,FC
enddo
31 close(unit=31)
! Print results ! Print results
if(debug) then if(debug) then
@ -137,15 +150,15 @@ subroutine read_F12_integrals(nBas,S,C,F,Y)
! Transform two-electron integrals ! Transform two-electron integrals
do mu=1,nBas ! do mu=1,nBas
do nu=1,nBas ! do nu=1,nBas
do la=1,nBas ! do la=1,nBas
do si=1,nBas ! do si=1,nBas
F(mu,nu,la,si) = (S(mu,la)*S(nu,si) - F(mu,nu,la,si))/ExpS ! F(mu,nu,la,si) = (S(mu,la)*S(nu,si) - F(mu,nu,la,si))/ExpS
Y(mu,nu,la,si) = (C(mu,nu,la,si) - Y(mu,nu,la,si))/ExpS ! Y(mu,nu,la,si) = (C(mu,nu,la,si) - Y(mu,nu,la,si))/ExpS
enddo ! enddo
enddo ! enddo
enddo ! enddo
enddo ! enddo
end subroutine read_F12_integrals end subroutine read_F12_integrals

View File

@ -1,5 +1,5 @@
subroutine read_methods(doHF,doMOM, & subroutine read_methods(doHF,doMOM, &
doMP2,doMP3, & doMP2,doMP3,doMP2F12, &
doCIS,doTDHF,doADC, & doCIS,doTDHF,doADC, &
doGF2,doGF3, & doGF2,doGF3, &
doG0W0,doevGW,doqsGW, & doG0W0,doevGW,doqsGW, &
@ -12,7 +12,7 @@ subroutine read_methods(doHF,doMOM, &
! Input variables ! Input variables
logical,intent(out) :: doHF,doMOM logical,intent(out) :: doHF,doMOM
logical,intent(out) :: doMP2,doMP3 logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCIS,doTDHF,doADC logical,intent(out) :: doCIS,doTDHF,doADC
logical,intent(out) :: doGF2,doGF3 logical,intent(out) :: doGF2,doGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0W0,doevGW,doqsGW
@ -31,8 +31,9 @@ subroutine read_methods(doHF,doMOM, &
doHF = .false. doHF = .false.
doMOM = .false. doMOM = .false.
doMP2 = .false. doMP2 = .false.
doMP3 = .false. doMP3 = .false.
doMP2F12 = .false.
doCIS = .false. doCIS = .false.
doTDHF = .false. doTDHF = .false.
@ -57,9 +58,10 @@ subroutine read_methods(doHF,doMOM, &
! Read MPn methods ! Read MPn methods
read(1,*) read(1,*)
read(1,*) answer1,answer2 read(1,*) answer1,answer2,answer3
if(answer1 == 'T') doMP2 = .true. if(answer1 == 'T') doMP2 = .true.
if(answer2 == 'T') doMP3 = .true. if(answer2 == 'T') doMP3 = .true.
if(answer3 == 'T') doMP2F12 = .true.
! Read excited state methods ! Read excited state methods