diff --git a/input/methods b/input/methods index c22322a..6331d53 100644 --- a/input/methods +++ b/input/methods @@ -2,14 +2,14 @@ T F F F # MP2* MP3 MP2-F12 F F F -# CCD DCD CCSD CCSD(T) - F F F F -# drCCD rCCD lCCD pCCD - T T F F +# CCD pCCD DCD CCSD CCSD(T) + F F F F F +# drCCD rCCD crCCD lCCD + T T T T # CIS* CIS(D) CID CISD FCI F F F F F -# RPA* RPAx* ppRPA - T T F +# RPA* RPAx* crRPA ppRPA + T T T T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/src/CC/crCCD.f90 b/src/CC/crCCD.f90 index 36cd10d..eb2e0f7 100644 --- a/src/CC/crCCD.f90 +++ b/src/CC/crCCD.f90 @@ -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 :: delta_OOVV(:,:,:,:) - double precision,allocatable :: OOOO(:,:,:,:) double precision,allocatable :: OOVV(:,:,:,:) 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 :: t2(:,:,:,:) @@ -66,7 +53,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF write(*,*) write(*,*)'**************************************' - write(*,*)'| crossed-ring CCD calculation |' + write(*,*)'| Crossed-ring CCD calculation |' write(*,*)'**************************************' write(*,*) @@ -105,13 +92,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! 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) - 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) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas) deallocate(dbERI) @@ -129,9 +113,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Initialization - allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(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)) + allocate(r2(nO,nO,nV,nV)) Conv = 1d0 nSCF = 0 @@ -159,23 +141,9 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! 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) - -! 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(:,:,:,:) + r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) ! Check convergence @@ -183,7 +151,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Update amplitudes - t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + t2(:,:,:,:) = t2(:,:,:,:) + r2(:,:,:,:)/delta_OOVV(:,:,:,:) ! Compute correlation energy @@ -227,10 +195,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)' crossed-ring CCD energy ' + write(*,*)' crossed-ring CCD energy ' write(*,*)'----------------------------------------------------' - write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD - write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD + write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD write(*,*)'----------------------------------------------------' write(*,*) diff --git a/src/CC/drCCD.f90 b/src/CC/drCCD.f90 index 0d21165..063c4b0 100644 --- a/src/CC/drCCD.f90 +++ b/src/CC/drCCD.f90 @@ -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 :: OVVO(:,:,:,:) - double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) 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 - 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) 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) @@ -135,7 +133,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! 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(:,:,:,:) diff --git a/src/CC/form_crossed_ring_r.f90 b/src/CC/form_crossed_ring_r.f90 new file mode 100644 index 0000000..0009e38 --- /dev/null +++ b/src/CC/form_crossed_ring_r.f90 @@ -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 diff --git a/src/CC/form_ring_r.f90 b/src/CC/form_ring_r.f90 index 2757650..8b75138 100644 --- a/src/CC/form_ring_r.f90 +++ b/src/CC/form_ring_r.f90 @@ -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 @@ -9,7 +9,6 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) integer,intent(in) :: nC,nO,nV,nR double precision,intent(in) :: t2(nO,nO,nV,nV) 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) ! 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 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 diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index 346f715..14dcb0c 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -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 :: OVVO(:,:,:,:) - double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) 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 - 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) 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) @@ -143,7 +141,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! 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(:,:,:,:) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 984550f..4646fb4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -9,10 +9,10 @@ program QuAcK logical :: dostab logical :: doKS logical :: doMP2,doMP3,doMP2F12 - logical :: doCCD,doDCD,doCCSD,doCCSDT - logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD + logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD logical :: doCIS,doCIS_D,doCID,doCISD,doFCI - logical :: doRPA,doRPAx,doppRPA + logical :: doRPA,doRPAx,docrRPA,doppRPA logical :: doADC logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW @@ -91,8 +91,6 @@ program QuAcK double precision :: start_CISD ,end_CISD ,t_CISD double precision :: start_FCI ,end_FCI ,t_FCI 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_GF2 ,end_GF2 ,t_GF2 double precision :: start_GF3 ,end_GF3 ,t_GF3 @@ -160,17 +158,17 @@ program QuAcK ! Which calculations do you want to do? - call read_methods(doRHF,doUHF,doKS,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & + call read_methods(doRHF,doUHF,doKS,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + doRPA,doRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doufG0W0,doufGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read options for methods @@ -645,6 +643,24 @@ program QuAcK 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 !------------------------------------------------------------------------ @@ -662,32 +678,11 @@ program QuAcK 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 !------------------------------------------------------------------------ - if(do_pCCD) then + if(dopCCD) then call cpu_time(start_CCD) 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 - call cpu_time(start_RPAx) + call cpu_time(start_RPA) if(unrestricted) then 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) end if - call cpu_time(end_RPAx) + call cpu_time(end_RPA) - t_RPAx = end_RPAx - start_RPAx - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds' + t_RPA = end_RPA - start_RPA + 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(*,*) end if @@ -812,12 +823,12 @@ program QuAcK 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 cpu_time(end_ppRPA) + call cpu_time(end_RPA) - t_ppRPA = end_ppRPA - start_ppRPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_ppRPA,' seconds' + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 3b77e7c..8ddd885 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,14 +1,14 @@ -subroutine read_methods(doRHF,doUHF,doKS,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & +subroutine read_methods(doRHF,doUHF,doKS,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + doRPA,doRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doufG0W0,doufGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read desired methods @@ -19,10 +19,10 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: doRHF,doUHF,doKS,doMOM logical,intent(out) :: doMP2,doMP3,doMP2F12 - logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT - logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD + logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD 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) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(out) :: doG0T0,doevGT,doqsGT @@ -48,14 +48,15 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doMP2F12 = .false. doCCD = .false. + dopCCD = .false. doDCD = .false. doCCSD = .false. doCCSDT = .false. do_drCCD = .false. do_rCCD = .false. + do_crCCD = .false. do_lCCD = .false. - do_pCCD = .false. doCIS = .false. doCIS_D = .false. @@ -65,6 +66,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doRPA = .false. doRPAx = .false. + docrRPA = .false. doppRPA = .false. doG0F2 = .false. @@ -105,19 +107,20 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read CC methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') doCCD = .true. - if(answer2 == 'T') doDCD = .true. - if(answer3 == 'T') doCCSD = .true. - if(answer4 == 'T') doCCSDT = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5 + if(answer1 == 'T') doCCD = .true. + if(answer2 == 'T') dopCCD = .true. + if(answer3 == 'T') doDCD = .true. + if(answer4 == 'T') doCCSD = .true. + if(answer5 == 'T') doCCSDT = .true. ! Read weird CC methods read(1,*) read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') do_drCCD = .true. if(answer2 == 'T') do_rCCD = .true. - if(answer3 == 'T') do_lCCD = .true. - if(answer4 == 'T') do_pCCD = .true. + if(answer3 == 'T') do_crCCD = .true. + if(answer4 == 'T') do_lCCD = .true. ! Read excited state methods @@ -131,10 +134,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & if(doCIS_D) doCIS = .true. read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') doRPA = .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 diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 index 3d6cb63..0dce1a7 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/RPA.f90 @@ -87,12 +87,12 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, endif - if(exchange_kernel) then +! if(exchange_kernel) then - EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(2) +! EcRPA(1) = 0.5d0*EcRPA(1) +! EcRPA(2) = 1.5d0*EcRPA(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/RPAx.f90 b/src/RPA/RPAx.f90 index 6b096e4..a7be80f 100644 --- a/src/RPA/RPAx.f90 +++ b/src/RPA/RPAx.f90 @@ -89,12 +89,12 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR endif - if(exchange_kernel) then +! if(exchange_kernel) then EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(2) = 1.5d0*EcRPAx(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/ehRPA.f90 b/src/RPA/crRPA.f90 similarity index 84% rename from src/RPA/ehRPA.f90 rename to src/RPA/crRPA.f90 index 103f5c3..38cb004 100644 --- a/src/RPA/ehRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -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) -! Perform random phase approximation calculation with exchange (aka TDHF) +! Crossed-ring channel of the random phase approximation implicit none include 'parameters.h' @@ -42,7 +42,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) write(*,*)'***********************************************************' - write(*,*)'| Random phase approximation calculation: eh channel |' + write(*,*)'| Random phase approximation calculation: cr channel |' 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, & 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)) 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, & 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)) endif - if(exchange_kernel) then +! if(exchange_kernel) then EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(2) = 1.5d0*EcRPAx(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (singlet) =',EcRPAx(1) - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (triplet) =',EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA 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 correlation energy (singlet) =',EcRPAx(1) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -109,7 +109,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! if(doACFDT) then ! write(*,*) '-------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of ehRPA correlation energy' +! write(*,*) 'Adiabatic connection version of crRPA correlation energy' ! write(*,*) '-------------------------------------------------------' ! write(*,*) @@ -127,4 +127,4 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! end if -end subroutine ehRPA +end subroutine crRPA