From b41fbba5152ade55c25c18afb067b70217bbec3b Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 16 Mar 2019 14:09:42 +0100 Subject: [PATCH] CCD is working --- input/basis | 29 ++++++++++++++++++++--------- input/int | 2 +- input/molecule | 4 ++-- input/weight | 29 ++++++++++++++++++++--------- src/MCQC/MCQC.f90 | 18 ++---------------- src/MCQC/read_geometry.f90 | 20 ++++++++++---------- src/MCQC/read_integrals.f90 | 33 +++++++++++++++++++-------------- src/MCQC/read_molecule.f90 | 10 +++++----- 8 files changed, 79 insertions(+), 66 deletions(-) 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..ed14d90 100644 --- a/input/int +++ b/input/int @@ -1,7 +1,7 @@ # Debuggin mode? F # Chemist notation for two-electron integral? - T + F # Exposant of the Slater geminal 1.0 # One-electron integrals: Ov Kin Nuc diff --git a/input/molecule b/input/molecule index 7d017f4..67ee9ac 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEl nCore nRyd - 1 2 0 0 + 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/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/MCQC/MCQC.f90 b/src/MCQC/MCQC.f90 index 9d70d88..5b45c63 100644 --- a/src/MCQC/MCQC.f90 +++ b/src/MCQC/MCQC.f90 @@ -14,7 +14,6 @@ program MCQC 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 :: EcCCD,EcCCSD,EcCCSDT double precision,allocatable :: ZNuc(:),rNuc(:,:),cHF(:,:),eHF(:),eG0W0(:),PHF(:,:) @@ -33,7 +32,6 @@ program MCQC double precision :: start_MOM ,end_MOM ,t_MOM double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CCSD ,end_CCSD ,t_CCSD - double precision :: start_CCSDT ,end_CCSDT ,t_CCSDT double precision :: start_CIS ,end_CIS ,t_CIS double precision :: start_TDHF ,end_TDHF ,t_TDHF double precision :: start_ADC ,end_ADC ,t_ADC @@ -202,7 +200,6 @@ program MCQC ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ - call chem_to_phys_ERI(nBas,ERI_AO_basis) call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) !------------------------------------------------------------------------ @@ -265,7 +262,7 @@ program MCQC if(doCCD) then call cpu_time(start_CCD) - call CCD(nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF,cHF,EcCCD) + call CCD(nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF,cHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -281,20 +278,9 @@ program MCQC if(doCCSD) then call cpu_time(start_CCSD) - call CCSD(doCCSDT,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF,cHF,EcCCSD) + call CCSD(doCCSDT,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF,cHF) call cpu_time(end_CCSD) -! if(doCCSDT) then - call cpu_time(start_CCSDT) -! call CCSDT(nBas,nEl,ERI_MO_basis,ENuc,ERHF,EcCCSD,eHF,cHF,EcCCSDT) - call cpu_time(end_CCSDT) - -! t_CCSDT = end_CCSDT - start_CCSDT -! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds' -! write(*,*) - -! end if - t_CCSD = end_CCSD - start_CCSD write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' write(*,*) diff --git a/src/MCQC/read_geometry.f90 b/src/MCQC/read_geometry.f90 index 60c60b8..a3a4cf9 100644 --- a/src/MCQC/read_geometry.f90 +++ b/src/MCQC/read_geometry.f90 @@ -1,4 +1,4 @@ -subroutine read_geometry(nAt,ZNuc,rA,ENuc) +subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) ! Read molecular geometry @@ -8,7 +8,7 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc) ! Ouput variables - integer,intent(in) :: nAt + integer,intent(in) :: nNuc ! Local variables @@ -19,7 +19,7 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc) ! Ouput variables - double precision,intent(out) :: ZNuc(NAt),rA(nAt,ncart),ENuc + double precision,intent(out) :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc ! Open file with geometry specification @@ -31,8 +31,8 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc) read(1,*) read(1,*) - do i=1,nAt - read(1,*) El,rA(i,1),rA(i,2),rA(i,3) + do i=1,nNuc + read(1,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3) ZNuc(i) = element_number(El) enddo @@ -40,9 +40,9 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc) ENuc = 0 - do i=1,nAt-1 - do j=i+1,nAt - RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2 + do i=1,nNuc-1 + do j=i+1,nNuc + RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2 ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB) enddo enddo @@ -54,10 +54,10 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc) write(*,'(A28)') '------------------' write(*,'(A28)') 'Molecular geometry' write(*,'(A28)') '------------------' - do i=1,NAt + do i=1,nNuc write(*,'(A28,1X,I16)') 'Atom n. ',i write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i) - write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,ncart) + write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rNuc(i,j),j=1,ncart) enddo write(*,*) write(*,'(A28)') '------------------' diff --git a/src/MCQC/read_integrals.f90 b/src/MCQC/read_integrals.f90 index fa435d1..e6f9dd6 100644 --- a/src/MCQC/read_integrals.f90 +++ b/src/MCQC/read_integrals.f90 @@ -13,6 +13,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) logical :: debug integer :: mu,nu,la,si double precision :: Ov,Kin,Nuc,ERI + double precision :: scale ! Output variables @@ -22,6 +23,8 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) debug = .false. + scale = 1d0 + open(unit=8 ,file='int/Ov.dat') open(unit=9 ,file='int/Kin.dat') open(unit=10,file='int/Nuc.dat') @@ -41,7 +44,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) T = 0d0 do read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Kin + T(mu,nu) = Kin/scale**2 enddo 9 close(unit=9) @@ -63,21 +66,23 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) G = 0d0 do read(11,*,end=11) mu,nu,la,si,ERI -! (12|34) + + ERI = ERI/scale +! <12|34> G(mu,nu,la,si) = ERI -! (21|34) - G(nu,mu,la,si) = ERI -! (12|43) - G(mu,nu,si,la) = ERI -! (21|43) - G(nu,mu,si,la) = ERI -! (34|12) +! <32|14> + G(la,nu,mu,si) = ERI +! <14|32> + G(mu,si,la,nu) = ERI +! <34|12> G(la,si,mu,nu) = ERI -! (43|12) - G(si,la,mu,nu) = ERI -! (34|21) - G(la,si,nu,mu) = ERI -! (43|21) +! <41|23> + G(si,mu,nu,la) = ERI +! <23|41> + G(nu,la,si,mu) = ERI +! <21|43> + G(nu,mu,si,la) = ERI +! <43|21> G(si,la,nu,mu) = ERI enddo 11 close(unit=11) diff --git a/src/MCQC/read_molecule.f90 b/src/MCQC/read_molecule.f90 index 0894e17..4625d03 100644 --- a/src/MCQC/read_molecule.f90 +++ b/src/MCQC/read_molecule.f90 @@ -1,6 +1,6 @@ -subroutine read_molecule(nAt,nEl,nO,nCore,nRyd) +subroutine read_molecule(nNuc,nEl,nO,nC,nR) -! Read number of atoms nAt and number of electrons nEl +! Read number of atoms and number of electrons implicit none @@ -8,7 +8,7 @@ subroutine read_molecule(nAt,nEl,nO,nCore,nRyd) ! Input variables - integer,intent(out) :: nAt,nEl,nO,nCore,nRyd + integer,intent(out) :: nNuc,nEl,nO,nC,nR ! Open file with geometry specification @@ -17,14 +17,14 @@ subroutine read_molecule(nAt,nEl,nO,nCore,nRyd) ! Read number of atoms and number of electrons read(1,*) - read(1,*) nAt,nEl,nCore,nRyd + read(1,*) nNuc,nEl,nC,nR nO = nEl/2 ! Print results write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of atoms',nAt + write(*,'(A28,1X,I16)') 'Number of atoms',nNuc write(*,'(A28)') '----------------------' write(*,*) write(*,'(A28)') '----------------------'