10
1
mirror of https://github.com/pfloos/quack synced 2024-07-04 02:16:06 +02:00

CCD is working

This commit is contained in:
Pierre-Francois Loos 2019-03-16 14:09:42 +01:00
parent 99acd24d99
commit b41fbba515
8 changed files with 79 additions and 66 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

@ -1,7 +1,7 @@
# Debuggin mode? # Debuggin mode?
F F
# Chemist notation for two-electron integral? # Chemist notation for two-electron integral?
T F
# 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

View File

@ -1,4 +1,4 @@
# nAt nEl nCore nRyd # nAt nEl nCore nRyd
1 2 0 0 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,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

@ -14,7 +14,6 @@ program MCQC
integer :: nNuc,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 :: EcCCD,EcCCSD,EcCCSDT
double precision,allocatable :: ZNuc(:),rNuc(:,:),cHF(:,:),eHF(:),eG0W0(:),PHF(:,:) 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_MOM ,end_MOM ,t_MOM
double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CCD ,end_CCD ,t_CCD
double precision :: start_CCSD ,end_CCSD ,t_CCSD 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_CIS ,end_CIS ,t_CIS
double precision :: start_TDHF ,end_TDHF ,t_TDHF double precision :: start_TDHF ,end_TDHF ,t_TDHF
double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_ADC ,end_ADC ,t_ADC
@ -202,7 +200,6 @@ program MCQC
! AO to MO integral transform for post-HF methods ! 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) call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -265,7 +262,7 @@ program MCQC
if(doCCD) then if(doCCD) then
call cpu_time(start_CCD) 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) call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD t_CCD = end_CCD - start_CCD
@ -281,20 +278,9 @@ program MCQC
if(doCCSD) then if(doCCSD) then
call cpu_time(start_CCSD) 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) 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 t_CCSD = end_CCSD - start_CCSD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds'
write(*,*) write(*,*)

View File

@ -1,4 +1,4 @@
subroutine read_geometry(nAt,ZNuc,rA,ENuc) subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
! Read molecular geometry ! Read molecular geometry
@ -8,7 +8,7 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc)
! Ouput variables ! Ouput variables
integer,intent(in) :: nAt integer,intent(in) :: nNuc
! Local variables ! Local variables
@ -19,7 +19,7 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc)
! Ouput variables ! 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 ! Open file with geometry specification
@ -31,8 +31,8 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc)
read(1,*) read(1,*)
read(1,*) read(1,*)
do i=1,nAt do i=1,nNuc
read(1,*) El,rA(i,1),rA(i,2),rA(i,3) read(1,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3)
ZNuc(i) = element_number(El) ZNuc(i) = element_number(El)
enddo enddo
@ -40,9 +40,9 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc)
ENuc = 0 ENuc = 0
do i=1,nAt-1 do i=1,nNuc-1
do j=i+1,nAt do j=i+1,nNuc
RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2 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) ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
enddo enddo
enddo enddo
@ -54,10 +54,10 @@ subroutine read_geometry(nAt,ZNuc,rA,ENuc)
write(*,'(A28)') '------------------' write(*,'(A28)') '------------------'
write(*,'(A28)') 'Molecular geometry' write(*,'(A28)') 'Molecular geometry'
write(*,'(A28)') '------------------' write(*,'(A28)') '------------------'
do i=1,NAt do i=1,nNuc
write(*,'(A28,1X,I16)') 'Atom n. ',i write(*,'(A28,1X,I16)') 'Atom n. ',i
write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(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 enddo
write(*,*) write(*,*)
write(*,'(A28)') '------------------' write(*,'(A28)') '------------------'

View File

@ -13,6 +13,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
logical :: debug logical :: debug
integer :: mu,nu,la,si integer :: mu,nu,la,si
double precision :: Ov,Kin,Nuc,ERI double precision :: Ov,Kin,Nuc,ERI
double precision :: scale
! Output variables ! Output variables
@ -22,6 +23,8 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
debug = .false. debug = .false.
scale = 1d0
open(unit=8 ,file='int/Ov.dat') open(unit=8 ,file='int/Ov.dat')
open(unit=9 ,file='int/Kin.dat') open(unit=9 ,file='int/Kin.dat')
open(unit=10,file='int/Nuc.dat') open(unit=10,file='int/Nuc.dat')
@ -41,7 +44,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
T = 0d0 T = 0d0
do do
read(9,*,end=9) mu,nu,Kin read(9,*,end=9) mu,nu,Kin
T(mu,nu) = Kin T(mu,nu) = Kin/scale**2
enddo enddo
9 close(unit=9) 9 close(unit=9)
@ -63,21 +66,23 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
G = 0d0 G = 0d0
do do
read(11,*,end=11) mu,nu,la,si,ERI read(11,*,end=11) mu,nu,la,si,ERI
! (12|34)
ERI = ERI/scale
! <12|34>
G(mu,nu,la,si) = ERI G(mu,nu,la,si) = ERI
! (21|34) ! <32|14>
G(nu,mu,la,si) = ERI G(la,nu,mu,si) = ERI
! (12|43) ! <14|32>
G(mu,nu,si,la) = ERI G(mu,si,la,nu) = ERI
! (21|43) ! <34|12>
G(nu,mu,si,la) = ERI
! (34|12)
G(la,si,mu,nu) = ERI G(la,si,mu,nu) = ERI
! (43|12) ! <41|23>
G(si,la,mu,nu) = ERI G(si,mu,nu,la) = ERI
! (34|21) ! <23|41>
G(la,si,nu,mu) = ERI G(nu,la,si,mu) = ERI
! (43|21) ! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI G(si,la,nu,mu) = ERI
enddo enddo
11 close(unit=11) 11 close(unit=11)

View File

@ -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 implicit none
@ -8,7 +8,7 @@ subroutine read_molecule(nAt,nEl,nO,nCore,nRyd)
! Input variables ! Input variables
integer,intent(out) :: nAt,nEl,nO,nCore,nRyd integer,intent(out) :: nNuc,nEl,nO,nC,nR
! Open file with geometry specification ! 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 number of atoms and number of electrons
read(1,*) read(1,*)
read(1,*) nAt,nEl,nCore,nRyd read(1,*) nNuc,nEl,nC,nR
nO = nEl/2 nO = nEl/2
! Print results ! Print results
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of atoms',nAt write(*,'(A28,1X,I16)') 'Number of atoms',nNuc
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'