diff --git a/input/basis b/input/basis index b2b2293..c7ec793 100644 --- a/input/basis +++ b/input/basis @@ -1,30 +1,38 @@ 1 6 S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 + 1 6665.0000000 0.0006920 + 2 1000.0000000 0.0053290 + 3 228.0000000 0.0270770 + 4 64.7100000 0.1017180 + 5 21.0600000 0.2747400 + 6 7.4950000 0.4485640 + 7 2.7970000 0.2850740 + 8 0.5215000 0.0152040 S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 + 1 6665.0000000 -0.0001460 + 2 1000.0000000 -0.0011540 + 3 228.0000000 -0.0057250 + 4 64.7100000 -0.0233120 + 5 21.0600000 -0.0639550 + 6 7.4950000 -0.1499810 + 7 2.7970000 -0.1272620 + 8 0.5215000 0.5445290 S 1 - 1 0.0280500 1.0000000 + 1 0.1596000 1.0000000 P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 + 1 9.4390000 0.0381090 + 2 2.0020000 0.2094800 + 3 0.5456000 0.5085570 P 1 - 1 0.0240300 1.0000000 + 1 0.1517000 1.0000000 D 1 - 1 0.1239000 1.0000000 - + 1 0.5500000 1.0000000 +2 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +P 1 + 1 0.7270000 1.0000000 diff --git a/input/methods b/input/methods index 70cc161..81c1b6d 100644 --- a/input/methods +++ b/input/methods @@ -1,20 +1,21 @@ # RHF UHF MOM F T F -# MP2 MP3 MP2-F12 +# MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD F F F F -# CIS CID CISD - T F F -# RPA RPAx ppRPA +# CIS* CID CISD + F F F +# RPA* RPAx* ppRPA F F F # G0F2 evGF2 G0F3 evGF3 F F F F -# G0W0 evGW qsGW - F F F +# G0W0* evGW* qsGW + T T F # G0T0 evGT qsGT F F F # MCMP2 F +# * unrestricted version available diff --git a/input/molecule b/input/molecule index 058d6dd..e2a4fd3 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 2 1 0 0 + 2 4 3 0 0 # Znuc x y z - Li 0.0 0.0 0.0 + C 0. 0. -0.16245872 + H 0. 0. 1.93436816 diff --git a/input/molecule.xyz b/input/molecule.xyz index c9a5a65..7a4f218 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,4 @@ - 1 + 2 - Li 0.0000000000 0.0000000000 0.0000000000 + C 0.0000000000 0.0000000000 -0.0859694585 + H 0.0000000000 0.0000000000 1.0236236215 diff --git a/input/options b/input/options index cb56d2c..0d3bf88 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # RHF: maxSCF thresh DIIS n_diis guess_type ortho_type - 64 0.0000001 T 5 1 1 + 64 0.00001 T 5 1 1 # MP: # CC: maxSCF thresh DIIS n_diis @@ -9,10 +9,10 @@ # 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 F 0.0 F F F F F + 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS - T F T + F F T # BSE: BSE dBSE dTDA evDyn - T T T T + F F T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4a63fe3..f7310a1 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -823,9 +823,18 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) + if(unrestricted) then + + call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & + ERHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,PHF,cHF,eHF,eG0W0) + + else + + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) + end if call cpu_time(end_evGW) t_evGW = end_evGW - start_evGW diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 5df2b99..c13a2a9 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -179,9 +179,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - eHF,eGW,EcBSE) + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,eHF,eGW,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 79f2dc2..d7d816e 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,6 +1,6 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H, & - ERI,PHF,cHF,eHF,eG0W0) +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI, & + PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -36,7 +36,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eG0W0(nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: H(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables diff --git a/src/QuAcK/evUGW.f90 b/src/QuAcK/evUGW.f90 new file mode 100644 index 0000000..3c8b2d8 --- /dev/null +++ b/src/QuAcK/evUGW.f90 @@ -0,0 +1,307 @@ +subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & + ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,PHF,cHF,eHF,eG0W0) + +! Perform self-consistent eigenvalue-only GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: COHSEX + logical,intent(in) :: BSE + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: G0W + logical,intent(in) :: GW0 + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + 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) :: eHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: PHF(nBas,nBas,nspin) + double precision,intent(in) :: eG0W0(nBas,nspin) + double precision,intent(in) :: Hc(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) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + logical :: linear_mixing + integer :: is + integer :: ispin + integer :: nSCF + integer :: n_diis + double precision :: rcond(nspin) + double precision :: Conv + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision :: alpha + double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: e_diis(:,:,:) + double precision,allocatable :: eGW(:,:) + double precision,allocatable :: eOld(:,:) + double precision,allocatable :: Z(:,:) + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent evGW calculation |' + write(*,*)'************************************************' + write(*,*) + +! COHSEX approximation + + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! GW0 + + if(GW0) then + write(*,*) 'GW0 scheme activated!' + write(*,*) + end if + +! G0W + + if(G0W) then + write(*,*) 'G0W scheme activated!' + write(*,*) + end if + +! Linear mixing + + linear_mixing = .false. + alpha = 0.2d0 + +! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc), & + XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), & + error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + +! Initialization + + nSCF = 0 + ispin = 1 + n_diis = 0 + Conv = 1d0 + e_diis(:,:,:) = 0d0 + error_diis(:,:,:) = 0d0 + eGW(:,:) = eG0W0(:,:) + eOld(:,:) = eGW(:,:) + Z(:,:) = 1d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Compute screening + + if(.not. GW0 .or. nSCF == 0) then + + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + + endif + + !----------------------! + ! Excitation densities ! + !----------------------! + + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + !------------------------------------------------! + ! Compute self-energy and renormalization factor ! + !------------------------------------------------! + + if(G0W) then + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + endif + + !-----------------------------------! + ! Solve the quasi-particle equation ! + !-----------------------------------! + + eGW(:,:) = eHF(:,:) + SigC(:,:) + + ! Convergence criteria + + Conv = maxval(abs(eGW(:,:) - eOld(:,:))) + + ! Print results + + call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) + + ! Linear mixing or DIIS extrapolation + + if(linear_mixing) then + + eGW(:,:) = alpha*eGW(:,:) + (1d0 - alpha)*eOld(:,:) + + else + + n_diis = min(n_diis+1,max_diis) + do is=1,nspin + call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis,error_diis(:,:,is), & + e_diis(:,:,is),eGW(:,is)-eOld(:,is),eGW(:,is)) + end do + +! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + + endif + + ! Save quasiparticles energy for next cycle + + eOld(:,:) = eGW(:,:) + + ! Increment + + nSCF = nSCF + 1 + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Plot stuff + +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis) + +! Perform BSE calculation + + if(BSE) then + + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,eGW,eGW,EcBSE) + +! if(exchange_kernel) then + +! EcBSE(1) = 0.5d0*EcBSE(1) +! EcBSE(2) = 1.5d0*EcBSE(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of BSE correlation energy' +! write(*,*) '------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + + endif + +end subroutine evUGW diff --git a/src/QuAcK/print_evUGW.f90 b/src/QuAcK/print_evUGW.f90 new file mode 100644 index 0000000..7c32002 --- /dev/null +++ b/src/QuAcK/print_evUGW.f90 @@ -0,0 +1,72 @@ +subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) + +! Print one-electron energies and other stuff for evGW + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nSCF + double precision,intent(in) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: EcRPA + double precision,intent(in) :: Conv + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: SigC(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGW(nBas,nspin) + + integer :: p + double precision :: HOMO + double precision :: LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = max(eGW(nO(1),1),eGW(nO(2),2)) + LUMO = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) + Gap = LUMO - HOMO + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation' + endif + 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,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',HOMO*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',LUMO*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*) + +end subroutine print_evUGW