LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00
Code Releases Activity

crossed ring CC and RPA working

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-11-10 14:47:26 +01:00
parent cec0fe9205
commit 3dfc2c3526
11 changed files with 197 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
input/methods
View File

@ -2,14 +2,14 @@
T F F F T F F F
# MP2* MP3 MP2-F12 # MP2* MP3 MP2-F12
F F F F F F
# CCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
F F F F F F F F F
# drCCD rCCD lCCD pCCD # drCCD rCCD crCCD lCCD
T T F F T T T T
# CIS* CIS(D) CID CISD FCI # CIS* CIS(D) CID CISD FCI
F F F F F F F F F F
# RPA* RPAx* ppRPA # RPA* RPAx* crRPA ppRPA
T T F T T T T
# G0F2* evGF2* qsGF2* G0F3 evGF3 # G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW # G0W0* evGW* qsGW* ufG0W0 ufGW

54
src/CC/crCCD.f90
View File

@ -38,21 +38,8 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
double precision,allocatable :: eV(:) double precision,allocatable :: eV(:)
double precision,allocatable :: delta_OOVV(:,:,:,:) double precision,allocatable :: delta_OOVV(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVOV(:,:,:,:) double precision,allocatable :: OVOV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VVVV(:,:,:,:)
double precision,allocatable :: X1(:,:,:,:)
double precision,allocatable :: X2(:,:)
double precision,allocatable :: X3(:,:)
double precision,allocatable :: X4(:,:,:,:)
double precision,allocatable :: u(:,:,:,:)
double precision,allocatable :: v(:,:,:,:)
double precision,allocatable :: r2r(:,:,:,:)
double precision,allocatable :: r2l(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:)
@ -66,7 +53,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*) write(*,*)
write(*,*)'**************************************' write(*,*)'**************************************'
write(*,*)'| crossed-ring CCD calculation |' write(*,*)'| Crossed-ring CCD calculation |'
write(*,*)'**************************************' write(*,*)'**************************************'
write(*,*) write(*,*)
@ -105,13 +92,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Create integral batches ! Create integral batches
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV),OVVO(nO,nV,nV,nO)) allocate(OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV))
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas)
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
deallocate(dbERI) deallocate(dbERI)
@ -129,9 +113,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Initialization ! Initialization
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) allocate(r2(nO,nO,nV,nV))
allocate(r2r(nO,nO,nV,nV),r2l(nO,nO,nV,nV))
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
Conv = 1d0 Conv = 1d0
nSCF = 0 nSCF = 0
@ -159,23 +141,9 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Compute residual ! Compute residual
! Form linear array call form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
! Form interemediate arrays
call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
! Form quadratic array
call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2r)
call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2l)
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) - r2r(:,:,:,:) - r2l(:,:,:,:)
! Check convergence ! Check convergence
@ -183,7 +151,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Update amplitudes ! Update amplitudes
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) t2(:,:,:,:) = t2(:,:,:,:) + r2(:,:,:,:)/delta_OOVV(:,:,:,:)
! Compute correlation energy ! Compute correlation energy
@ -227,10 +195,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*)' crossed-ring CCD energy ' write(*,*)' crossed-ring CCD energy '
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD
write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)

6
src/CC/drCCD.f90
View File

@ -39,7 +39,6 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VOOV(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:)
@ -85,11 +84,10 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Create integral batches ! Create integral batches
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
VOOV(:,:,:,:) = sERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas)
deallocate(sERI) deallocate(sERI)
@ -135,7 +133,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
! Compute residual ! Compute residual
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)

70
src/CC/form_crossed_ring_r.f90 Normal file
View File

