diff --git a/input/methods b/input/methods index 2c30f9f..52dabfe 100644 --- a/input/methods +++ b/input/methods @@ -1,11 +1,11 @@ # RHF UHF MOM - F T F -# MP2* MP3 MP2-F12 T F F +# MP2* MP3 MP2-F12 + F F F # CCD CCSD CCSD(T) - F F F + T F F # drCCD rCCD lCCD pCCD - F F F F + F F F T # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 8abb1c0..be544b7 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + T T T T T # 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 @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS T F T # BSE: BSE dBSE dTDA evDyn - T F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/CC/form_delta_OV.f90 b/src/CC/form_delta_OV.f90 index 02b737d..a4a671a 100644 --- a/src/CC/form_delta_OV.f90 +++ b/src/CC/form_delta_OV.f90 @@ -7,8 +7,8 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta) ! Input variables integer,intent(in) :: nC,nO,nV,nR - double precision,intent(in) :: eO(nO) - double precision,intent(in) :: eV(nV) + double precision,intent(in) :: eO(nO-nC) + double precision,intent(in) :: eV(nV-nR) ! Local variables @@ -18,7 +18,7 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta) double precision,intent(out) :: delta(nO,nV) - do i=nC+1,nO + do i=1,nO-nC do a=1,nV-nR delta(i,a) = eV(a) - eO(i) enddo diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index dcfe7ab..29e4629 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -23,7 +23,9 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) double precision :: Conv double precision :: ECCD,EcCCD - double precision,allocatable :: delta_OOVV(:,:) + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OV(:,:) double precision,allocatable :: OOOO(:,:) double precision,allocatable :: OOVV(:,:) @@ -40,6 +42,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) double precision :: rcond double precision,allocatable :: error_diis(:,:) double precision,allocatable :: t_diis(:,:) + double precision,external :: trace_matrix ! Hello world @@ -51,36 +54,28 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Form energy denominator - allocate(delta_OOVV(nO,nV)) + allocate(eO(nO-nC),eV(nV-nR),delta_OV(nO-nC,nV-nR)) - delta_OOVV(:,:) = 0d0 + eO(:) = eHF(nC+1:nO) + eV(:) = eHF(nO+1:nBas-nR) - do i=nC+1,nO - do a=1,nV-nR - delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i)) - enddo - enddo + call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) ! Create integral batches - allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV)) + allocate(OOOO(nO-nC,nO-nC),OOVV(nO-nC,nV-nR),OVOV(nO-nC,nV-nR),OVVO(nO-nC,nV-nR),VVVV(nV-nR,nV-nR)) - OOOO(:,:) = 0d0 - OOVV(:,:) = 0d0 - OVVO(:,:) = 0d0 - VVVV(:,:) = 0d0 - - do i=nC+1,nO - do j=nC+1,nO - OOOO(i,j) = ERI(i,i,j,j) + do i=1,nO-nC + do j=1,nO-nC + OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j) end do end do - do i=nC+1,nO + do i=1,nO-nC do a=1,nV-nR - OOVV(i,a) = ERI(i,i,nO+a,nO+a) - OVOV(i,a) = ERI(i,nO+a,i,nO+a) - OVVO(i,a) = ERI(i,nO+a,nO+a,i) + OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a) + OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a) + OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i) end do end do @@ -92,17 +87,17 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! MP2 guess amplitudes - allocate(t(nO,nV)) + allocate(t(nO-nC,nV-nR)) - t(:,:) = - OOVV(:,:)/delta_OOVV(:,:) + t(:,:) = -0.5d0*OOVV(:,:)/delta_OV(:,:) ! Memory allocation for DIIS - allocate(error_diis(nO*nV,max_diis),t_diis(nO*nV,max_diis)) + allocate(error_diis((nO-nC)*(nV-nR),max_diis),t_diis((nO-nC)*(nV-nR),max_diis)) ! Initialization - allocate(r(nO,nV),y(nO,nO)) + allocate(r(nO-nC,nV-nR),y(nO-nC,nO-nC)) Conv = 1d0 nSCF = 0 @@ -130,24 +125,17 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Form intermediate array - y(:,:) = 0d0 - do i=nC+1,nO - do j=nC+1,nO - do b=1,nV-nR - y(i,j) = y(i,j) + OOVV(j,b)*t(i,b) - end do - end do - end do + y(:,:) = matmul(t,transpose(OOVV)) ! Compute residual - do i=nC+1,nO + do i=1,nO-nC do a=1,nV-nR - r(i,a) = OOVV(i,a) + delta_OOVV(i,a)*t(i,a) & + r(i,a) = OOVV(i,a) + 2d0*delta_OV(i,a)*t(i,a) & - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a) - do j=nC+1,nO + do j=1,nO-nC r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) end do @@ -164,16 +152,11 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Update amplitudes - t(:,:) = t(:,:) - r(:,:)/delta_OOVV(:,:) + t(:,:) = t(:,:) - 0.5d0*r(:,:)/delta_OV(:,:) ! Compute correlation energy - EcCCD = 0d0 - do i=nC+1,nO - do a=1,nV-nR - EcCCD = EcCCD + t(i,a)*OOVV(i,a) - end do - end do + EcCCD = trace_matrix(nO,matmul(t,transpose(OOVV))) ! Dump results @@ -182,7 +165,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t) + call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-0.5d0*r/delta_OV,t) ! Reset DIIS if required diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index a7a976d..a7a2f0f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -268,6 +268,14 @@ program QuAcK if(doRHF) then + ! Check that RHF calculation is worth doing... + + if(nO(1) /= nO(2)) then + write(*,*) ' !!! The system does not appear to be closed shell !!!' + write(*,*) + stop + end if + call cpu_time(start_HF) call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,ERI_AO,dipole_int,X,ERHF,eHF,cHF,PHF) diff --git a/src/modules/.gitignore b/src/modules/.gitignore index e69de29..63a7748 100644 --- a/src/modules/.gitignore +++ b/src/modules/.gitignore @@ -0,0 +1 @@ +*.mod