diff --git a/input/methods b/input/methods index d1b5eb7..af2bd81 100644 --- a/input/methods +++ b/input/methods @@ -1,9 +1,9 @@ # RHF UHF KS MOM - T F F F + F T F F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) - T F F + F F F # drCCD rCCD lCCD pCCD F F F F # CIS* CIS(D) CID CISD @@ -12,8 +12,8 @@ F F F # G0F2 evGF2 G0F3 evGF3 F F F F -# G0W0* evGW* qsGW - F F F +# G0W0* evGW* qsGW* + F F T # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 0094fa0..93da02b 100644 --- a/input/options +++ b/input/options @@ -5,11 +5,11 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T F T T + F 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 - 256 0.00001 T 5 T 0.0 F F F F F + 128 0.00001 T 5 T 0.00367493 F F T F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/mol/ch.xyz b/mol/ch.xyz index f2e0bda..9fdbfc8 100644 --- a/mol/ch.xyz +++ b/mol/ch.xyz @@ -1,4 +1,4 @@ 2 -\ce{CH},^2\Pi,CC3,aug-cc-pVTZ + C 0.00000000 0.00000000 -0.08596945 -H 0.00000000 0.00000000 1.02362355 \ No newline at end of file +H 0.00000000 0.00000000 1.02362355 diff --git a/mol/h2.xyz b/mol/h2.xyz index 494db3f..a302682 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 5.0 +H 0.0 0.0 1.4 diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 165330a..558bafd 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -174,8 +174,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T n_diis = min(n_diis+1,max_diis) do ispin=1,nspin - if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin), & - err(:,:,ispin),F(:,:,ispin)) + if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), & + F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin)) end do ! Reset DIIS if required diff --git a/src/MBPT/evUGW.f90 b/src/MBPT/evUGW.f90 index afae05e..a4f1529 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/MBPT/evUGW.f90 @@ -153,8 +153,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - - endif !----------------------! @@ -203,8 +201,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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)) + call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis,error_diis(:,1:n_diis,is), & + e_diis(:,1:n_diis,is),eGW(:,is)-eOld(:,is),eGW(:,is)) end do ! Reset DIIS if required diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index d57fc00..65fdcca 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW) +subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EqsGW) ! Print one-electron energies and other stuff for qsGW @@ -9,7 +9,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z ! Input variables integer,intent(in) :: nBas,nO,nSCF - double precision,intent(in) :: ENuc,EcRPA,EcGM,Conv,thresh + double precision,intent(in) :: ENuc,EcRPA,Conv,thresh double precision,intent(in) :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) @@ -35,8 +35,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z EV = trace_matrix(nBas,matmul(P,V)) EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) -! Ec = 0.5d0*trace_matrix(nBas,matmul(P,SigC)) - Ec = 0d0 + Ec = 0.5d0*trace_matrix(nBas,matmul(P,SigC)) EqsGW = ET + EV + EJ + Ex + Ec ! Dump results @@ -57,7 +56,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z '|',x,'|',eHF(x)*HaToeV,'|',(eGW(x)-eHF(x))*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' enddo - write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv @@ -67,11 +65,9 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z write(*,'(2X,A30,F15.6)') 'qsGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------' write(*,'(2X,A30,F15.6)') 'qsGW total energy =',EqsGW + ENuc - write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA - write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 new file mode 100644 index 0000000..38aa146 --- /dev/null +++ b/src/MBPT/print_qsUGW.f90 @@ -0,0 +1,204 @@ +subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,dipole) + +! Print one-electron energies and other stuff for qsUGW + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: Ov(nBas,nBas) + 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) :: thresh + double precision,intent(in) :: eGW(nBas,nspin) + double precision,intent(in) :: cGW(nBas,nBas,nspin) + double precision,intent(in) :: PGW(nBas,nBas,nspin) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: J(nBas,nBas,nspin) + double precision,intent(in) :: K(nBas,nBas,nspin) + double precision,intent(in) :: SigC(nBas,nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: p + integer :: ispin,ixyz + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + double precision :: ET(nspin) + double precision :: EV(nspin) + double precision :: EJ(nsp) + double precision :: Ex(nspin) + double precision :: Ec(nsp) + double precision :: EqsGW + double precision :: S_exact,S2_exact + double precision :: S,S2 + double precision,external :: trace_matrix + +! HOMO and LUMO + + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGW(nO(ispin),ispin) + LUMO(ispin) = eGW(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = eGW(1,ispin) + Gap(ispin) = 0d0 + end if + end do + + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2 = S2_exact + nO(2) - sum(matmul(transpose(cGW(:,1:nO(1),1)),matmul(Ov,cGW(:,1:nO(2),2)))**2) + + S_exact = 0.5d0*dble(nO(1) - nO(2)) + S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) + +!------------------------------------------------------------------------ +! Compute total energy +!------------------------------------------------------------------------ + +! Kinetic energy + + do ispin=1,nspin + ET(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),T(:,:))) + end do + +! Potential energy + + do ispin=1,nspin + EV(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),V(:,:))) + end do + +! Coulomb energy + + EJ(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,1))) + EJ(2) = trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,2))) + EJ(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),J(:,:,2))) + +! Exchange energy + + do ispin=1,nspin + Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,ispin),K(:,:,ispin))) + end do + +! Correlation energy + + Ec(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,1))) + Ec(2) = trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,2))) + Ec(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),SigC(:,:,2))) + +! Total energy + + EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',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,'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,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,A19,F15.5)')'max(|FPS - SPF|) = ',Conv + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'qsGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'qsGW LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'qsGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'qsGW total energy =',EqsGW + ENuc + write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',sum(Ex(:)) + write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',sum(Ec(:)) + write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*) + +! Dump results for final iteration + + if(Conv < thresh) then + + write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40)') ' Summary ' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',sum(EJ(:)) + sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',EqsGW,' au' + write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A40,1X,F16.10,A3)') ' qsUGW energy: ',EqsGW + ENuc,' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,F13.6)') ' S (exact) :',2d0*S_exact + 1d0 + write(*,'(A40,F13.6)') ' S :',2d0*S + 1d0 + write(*,'(A40,F13.6)') ' (exact) :',S2_exact + write(*,'(A40,F13.6)') ' :',S2 + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A45)') ' Dipole moment (Debye) ' + write(*,'(19X,4A10)') 'X','Y','Z','Tot.' + write(*,'(19X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + + endif + +end subroutine print_qsUGW diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index bb594ef..30d8119 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -35,7 +35,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBAs) + double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) @@ -131,7 +131,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Initialization - nSCF = 0 + nSCF = -1 n_diis = 0 ispin = 1 Conv = 1d0 @@ -147,6 +147,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE do while(Conv > thresh .and. nSCF <= maxSCF) + ! Increment + + nSCF = nSCF + 1 + ! Buid Coulomb matrix call Coulomb_matrix_AO_basis(nBas,P,ERI_AO_basis,J) @@ -223,11 +227,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Print results ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW) - - ! Increment - - nSCF = nSCF + 1 + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EqsGW) enddo !------------------------------------------------------------------------ diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 new file mode 100644 index 0000000..56b1fc9 --- /dev/null +++ b/src/MBPT/qsUGW.f90 @@ -0,0 +1,386 @@ +subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int,dipole_int_aa, & + dipole_int_bb,PHF,cHF,eHF) + +! Perform a quasiparticle self-consistent GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: COHSEX + logical,intent(in) :: SOSEX + 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) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + 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) :: EUHF + 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) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + logical :: doGWPT = .false. + + integer :: nSCF + integer :: nBasSq + integer :: ispin + integer :: is + integer :: n_diis + integer :: nS_aa,nS_bb,nS_sc + double precision :: dipole(ncart) + + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: Conv + double precision :: rcond(nspin) + double precision,external :: trace_matrix + double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: F_diis(:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + double precision,allocatable :: c(:,:,:) + double precision,allocatable :: cp(:,:,:) + double precision,allocatable :: eGW(:,:) + double precision,allocatable :: P(:,:,:) + double precision,allocatable :: F(:,:,:) + double precision,allocatable :: Fp(:,:,:) + double precision,allocatable :: J(:,:,:) + double precision,allocatable :: K(:,:,:) + double precision,allocatable :: SigC(:,:,:) + double precision,allocatable :: SigCp(:,:,:) + double precision,allocatable :: SigCm(:,:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: error(:,:,:) + +! Hello world + + write(*,*) + write(*,*)'*************************************************' + write(*,*)'| Self-consistent unrestricted qsGW calculation |' + write(*,*)'*************************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsUGW !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + +! SOSEX correction + + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if + +! 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 + +! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & + J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),SigCm(nBas,nBas,nspin), & + Z(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(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:,:) = PHF(:,:,:) + eGW(:,:) = eHF(:,:) + c(:,:,:) = cHF(:,:,:) + F_diis(:,:,:) = 0d0 + error_diis(:,:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Buid Coulomb matrix + + do is=1,nspin + call Coulomb_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) + end do + + ! Compute exchange part of the self-energy + + do is=1,nspin + call exchange_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),K(:,:,is)) + end do + + !-------------------------------------------------- + ! AO to MO transformation of two-electron integrals + !-------------------------------------------------- + + ! 4-index transform for (aa|aa) block + + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_aaaa) + + ! 4-index transform for (aa|bb) block + + call AOtoMO_integral_transform(1,1,2,2,nBas,c,ERI_AO,ERI_aabb) + + ! 4-index transform for (bb|bb) block + + call AOtoMO_integral_transform(2,2,2,2,nBas,c,ERI_AO,ERI_bbbb) + + ! Compute linear response + + 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,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 + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + do is=1,nspin + SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) + SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is))) + end do + + do is=1,nspin + call MOtoAO_transform(nBas,S,c(:,:,is),SigCp(:,:,is)) + end do + + ! Solve the quasi-particle equation + + do is=1,nspin + F(:,:,is) = Hc(:,:) + J(:,:,is) + J(:,:,mod(is,2)+1) + K(:,:,is) + SigCp(:,:,is) + end do + + ! Check convergence + + do is=1,nspin + error(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is)) + end do + + if(nSCF > 1) conv = maxval(abs(error(:,:,:))) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + do is=1,nspin + if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) + end do + + ! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + + ! Transform Fock matrix in orthogonal basis + + do is=1,nspin + Fp(:,:,is) = matmul(transpose(X(:,:)),matmul(F(:,:,is),X(:,:))) + end do + + ! Diagonalize Fock matrix to get eigenvectors and eigenvalues + + cp(:,:,:) = Fp(:,:,:) + do is=1,nspin + call diagonalize_matrix(nBas,cp(:,:,is),eGW(:,is)) + end do + + ! Back-transform eigenvectors in non-orthogonal basis + + do is=1,nspin + c(:,:,is) = matmul(X(:,:),cp(:,:,is)) + end do + + ! Compute density matrix + + do is=1,nspin + P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) + end do + + ! Print results + + call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_qsUGW(nBas,nO,S,nSCF,Conv,thresh,eGW,c,P,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,dipole) + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Compute second-order correction of the Hermitization error + +!if(doGWPT) call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm) + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_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, & + S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,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@qsUGW correlation energy (spin-conserved) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW correlation energy (spin-flip) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsUGW total energy =',ENuc + EUHF + 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@qsUGW correlation energy ' + write(*,*) '--------------------------------------------------------------' + write(*,*) + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + + call unrestricted_ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcAC) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy (spin-conserved) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy (spin-flip) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@qsUGW total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end if + +end subroutine qsUGW diff --git a/src/MBPT/unrestricted_self_energy_correlation.f90 b/src/MBPT/unrestricted_self_energy_correlation.f90 new file mode 100644 index 0000000..1fa2f5b --- /dev/null +++ b/src/MBPT/unrestricted_self_energy_correlation.f90 @@ -0,0 +1,94 @@ +subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,j,a,b,p,q,jb + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas,nspin) + +! Initialize + + SigC(:,:,:) = 0d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do jb=1,nSt + eps = e(p,1) - e(i,1) + Omega(jb) + SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(p,1) - e(a,1) - Omega(jb) + SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do jb=1,nSt + eps = e(p,2) - e(i,2) + Omega(jb) + SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(p,2) - e(a,2) - Omega(jb) + SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + +end subroutine unrestricted_self_energy_correlation diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 540d625..5e94a58 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -884,9 +884,22 @@ program QuAcK if(doqsGW) then call cpu_time(start_qsGW) - call qsGW(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,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) + + if(unrestricted) then + + call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int,dipole_int_aa, & + dipole_int_bb,PHF,cHF,eHF) + + else + + call qsGW(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,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF) + + end if + call cpu_time(end_qsGW) t_qsGW = end_qsGW - start_qsGW