@ -0,0 +1,70 @@
subroutine form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
! Form residuals for crossed-ring CCD
implicit none
! Input variables
integer,intent(in) :: nC,nO,nV,nR
double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: OVOV(nO,nV,nO,nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
double precision,allocatable :: y(:,:,:,:)
! Output variables
double precision,intent(out) :: r2(nO,nO,nV,nV)
r2(:,:,:,:) = 0d0
allocate(y(nV,nO,nO,nV))
y(:,:,:,:) = 0d0
do i=nC+1,nO
do b=1,nV-nR
do l=nC+1,nO
do d=1,nV-nR
do k=nC+1,nO
do c=1,nV-nR
y(b,i,l,d) = y(b,i,l,d) + t2(i,k,c,b)*OOVV(k,l,c,d)
end do
end do
end do
end do
end do
end do
do i=nC+1,nO
do j=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
do k=nC+1,nO
do c=1,nV-nR
r2(i,j,a,b) = r2(i,j,a,b) - OVOV(k,b,i,c)*t2(k,j,a,c) - OVOV(k,a,j,c)*t2(i,k,c,b)
end do
end do
do l=nC+1,nO
do d=1,nV-nR
r2(i,j,a,b) = r2(i,j,a,b) - y(b,i,l,d)*t2(l,j,a,d)
end do
end do
end do
end do
end do
end do
end subroutine form_crossed_ring_r

5
src/CC/form_ring_r.f90
View File

@ -1,4 +1,4 @@
subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
! Form residuals for ring CCD ! Form residuals for ring CCD
@ -9,7 +9,6 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
integer,intent(in) :: nC,nO,nV,nR integer,intent(in) :: nC,nO,nV,nR
double precision,intent(in) :: t2(nO,nO,nV,nV) double precision,intent(in) :: t2(nO,nO,nV,nV)
double precision,intent(in) :: OVVO(nO,nV,nV,nO) double precision,intent(in) :: OVVO(nO,nV,nV,nO)
double precision,intent(in) :: VOOV(nV,nO,nO,nV)
double precision,intent(in) :: OOVV(nO,nO,nV,nV) double precision,intent(in) :: OOVV(nO,nO,nV,nV)
! Local variables ! Local variables
@ -50,7 +49,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2)
do k=nC+1,nO do k=nC+1,nO
do c=1,nV-nR do c=1,nV-nR
r2(i,j,a,b) = r2(i,j,a,b) + VOOV(a,k,i,c)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c) r2(i,j,a,b) = r2(i,j,a,b) + OVVO(k,a,c,i)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c)
end do end do
end do end do

6
src/CC/rCCD.f90
View File

@ -40,7 +40,6 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VOOV(:,:,:,:)
double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: r2(:,:,:,:)
double precision,allocatable :: t2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:)
@ -93,11 +92,10 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
! Create integral batches ! Create integral batches
allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO))
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
VOOV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas)
deallocate(dbERI) deallocate(dbERI)
@ -143,7 +141,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
! Compute residual ! Compute residual
call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)

103
src/QuAcK/QuAcK.f90
View File

