diff --git a/input/basis b/input/basis index b9ca7b5..c12e4d6 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,20 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +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 diff --git a/input/int b/input/int index 86aa017..6f9394e 100644 --- a/input/int +++ b/input/int @@ -5,10 +5,10 @@ # Exposant of the Slater geminal 1.0 # One-electron integrals: Ov Kin Nuc - T T T + F F F # Two-electron integrals: ERI F12 Yuk Erf T F F F # Three-electron integrals: Type1 Type2 Type3 - F F F + T F F # Four-electron integrals: Type1 Type2 Type3 F F F diff --git a/input/methods b/input/methods index 183d9ac..efe6a18 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # HF MOM T F -# MP2 MP3 - F F +# MP2 MP3 MP2-F12 + F F T # CIS TDHF ADC F F F # GF2 GF3 diff --git a/input/molecule b/input/molecule index 375e1fb..67ee9ac 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ -# nAt nEl nEla nElb - 1 2 1 1 +# nAt nEl nCore nRyd + 1 4 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Be 0.0 0.0 0.0 diff --git a/input/options b/input/options index 5549cb1..6d94128 100644 --- a/input/options +++ b/input/options @@ -1,7 +1,7 @@ # RHF: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.0000001 T 5 1 1 -# MPn: - +# MP: MP2-F12 + T # CIS/TDHF: singlet triplet T F # GF: maxSCF thresh DIIS n_diis renormalization diff --git a/input/weight b/input/weight index b9ca7b5..c12e4d6 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,20 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +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 diff --git a/src/IntPak/ComputeOv.f90 b/src/IntPak/ComputeOv.f90 index 3aceb65..9c8c5d8 100644 --- a/src/IntPak/ComputeOv.f90 +++ b/src/IntPak/ComputeOv.f90 @@ -1,4 +1,4 @@ -subroutine ComputeOv(debug,NBasis,nShell, & +subroutine ComputeOv(debug,nBas,nShell, & CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & npOv,nSigpOv,ncOv,nSigcOv,S) @@ -11,7 +11,7 @@ subroutine ComputeOv(debug,NBasis,nShell, & ! Input variables logical,intent(in) :: debug - integer,intent(in) :: NBasis,nShell + integer,intent(in) :: nBas,nShell double precision,intent(in) :: CenterShell(maxShell,3) integer,intent(in) :: TotAngMomShell(maxShell),KShell(maxShell) double precision,intent(in) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) @@ -39,7 +39,7 @@ subroutine ComputeOv(debug,NBasis,nShell, & ! Output variables 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 diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90 index 1f51110..fb23ec2 100644 --- a/src/IntPak/IntPak.f90 +++ b/src/IntPak/IntPak.f90 @@ -19,7 +19,7 @@ program IntPak logical :: do4eInt(n4eInt) integer :: nNuc,nBas,iType - integer :: nEl,nO,nV + integer :: nEl,nO,nV,nCore,nRyd double precision :: ExpS double precision :: ENuc integer :: KG @@ -69,7 +69,7 @@ program IntPak ! nBas = number of basis functions (see below) ! = 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)) @@ -346,13 +346,13 @@ program IntPak if(do3eInt(1)) then iType = 1 - KG = 1 -! KG = 6 +! KG = 1 + KG = 6 allocate(DG(1:KG),ExpG(1:KG)) - DG = (/ 1d0 /) - ExpG = (/ 1d0 /) -! 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 /) +! DG = (/ 1d0 /) +! ExpG = (/ 1d0 /) + 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 /) call cpu_time(start_3eInt(iType)) call Compute3eInt(debug,iType,nShell, & diff --git a/src/IntPak/S3eInt.f90 b/src/IntPak/S3eInt.f90 index faccdaf..4c45eb3 100644 --- a/src/IntPak/S3eInt.f90 +++ b/src/IntPak/S3eInt.f90 @@ -34,17 +34,35 @@ subroutine S3eInt(debug,iType,np3eInt,nSigp3eInt, & double precision :: p3eInt p3eInt = 0d0 - do k=1,KG - do l=1,KG + + if(iType == 1) then + + do k=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) - enddo - enddo + p3eInt = p3eInt & + + DG(k)*G3eInt(debug,iType, & + ExpSG13,ExpSG23, & + ExpBra,CenterBra,AngMomBra, & + ExpKet,CenterKet,AngMomKet) + end do + + 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 diff --git a/src/MCQC/MCQC.f90 b/src/MCQC/MCQC.f90 index 42c8f3a..fb9315f 100644 --- a/src/MCQC/MCQC.f90 +++ b/src/MCQC/MCQC.f90 @@ -10,11 +10,11 @@ program MCQC logical :: doG0W0,doevGW,doqsGW logical :: doMCMP2,doMinMCMP2 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 :: 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,allocatable :: TotAngMomShell(:),KShell(:) @@ -25,7 +25,7 @@ program MCQC double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:) 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_MOM ,end_MOM ,t_MOM @@ -77,7 +77,7 @@ program MCQC ! Which calculations do you want to do? call read_methods(doHF,doMOM, & - doMP2,doMP3, & + doMP2,doMP3,doMP2F12, & doCIS,doTDHF,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & @@ -95,7 +95,6 @@ program MCQC doeNCusp = .false. doMinMCMP2 = .false. - doMP2F12 = .false. !------------------------------------------------------------------------ ! Read input information @@ -109,12 +108,12 @@ program MCQC ! nBas = number of basis functions (see below) ! = nO + nV - call read_molecule(nAt,nEl,nC,nO,nR) - allocate(ZNuc(nAt),rAt(nAt,3)) + call read_molecule(nNuc,nEl,nC,nO,nR) + allocate(ZNuc(nNuc),rNuc(nNuc,3)) ! Read geometry - call read_geometry(nAt,ZNuc,rAt,ENuc) + call read_geometry(nNuc,ZNuc,rNuc,ENuc) allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), & DShell(maxShell,maxK),ExpShell(maxShell,maxK)) @@ -123,14 +122,14 @@ program MCQC ! 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) !------------------------------------------------------------------------ ! 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 @@ -241,10 +240,10 @@ program MCQC call cpu_time(start_MP2F12) ! 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 -! call ReadF12Ints(nBas,S,G,F12,Yuk) -! call MP2F12(nBas,nC,nO,nV,nA,G,F12,Yuk,ERHF,EcMP2(1),c,c,EcMP2F12) + call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC) + call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,EcMP2(1),cHF,cHF,EcMP2F12) call cpu_time(end_MP2F12) t_MP2F12 = end_MP2F12 - start_MP2F12 diff --git a/src/MCQC/MP2F12.f90 b/src/MCQC/MP2F12.f90 index 27cf274..5fda128 100644 --- a/src/MCQC/MP2F12.f90 +++ b/src/MCQC/MP2F12.f90 @@ -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 ! 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) :: c(nBas,nBas),cA(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) :: c(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 - double precision,allocatable :: ooCoo(:,:,:,:),ooFoo(:,:,:,:),ooYoo(:,:,:,:) - double precision,allocatable :: ooCoa(:,:,:,:),ooFoa(:,:,:,:) - double precision,allocatable :: ooCvv(:,:,:,:),ooFvv(:,:,:,:) + double precision,allocatable :: ooCoo(:,:,:,:) + double precision,allocatable :: ooFoo(:,:,:,:) + double precision,allocatable :: ooYoo(:,:,:,:) + double precision,allocatable :: ooCvv(:,:,:,:) + double precision,allocatable :: ooFvv(:,:,:,:) + double precision,allocatable :: oooFCooo(:,:,:,:,:,:) double precision,allocatable :: cO(:,:),cV(:,:) 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 @@ -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 - allocate(ooCoa(nO,nO,nO,nA),ooFoa(nO,nO,nO,nA)) - call AOtoMO_oooa(nBas,nO,nA,cO,cA,ERI,ooCoa) - call AOtoMO_oooa(nBas,nO,nA,cO,cA,F12,ooFoa) + allocate(oooFCooo(nO,nO,nO,nO,nO,nO)) E3a = 0d0 E3b = 0d0 do i=1,nO do j=1,nO do k=1,nO - do x=1,nA - E3a = E3a + ooCoa(i,j,k,x)*ooFoa(j,i,k,x) - E3b = E3b + ooCoa(i,j,k,x)*ooFoa(i,j,k,x) - enddo + E3a = E3a + oooFCooo(i,j,k,k,j,i) + E3b = E3b + oooFCooo(i,j,k,k,i,j) enddo enddo enddo - deallocate(ooCoa,ooFoa) + deallocate(oooFCooo) ! Compute the four-electron part of the MP2-F12 energy diff --git a/src/MCQC/read_F12_integrals.f90 b/src/MCQC/read_F12_integrals.f90 index 3ef6461..0846be2 100644 --- a/src/MCQC/read_F12_integrals.f90 +++ b/src/MCQC/read_F12_integrals.f90 @@ -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 @@ -12,12 +12,15 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) ! Local variables logical :: debug - integer :: mu,nu,la,si + integer :: mu,nu,la,si,ka,ta double precision :: ERI,F12,Yuk,ExpS ! 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. @@ -26,8 +29,9 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) open(unit=21,file='int/ERI.dat') open(unit=22,file='int/F12.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 do @@ -51,7 +55,7 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) enddo 21 close(unit=21) -! Read F12 integrals +! Read f12 integrals F = 0d0 do @@ -74,7 +78,8 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) F(si,la,nu,mu) = F12 enddo 22 close(unit=22) -! Read electron repulsion integrals + +! Read f12/r12 integrals Y = 0d0 do @@ -98,6 +103,14 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) enddo 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 if(debug) then @@ -137,15 +150,15 @@ subroutine read_F12_integrals(nBas,S,C,F,Y) ! Transform two-electron integrals - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - 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 - enddo - enddo - enddo - enddo +! do mu=1,nBas +! do nu=1,nBas +! do la=1,nBas +! do si=1,nBas +! 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 +! enddo +! enddo +! enddo +! enddo end subroutine read_F12_integrals diff --git a/src/MCQC/read_methods.f90 b/src/MCQC/read_methods.f90 index 6db5635..f4389cb 100644 --- a/src/MCQC/read_methods.f90 +++ b/src/MCQC/read_methods.f90 @@ -1,5 +1,5 @@ subroutine read_methods(doHF,doMOM, & - doMP2,doMP3, & + doMP2,doMP3,doMP2F12, & doCIS,doTDHF,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & @@ -12,7 +12,7 @@ subroutine read_methods(doHF,doMOM, & ! Input variables 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) :: doGF2,doGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW @@ -31,8 +31,9 @@ subroutine read_methods(doHF,doMOM, & doHF = .false. doMOM = .false. - doMP2 = .false. - doMP3 = .false. + doMP2 = .false. + doMP3 = .false. + doMP2F12 = .false. doCIS = .false. doTDHF = .false. @@ -57,9 +58,10 @@ subroutine read_methods(doHF,doMOM, & ! Read MPn methods read(1,*) - read(1,*) answer1,answer2 - if(answer1 == 'T') doMP2 = .true. - if(answer2 == 'T') doMP3 = .true. + read(1,*) answer1,answer2,answer3 + if(answer1 == 'T') doMP2 = .true. + if(answer2 == 'T') doMP3 = .true. + if(answer3 == 'T') doMP2F12 = .true. ! Read excited state methods