diff --git a/input/methods b/input/methods index db30272..f31cdf6 100644 --- a/input/methods +++ b/input/methods @@ -2,10 +2,10 @@ T F F F # MP2* MP3 MP2-F12 F F F -# CCD CCSD CCSD(T) - F F F +# CCD DCD CCSD CCSD(T) + T T F F # drCCD rCCD lCCD pCCD - F F F T + F F F F # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA diff --git a/input/options b/input/options index 4e7d31d..c63b7ab 100644 --- a/input/options +++ b/input/options @@ -9,9 +9,9 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.0 F F F F F + 256 0.00001 T 5 T 0.00367493 F F F F F # ACFDT: AC Kx XBS - F F T + T F T # BSE: BSE dBSE dTDA evDyn F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/CC/DCD.f90 b/src/CC/DCD.f90 new file mode 100644 index 0000000..aea4af2 --- /dev/null +++ b/src/CC/DCD.f90 @@ -0,0 +1,283 @@ +subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! DCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: EDCD,EcDCD + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVOV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: xO(:,:) + double precision,allocatable :: xV(:,:) + double precision,allocatable :: xOV(:,:,:,:) + + double precision,allocatable :: r(:,:,:,:) + double precision,allocatable :: t(:,:,:,:) + double precision,allocatable :: tt(:,:,:,:) + + integer :: n_diis + integer :: i,j,k,l + integer :: a,b,c,d + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| DCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Form energy denominator + + allocate(eO(nO-nC),eV(nV-nR)) + allocate(delta(nO-nC,nO-nC,nV-nR,nV-nR)) + + eO(:) = eHF(nC+1:nO) + eV(:) = eHF(nO+1:nBas-nR) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta) + +! Create integral batches + + allocate(OOOO(nO-nC,nO-nC,nO-nC,nO-nC), & + OOVV(nO-nC,nO-nC,nV-nR,nV-nR), & + OVOV(nO-nC,nV-nR,nO-nC,nV-nR), & + OVVO(nO-nC,nV-nR,nV-nR,nO-nC), & + VVVV(nV-nR,nV-nR,nV-nR,nV-nR)) + + OOOO(:,:,:,:) = ERI(nC+1:nO ,nC+1:nO ,nC+1:nO ,nC+1:nO ) + OOVV(:,:,:,:) = ERI(nC+1:nO ,nC+1:nO ,nO+1:nBas-nR,nO+1:nBas-nR) + OVOV(:,:,:,:) = ERI(nC+1:nO ,nO+1:nBas-nR,nC+1:nO ,nO+1:nBas-nR) + OVVO(:,:,:,:) = ERI(nC+1:nO ,nO+1:nBas-nR,nO+1:nBas-nR,nC+1:nO ) + VVVV(:,:,:,:) = ERI(nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR) + +! MP2 guess amplitudes + + allocate(t(nO-nC,nO-nC,nV-nR,nV-nR)) + + t(:,:,:,:) = -OOVV(:,:,:,:)/delta(:,:,:,:) + +! Memory allocation for DIIS + + allocate(error_diis((nO-nR)**2*(nV-nR)**2,max_diis), & + t_diis((nO-nR)**2*(nV-nR)**2,max_diis)) + +! Initialization + + allocate(tt(nO-nC,nO-nC,nV-nR,nV-nR), & + r(nO-nC,nO-nC,nV-nR,nV-nR), & + xO(nO-nC,nO-nC),xV(nV-nR,nV-nR), & + xOV(nO-nC,nO-nC,nV-nR,nV-nR)) + + Conv = 1d0 + nSCF = 0 + + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| DCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(DCD)','|','Ec(DCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Compute correlation energy + + EcDCD = 0d0 + do i=1,nO-nC + do j=1,nO-nC + do a=1,nV-nR + do b=1,nV-nR + + EcDCD = EcDCD + (2d0*OOVV(i,j,a,b) - OOVV(i,j,b,a))*t(i,j,a,b) + + enddo + enddo + enddo + enddo + + ! Dump results + + EDCD = ERHF + EcDCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',EDCD+ENuc,'|',EcDCD,'|',Conv,'|' + + ! Increment + + nSCF = nSCF + 1 + + ! Form interemediate arrays + + do i=1,nO-nC + do j=1,nO-nC + do a=1,nV-nR + do b=1,nV-nR + tt(i,j,a,b) = 2d0*t(i,j,a,b) - t(i,j,b,a) + enddo + enddo + enddo + enddo + + xV(:,:) = 0d0 + do a=1,nV-nR + do b=1,nV-nR + do k=1,nO-nC + do l=1,nO-nC + do c=1,nV-nR + xV(a,b) = xV(a,b) + tt(k,l,a,c)*OOVV(l,k,c,b) + end do + end do + end do + end do + end do + + xO(:,:) = 0d0 + do i=1,nO-nC + do j=1,nO-nC + do k=1,nO-nC + do c=1,nV-nR + do d=1,nV-nR + xO(i,j) = xO(i,j) + tt(i,k,c,d)*OOVV(k,j,d,c) + end do + end do + end do + end do + end do + + xOV(:,:,:,:) = 0d0 + do i=1,nO-nC + do j=1,nO-nC + do a=1,nV-nR + do b=1,nV-nR + do k=1,nO-nC + do c=1,nV-nR + xOV(i,j,a,b) = xOV(i,j,a,b) + tt(i,k,a,c)*OOVV(k,j,c,b) + end do + end do + end do + end do + end do + end do + + ! Compute residual + + r(:,:,:,:) = OOVV(:,:,:,:) + delta(:,:,:,:)*t(:,:,:,:) + + do i=1,nO-nC + do j=1,nO-nC + do a=1,nV-nR + do b=1,nV-nR + + do c=1,nV-nR + r(i,j,a,b) = r(i,j,a,b) - 0.5d0*xV(a,c)*t(i,j,c,b) & + - 0.5d0*xV(b,c)*t(j,i,c,a) + do d=1,nV-nR + r(i,j,a,b) = r(i,j,a,b) + VVVV(a,b,c,d)*t(i,j,c,d) + end do + end do + + do k=1,nO-nC + r(i,j,a,b) = r(i,j,a,b) - 0.5d0*xO(k,i)*t(k,j,a,b) & + - 0.5d0*xO(k,j)*t(k,i,b,a) + do l=1,nO-nC + r(i,j,a,b) = r(i,j,a,b) + OOOO(i,j,k,l)*t(k,l,a,b) + end do + end do + + do k=1,nO-nC + do c=1,nV-nR + + r(i,j,a,b) = r(i,j,a,b) + xOV(i,k,a,c)*tt(k,j,c,b) & + - OVOV(k,a,i,c)*t(k,j,c,b) & + - OVOV(k,b,i,c)*t(k,j,a,c) & + - OVOV(k,b,j,c)*t(k,i,c,a) & + - OVOV(k,a,j,c)*t(k,i,b,c) & + + tt(i,k,a,c)*OVVO(k,b,c,j) & + + tt(j,k,b,c)*OVVO(k,a,c,i) + end do + end do + + end do + end do + end do + end do + + ! Check convergence + + Conv = maxval(abs(r(:,:,:,:))) + + ! Update amplitudes + + t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta(:,:,:,:) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,(nO-nC)**2*(nV-nR)**2,(nO-nC)**2*(nV-nR)**2,n_diis,error_diis,t_diis,-r/delta,t) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +end subroutine DCD diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 395f3e9..e0a75d6 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -8,7 +8,7 @@ program QuAcK logical :: doRHF,doUHF,doMOM logical :: doKS logical :: doMP2,doMP3,doMP2F12 - logical :: doCCD,doCCSD,doCCSDT + logical :: doCCD,doDCD,doCCSD,doCCSDT logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical :: doCIS,doCIS_D,doCID,doCISD logical :: doRPA,doRPAx,doppRPA @@ -81,6 +81,7 @@ program QuAcK double precision :: start_MOM ,end_MOM ,t_MOM double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO double precision :: start_CCD ,end_CCD ,t_CCD + double precision :: start_DCD ,end_DCD ,t_DCD double precision :: start_CCSD ,end_CCSD ,t_CCSD double precision :: start_CIS ,end_CIS ,t_CIS double precision :: start_CID ,end_CID ,t_CID @@ -156,7 +157,7 @@ program QuAcK call read_methods(doRHF,doUHF,doKS,doMOM, & doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & + doCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & @@ -528,6 +529,23 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform DCD calculation +!------------------------------------------------------------------------ + + if(doDCD) then + + call cpu_time(start_DCD) + call DCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & + ERI_MO,ENuc,ERHF,eHF) + call cpu_time(end_DCD) + + t_DCD = end_DCD - start_DCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_DCD,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Perform CCSD or CCSD(T) calculation !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 7914233..ab0d2d6 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,6 +1,6 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & + doCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & @@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: doRHF,doUHF,doKS,doMOM logical,intent(out) :: doMP2,doMP3,doMP2F12 - logical,intent(out) :: doCCD,doCCSD,doCCSDT + logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD logical,intent(out) :: doRPA,doRPAx,doppRPA @@ -46,6 +46,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doMP2F12 = .false. doCCD = .false. + doDCD = .false. doCCSD = .false. doCCSDT = .false. @@ -98,10 +99,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read CC methods read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') doCCD = .true. - if(answer2 == 'T') doCCSD = .true. - if(answer3 == 'T') doCCSDT = .true. + if(answer2 == 'T') doDCD = .true. + if(answer3 == 'T') doCCSD = .true. + if(answer4 == 'T') doCCSDT = .true. ! Read weird CC methods read(1,*)