@ -9,10 +9,10 @@ program QuAcK
logical :: dostab logical :: dostab
logical :: doKS logical :: doKS
logical :: doMP2,doMP3,doMP2F12 logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doDCD,doCCSD,doCCSDT logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
logical :: doCIS,doCIS_D,doCID,doCISD,doFCI logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical :: doRPA,doRPAx,doppRPA logical :: doRPA,doRPAx,docrRPA,doppRPA
logical :: doADC logical :: doADC
logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
@ -91,8 +91,6 @@ program QuAcK
double precision :: start_CISD ,end_CISD ,t_CISD double precision :: start_CISD ,end_CISD ,t_CISD
double precision :: start_FCI ,end_FCI ,t_FCI double precision :: start_FCI ,end_FCI ,t_FCI
double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_RPAx ,end_RPAx ,t_RPAx
double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA
double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_ADC ,end_ADC ,t_ADC
double precision :: start_GF2 ,end_GF2 ,t_GF2 double precision :: start_GF2 ,end_GF2 ,t_GF2
double precision :: start_GF3 ,end_GF3 ,t_GF3 double precision :: start_GF3 ,end_GF3 ,t_GF3
@ -160,17 +158,17 @@ program QuAcK
! Which calculations do you want to do? ! Which calculations do you want to do?
call read_methods(doRHF,doUHF,doKS,doMOM, & call read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, & doRPA,doRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, & doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, & doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, & doG0T0,doevGT,doqsGT, &
doMCMP2) doMCMP2)
! Read options for methods ! Read options for methods
@ -645,6 +643,24 @@ program QuAcK
end if end if
!------------------------------------------------------------------------
! Perform crossed-ring CCD calculation
!------------------------------------------------------------------------
if(do_crCCD) then
call cpu_time(start_CCD)
call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Perform ladder CCD calculation ! Perform ladder CCD calculation
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -662,32 +678,11 @@ program QuAcK
end if end if
!------------------------------------------------------------------------
! Perform crossed-ring CCD calculation
!------------------------------------------------------------------------
do_crCCD = .false.
if(do_crCCD) then
call cpu_time(start_CCD)
call ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Perform pair CCD calculation ! Perform pair CCD calculation
!------------------------------------------------------------------------ !------------------------------------------------------------------------
if(do_pCCD) then if(dopCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
@ -787,7 +782,7 @@ program QuAcK
if(doRPAx) then if(doRPAx) then
call cpu_time(start_RPAx) call cpu_time(start_RPA)
if(unrestricted) then if(unrestricted) then
call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
@ -798,10 +793,26 @@ program QuAcK
call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
end if end if
call cpu_time(end_RPAx) call cpu_time(end_RPA)
t_RPAx = end_RPAx - start_RPAx t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute cr-RPA excitations
!------------------------------------------------------------------------
if(docrRPA) then
call cpu_time(start_RPA)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if
@ -812,12 +823,12 @@ program QuAcK
if(doppRPA) then if(doppRPA) then
call cpu_time(start_ppRPA) call cpu_time(start_RPA)
call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA) call cpu_time(end_RPA)
t_ppRPA = end_ppRPA - start_ppRPA t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_ppRPA,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*) write(*,*)
end if end if

52
src/QuAcK/read_methods.f90
View File

@ -1,14 +1,14 @@
subroutine read_methods(doRHF,doUHF,doKS,doMOM, & subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, & doRPA,doRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, & doG0W0,doevGW,doqsGW, &
doufG0W0,doufGW, & doufG0W0,doufGW, &
doG0T0,doevGT,doqsGT, & doG0T0,doevGT,doqsGT, &
doMCMP2) doMCMP2)
! Read desired methods ! Read desired methods
@ -19,10 +19,10 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
logical,intent(out) :: doRHF,doUHF,doKS,doMOM logical,intent(out) :: doRHF,doUHF,doKS,doMOM
logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(out) :: doRPA,doRPAx,doppRPA logical,intent(out) :: doRPA,doRPAx,docrRPA,doppRPA
logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(out) :: doG0T0,doevGT,doqsGT logical,intent(out) :: doG0T0,doevGT,doqsGT
@ -48,14 +48,15 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2F12 = .false. doMP2F12 = .false.
doCCD = .false. doCCD = .false.
dopCCD = .false.
doDCD = .false. doDCD = .false.
doCCSD = .false. doCCSD = .false.
doCCSDT = .false. doCCSDT = .false.
do_drCCD = .false. do_drCCD = .false.
do_rCCD = .false. do_rCCD = .false.
do_crCCD = .false.
do_lCCD = .false. do_lCCD = .false.
do_pCCD = .false.
doCIS = .false. doCIS = .false.
doCIS_D = .false. doCIS_D = .false.
@ -65,6 +66,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doRPA = .false. doRPA = .false.
doRPAx = .false. doRPAx = .false.
docrRPA = .false.
doppRPA = .false. doppRPA = .false.
doG0F2 = .false. doG0F2 = .false.
@ -105,19 +107,20 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Read CC methods ! Read CC methods
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4 read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') doCCD = .true. if(answer1 == 'T') doCCD = .true.
if(answer2 == 'T') doDCD = .true. if(answer2 == 'T') dopCCD = .true.
if(answer3 == 'T') doCCSD = .true. if(answer3 == 'T') doDCD = .true.
if(answer4 == 'T') doCCSDT = .true. if(answer4 == 'T') doCCSD = .true.
if(answer5 == 'T') doCCSDT = .true.
! Read weird CC methods ! Read weird CC methods
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4 read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') do_drCCD = .true. if(answer1 == 'T') do_drCCD = .true.
if(answer2 == 'T') do_rCCD = .true. if(answer2 == 'T') do_rCCD = .true.
if(answer3 == 'T') do_lCCD = .true. if(answer3 == 'T') do_crCCD = .true.
if(answer4 == 'T') do_pCCD = .true. if(answer4 == 'T') do_lCCD = .true.
! Read excited state methods ! Read excited state methods
@ -131,10 +134,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
if(doCIS_D) doCIS = .true. if(doCIS_D) doCIS = .true.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3 read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doRPA = .true. if(answer1 == 'T') doRPA = .true.
if(answer2 == 'T') doRPAx = .true. if(answer2 == 'T') doRPAx = .true.
if(answer3 == 'T') doppRPA = .true. if(answer3 == 'T') docrRPA = .true.
if(answer4 == 'T') doppRPA = .true.
! Read Green function methods ! Read Green function methods

