diff --git a/input/methods b/input/methods index f719fe3..4673c55 100644 --- a/input/methods +++ b/input/methods @@ -6,12 +6,12 @@ F F F F # drCCD rCCD lCCD pCCD F F F F -# CIS* CIS(D) CID CISD - F F F F +# CIS* CIS(D) CID CISD FCI + F F F F F # RPA* RPAx* ppRPA F F F -# G0F2 evGF2 qsGF2* G0F3 evGF3 - F F T F F +# G0F2* evGF2 qsGF2* G0F3 evGF3 + T F F F F # G0W0* evGW* qsGW* F F F # G0T0 evGT qsGT diff --git a/input/options b/input/options index f64a274..d2044fb 100644 --- a/input/options +++ b/input/options @@ -7,7 +7,7 @@ # spin: TDA singlet triplet spin_conserved spin_flip F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm - 256 0.00001 T 5 T 0.001 3 + 256 0.00001 T 5 T 0.001 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 # ACFDT: AC Kx XBS diff --git a/mol/butadiene.xyz b/mol/butadiene.xyz index f2069d8..8bbbc73 100644 --- a/mol/butadiene.xyz +++ b/mol/butadiene.xyz @@ -1,5 +1,5 @@ 10 -Butadiene,^1A_g,CC3,aug-cc-pVTZ + C 0.60673471 0.00000000 0.39936380 C -0.60673471 0.00000000 -0.39936380 C 1.84223863 0.00000000 -0.11897388 @@ -9,4 +9,4 @@ H -0.48033933 0.00000000 -1.47579018 H 1.99778057 0.00000000 -1.19009558 H -1.99778057 0.00000000 1.19009558 H 2.71819794 0.00000000 0.51257105 -H -2.71819794 0.00000000 -0.51257105 \ No newline at end of file +H -2.71819794 0.00000000 -0.51257105 diff --git a/src/CI/FCI.f90 b/src/CI/FCI.f90 new file mode 100644 index 0000000..1fcabbd --- /dev/null +++ b/src/CI/FCI.f90 @@ -0,0 +1,31 @@ +subroutine FCI(nBas,nC,nO,nV,nR,ERI,e) + +! Perform a full configuration interaction calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + +! Hello world + + write(*,*) + write(*,*)'**********************************' + write(*,*)'| Full Configuration Interaction |' + write(*,*)'**********************************' + write(*,*) + +! Form FCI vector + +! Form FCI matrix + +! Diagonalize FCI matrix + + +end subroutine FCI diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 new file mode 100644 index 0000000..8f64d4e --- /dev/null +++ b/src/GF/UG0F2.f90 @@ -0,0 +1,122 @@ +subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF) + +! Perform unrestricted G0W0 calculation + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: BSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + logical,intent(in) :: linearize + double precision,intent(in) :: eta + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas,nspin) + +! Local variables + + integer :: is + integer :: ispin + double precision :: Ec(nsp) + double precision :: EcBSE(nspin) + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: Z(:,:) + integer :: nS_aa,nS_bb,nS_sc + + double precision,allocatable :: eGF2lin(:,:) + double precision,allocatable :: eGF2(:,:) + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot G0F2 calculation |' + write(*,*)'| *** Unrestricted version *** |' + write(*,*)'************************************************' + write(*,*) + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(SigC(nBas,nspin),Z(nBas,nspin),eGF2(nBas,nspin),eGF2lin(nBas,nspin)) + +!---------------------! +! Compute self-energy ! +!---------------------! + + call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! + + eGF2lin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGF2(:,:) = eGF2lin(:,:) + + else + + ! Find graphical solution of the QP equation + + print*,'!!! Graphical solution NYI for UG0F2 !!!' + stop + + end if + +! Compute MP2 correlation energy + + call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eHF,Ec) + +! Dump results + + call print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec) + +! Perform BSE calculation + + if(BSE) then + + print*,'!!! BSE2 NYI for UG0F2 !!!' + + end if + +end subroutine UG0F2 diff --git a/src/GF/print_UG0F2.f90 b/src/GF/print_UG0F2.f90 new file mode 100644 index 0000000..1495401 --- /dev/null +++ b/src/GF/print_UG0F2.f90 @@ -0,0 +1,73 @@ +subroutine print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec) + +! Print one-electron energies and other stuff for G0W0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: Ec(nsp) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: SigC(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGF2(nBas,nspin) + + integer :: p + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + +! HOMO and LUMO + + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGF2(nO(ispin),ispin) + LUMO(ispin) = eGF2(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = eGF2(1,ispin) + Gap(ispin) = 0d0 + end if + end do + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*)' Unrestricted one-shot G0F2 calculation (eV)' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & + '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'UG0F2 HOMO energy:',maxval(HOMO(:))*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'UG0F2 LUMO energy:',minval(LUMO(:))*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'UG0F2 HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') ' UG0F2 total energy :',ENuc + EUHF + sum(Ec(:)),' au' + write(*,'(2X,A30,F15.6,A3)') ' UG0F2 correlation energy:',sum(Ec(:)),' au' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*) + +end subroutine print_UG0F2 + + diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index eaa0dbb..70a8e68 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -306,4 +306,12 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) +! Perform BSE calculation + + if(BSE) then + + print*,'!!! BSE2 NYI for qsUGF2 !!!' + + end if + end subroutine qsUGF2 diff --git a/src/GF/unrestricted_self_energy_GF2.f90 b/src/GF/unrestricted_self_energy_GF2.f90 index bbd9e72..7bc83dd 100644 --- a/src/GF/unrestricted_self_energy_GF2.f90 +++ b/src/GF/unrestricted_self_energy_GF2.f90 @@ -35,6 +35,8 @@ subroutine unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_b ! Compute self-energy | !---------------------! + SigC(:,:,:) = 0d0 + !----------------! ! Spin-up sector !----------------! diff --git a/src/GF/unrestricted_self_energy_GF2_diag.f90 b/src/GF/unrestricted_self_energy_GF2_diag.f90 index d88c1eb..6f06010 100644 --- a/src/GF/unrestricted_self_energy_GF2_diag.f90 +++ b/src/GF/unrestricted_self_energy_GF2_diag.f90 @@ -35,6 +35,8 @@ subroutine unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab, ! Compute self-energy | !---------------------! + SigC(:,:) = 0d0 + !----------------! ! Spin-up sector !----------------! diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index db7f56d..7802f04 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -207,6 +207,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) + call matout(nBas,nBAs,SigCp) + ! Compute commutator and convergence criteria error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 64e4819..1b191d4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -11,7 +11,7 @@ program QuAcK logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doDCD,doCCSD,doCCSDT logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD - logical :: doCIS,doCIS_D,doCID,doCISD + logical :: doCIS,doCIS_D,doCID,doCISD,doFCI logical :: doRPA,doRPAx,doppRPA logical :: doADC logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 @@ -89,6 +89,7 @@ program QuAcK double precision :: start_CIS ,end_CIS ,t_CIS double precision :: start_CID ,end_CID ,t_CID 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 @@ -162,7 +163,7 @@ program QuAcK doMP2,doMP3,doMP2F12, & doCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCIS_D,doCID,doCISD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doqsGF2, & doG0F3,doevGF3, & @@ -812,8 +813,19 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + + if(unrestricted) then + + call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF) + + else + + call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & + eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + + end if + call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -1176,6 +1188,22 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Compute FCI +!------------------------------------------------------------------------ + + if(doFCI) then + + call cpu_time(start_FCI) + call FCI(nBas,nC,nO,nV,nR,ERI_MO,eHF) + call cpu_time(end_FCI) + + t_FCI = end_FCI - start_FCI + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_FCI,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! End of QuAcK !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 1401cc9..5228e17 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -2,7 +2,7 @@ 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, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doqsGF2, & doG0F3,doevGF3, & @@ -20,7 +20,7 @@ subroutine read_methods(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) :: doCIS,doCIS_D,doCID,doCISD + logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI logical,intent(out) :: doRPA,doRPAx,doppRPA logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW @@ -60,6 +60,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doCIS_D = .false. doCID = .false. doCISD = .false. + doFCI = .false. doRPA = .false. doRPAx = .false. @@ -118,11 +119,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read excited state methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 + read(1,*) answer1,answer2,answer3,answer4,answer5 if(answer1 == 'T') doCIS = .true. if(answer2 == 'T') doCIS_D = .true. if(answer3 == 'T') doCID = .true. if(answer4 == 'T') doCISD = .true. + if(answer5 == 'T') doFCI = .true. if(doCIS_D) doCIS = .true. read(1,*)