8
src/RPA/RPA.f90
View File

@ -87,12 +87,12 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,
endif endif
if(exchange_kernel) then ! if(exchange_kernel) then
EcRPA(1) = 0.5d0*EcRPA(1) ! EcRPA(1) = 0.5d0*EcRPA(1)
EcRPA(2) = 1.5d0*EcRPA(2) ! EcRPA(2) = 1.5d0*EcRPA(2)
end if ! end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

4
src/RPA/RPAx.f90
View File

@ -89,12 +89,12 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR
endif endif
if(exchange_kernel) then ! if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2) EcRPAx(2) = 1.5d0*EcRPAx(2)
end if ! end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'

26
src/RPA/ehRPA.f90 → src/RPA/crRPA.f90
View File

@ -1,7 +1,7 @@
subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI,dipole_int,eHF) ERI,dipole_int,eHF)
! Perform random phase approximation calculation with exchange (aka TDHF) ! Crossed-ring channel of the random phase approximation
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -42,7 +42,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
write(*,*) write(*,*)
write(*,*)'***********************************************************' write(*,*)'***********************************************************'
write(*,*)'| Random phase approximation calculation: eh channel |' write(*,*)'| Random phase approximation calculation: cr channel |'
write(*,*)'***********************************************************' write(*,*)'***********************************************************'
write(*,*) write(*,*)
@ -70,7 +70,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
@ -83,24 +83,24 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, &
EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif endif
if(exchange_kernel) then ! if(exchange_kernel) then
EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(1) = 0.5d0*EcRPAx(1)
EcRPAx(2) = 1.5d0*EcRPAx(2) EcRPAx(2) = 1.5d0*EcRPAx(2)
end if ! end if
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (singlet) =',EcRPAx(1) write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPAx(1)
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (triplet) =',EcRPAx(2) write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy =',EcRPAx(1) + EcRPAx(2) write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2)
write(*,'(2X,A50,F20.10)') 'Tr@ehRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,*) write(*,*)
@ -109,7 +109,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
! if(doACFDT) then ! if(doACFDT) then
! write(*,*) '-------------------------------------------------------' ! write(*,*) '-------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of ehRPA correlation energy' ! write(*,*) 'Adiabatic connection version of crRPA correlation energy'
! write(*,*) '-------------------------------------------------------' ! write(*,*) '-------------------------------------------------------'
! write(*,*) ! write(*,*)
@ -127,4 +127,4 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n
! end if ! end if
end subroutine ehRPA end subroutine crRPA