new version of SOSEX (need serious debugging)

include/parameters.h

@@ -1,7 +1,7 @@
   integer,parameter              :: ncart = 3
   integer,parameter              :: nspin = 2
   integer,parameter              :: nsp = 3
-  integer,parameter              :: maxEns = 3
+  integer,parameter              :: maxEns = 4
   integer,parameter              :: maxShell = 512
   integer,parameter              :: maxL = 7
   integer,parameter              :: n1eInt = 3

input/dft

@@ -6,31 +6,34 @@
 # GGA          = 2: B88,G96,PBE
 # MGGA         = 3: 
 # Hybrid       = 4: HF,B3LYP,PBE
-  4 HF
+  1 S51
 # correlation rung: 
 # Hartree      = 0: H
 # LDA          = 1: PW92,VWN3,VWN5,eVWN5
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
 # Hybrid       = 4: HF,B3LYP,PBE
- 0 H
+ 1 VWN5
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
-  1
+  4
 # occupation numbers 
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
- 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+
+ 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
- 0.0 0.0
+ 0.25 0.25 0.25
 # N-centered?
- F   
+ T   
 # Parameters for CC weight-dependent exchange functional
 0.0 0.0 0.0
 0.0 0.0 0.0

input/methods

@@ -13,9 +13,9 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F     F      F    F
 # G0W0* evGW* qsGW* ufG0W0 ufGW
-  F    F    F       F      F
+  T    F    F       F      F
 # G0T0 evGT qsGT 
-  T    F    F
+  F    F    F
 # MCMP2
   F
 # * unrestricted version available

input/options

@@ -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     T  0.0  F      F     F     F   F  
+         256   0.00001   T    5     T  0.0  F      T     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    T    T    F
+        F    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

src/MBPT/G0W0.f90

@@ -1,4 +1,4 @@
-subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 
+subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & 
                 dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,            & 
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
 
@@ -14,7 +14,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
   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
@@ -76,13 +75,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
 
   EcRPA = 0d0
 
-! SOSEX correction
-
-  if(SOSEX) then 
-    write(*,*) 'SOSEX correction activated but BUG!'
-    stop
-  end if
-
 ! COHSEX approximation
 
   if(COHSEX) then 

src/MBPT/G0W0_SOSEX.f90

@@ -0,0 +1,196 @@
+subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, & 
+                      nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eSOSEX)
+
+! Perform the SOSEX extension of G0W0
+
+  implicit none
+  include 'parameters.h'
+  include 'quadrature.h'
+
+! Input variables
+
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: exchange_kernel
+  logical,intent(in)            :: doXBS
+  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)            :: singlet
+  logical,intent(in)            :: triplet
+  double precision,intent(in)   :: eta
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: ERI_AO(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_MO(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: Vxc(nBas)
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: cHF(nBas,nBas)
+  double precision,intent(in)   :: PHF(nBas,nBas)
+
+! Local variables
+
+  logical                       :: print_W = .true.
+  integer                       :: ispin
+  double precision              :: EcRPA(nspin)
+  double precision              :: EcBSE(nspin)
+  double precision              :: EcAC(nspin)
+  double precision              :: EcGM
+  double precision,allocatable  :: SigX(:)
+  double precision,allocatable  :: SigC(:)
+  double precision,allocatable  :: Z(:)
+  double precision,allocatable  :: OmRPA(:,:)
+  double precision,allocatable  :: XpY_RPA(:,:,:)
+  double precision,allocatable  :: XmY_RPA(:,:,:)
+  double precision,allocatable  :: rho_RPA(:,:,:,:)
+
+! Output variables
+
+  double precision              :: eSOSEX(nBas)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'|          One-shot SOSEX calculation          |'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Initialization
+
+  EcRPA = 0d0
+
+! 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
+
+  allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS,nspin),XpY_RPA(nS,nS,nspin),XmY_RPA(nS,nS,nspin),rho_RPA(nBas,nBas,nS,nspin))
+
+!-------------------!
+! Compute screening !
+!-------------------!
+
+  do ispin=1,nspin
+    call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & 
+                         OmRPA(:,ispin),rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
+    if(print_W) call print_excitation('RPA@HF      ',ispin,nS,OmRPA)
+  end do
+
+!--------------------------!
+! Compute spectral weights !
+!--------------------------!
+
+  call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
+
+!------------------------!
+! Compute GW self-energy !
+!------------------------!
+
+  call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX)
+  call self_energy_correlation_SOSEX_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC)
+
+!--------------------------------!
+! Compute renormalization factor !
+!--------------------------------!
+
+  call renormalization_factor_SOSEX(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z)
+
+!-----------------------------------!
+! Solve the quasi-particle equation !
+!-----------------------------------!
+
+  eSOSEX(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:))
+
+! Compute the RPA correlation energy
+
+  do ispin=1,nspin
+    call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO,OmRPA(:,ispin), & 
+                         rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin))
+  end do
+
+!--------------!
+! Dump results !
+!--------------!
+
+  call print_SOSEX(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eSOSEX,EcRPA,EcGM)
+
+! Deallocate memory
+
+  deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA)
+
+! Perform BSE calculation
+
+  if(BSE) then
+
+    call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eSOSEX,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@SOSEX correlation energy (singlet) =',EcBSE(1)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX correlation energy (triplet) =',EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX correlation energy           =',EcBSE(1) + EcBSE(2)
+    write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX 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@SOSEX 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_MO,eHF,eSOSEX,EcAC)
+
+      write(*,*)
+      write(*,*)'-------------------------------------------------------------------------------'
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy (singlet) =',EcAC(1)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy (triplet) =',EcAC(2)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy           =',EcAC(1) + EcAC(2)
+      write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2)
+      write(*,*)'-------------------------------------------------------------------------------'
+      write(*,*)
+
+    end if
+
+  end if
+
+end subroutine G0W0_SOSEX

src/MBPT/evGW.f90

@@ -1,4 +1,4 @@
-subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,  & 
+subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,  & 
                 G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,        &
                 ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0)
 
@@ -18,7 +18,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   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
@@ -78,13 +77,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   write(*,*)'************************************************'
   write(*,*)
 
-! SOSEX correction
-
-  if(SOSEX) then 
-    write(*,*) 'SOSEX correction activated but BUG!'
-    stop
-  end if
-
 ! COHSEX approximation
 
   if(COHSEX) then 

src/MBPT/excitation_density_SOSEX.f90

@@ -1,31 +1,63 @@
-subroutine excitation_density_SOSEX(nBas,nC,nO,nR,nS,G,XpY,rho)
+subroutine excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY,rho)
 
-! Compute excitation densities
+! Compute excitation densities for SOSEX
 
   implicit none
+  include 'parameters.h'
 
 ! Input variables
 
-  integer,intent(in)            :: nBas,nC,nO,nR,nS
+  integer,intent(in)            :: nBas
-  double precision,intent(in)   :: G(nBas,nBas,nBas,nBas),XpY(nS,nS)
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: XpY(nS,nS,nspin)
 
 ! Local variables
 
-  integer                       :: ia,jb,x,y,j,b
+  integer                       :: ispin
+  integer                       :: p,q
+  integer                       :: i,a,j,b
+  integer                       :: ia,jb
 
 ! Output variables
 
-  double precision,intent(out)  :: rho(nBas,nBas,nS)
+  double precision,intent(out)  :: rho(nBas,nBas,nS,nspin)
+
+  rho(:,:,:,:) = 0d0
+
+! Singlet part
+
+  ispin = 1
 
-  rho(:,:,:) = 0d0
   do ia=1,nS
-    do x=nC+1,nBas-nR
+    do p=nC+1,nBas-nR
-      do y=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
         jb = 0
         do j=nC+1,nO
           do b=nO+1,nBas-nR
             jb = jb + 1
-            rho(x,y,ia) = rho(x,y,ia) + G(x,y,b,j)*XpY(ia,jb)
+            rho(p,q,ia,ispin) = rho(p,q,ia,ispin) + ERI(p,j,q,b)*XpY(ia,jb,ispin)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+! Triplet part
+
+  ispin = 2
+
+  do ia=1,nS
+    do p=nC+1,nBas-nR
+      do q=nC+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            rho(p,q,ia,ispin) = rho(p,q,ia,ispin) + (ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb,ispin)
           enddo
         enddo
       enddo

src/MBPT/print_SOSEX.f90

@@ -0,0 +1,55 @@
+subroutine print_SOSEX(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eSOSEX,EcRPA,EcGM)
+
+! Print one-electron energies and other stuff for SOSEX
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas,nO
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ERHF
+  double precision,intent(in)        :: EcRPA
+  double precision,intent(in)        :: EcGM
+  double precision,intent(in)        :: eHF(nBas)
+  double precision,intent(in)        :: SigC(nBas)
+  double precision,intent(in)        :: Z(nBas)
+  double precision,intent(in)        :: eSOSEX(nBas)
+
+  integer                            :: p,HOMO,LUMO
+  double precision                   :: Gap
+
+! HOMO and LUMO
+
+  HOMO = nO
+  LUMO = HOMO + 1
+  Gap = eSOSEX(LUMO)-eSOSEX(HOMO)
+
+! Dump results
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)'  One-shot SOSEX calculation  '
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+            '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
+  write(*,*)'-------------------------------------------------------------------------------'
+
+  do p=1,nBas
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eSOSEX(p)*HaToeV,'|'
+  enddo
+
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') 'SOSEX HOMO      energy:',eSOSEX(HOMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'SOSEX LUMO      energy:',eSOSEX(LUMO)*HaToeV,' eV'
+  write(*,'(2X,A30,F15.6,A3)') 'SOSEX HOMO-LUMO gap   :',Gap*HaToeV,' eV'
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A30,F15.6,A3)') 'RPA@SOSEX total energy      :',ENuc + ERHF + EcRPA,' au'
+  write(*,'(2X,A30,F15.6,A3)') 'RPA@SOSEX correlation energy:',EcRPA,' au'
+  write(*,'(2X,A30,F15.6,A3)') 'GM@SOSEX  total energy      :',ENuc + ERHF + EcGM,' au'
+  write(*,'(2X,A30,F15.6,A3)') 'GM@SOSEX  correlation energy:',EcGM,' au'
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+end subroutine print_SOSEX
+
+

src/MBPT/qsGW.f90

@@ -1,4 +1,4 @@
-subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,          & 
+subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,          & 
                 G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
                 S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
@@ -16,7 +16,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   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
@@ -113,13 +112,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
   nBasSq = nBas*nBas
 
-! SOSEX correction
-
-  if(SOSEX) then 
-    write(*,*) 'SOSEX correction activated but BUG!'
-    stop
-  end if
-
 ! COHSEX approximation
 
   if(COHSEX) then 

src/MBPT/qsUGW.f90

@@ -1,4 +1,4 @@
-subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,    & 
+subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,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_AO,dipole_int_aa, &
                 dipole_int_bb,PHF,cHF,eHF)
@@ -17,7 +17,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
   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
@@ -121,13 +120,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS
 
   nBasSq = nBas*nBas
 
-! SOSEX correction
-
-  if(SOSEX) then 
-    write(*,*) 'SOSEX correction activated but BUG!'
-    stop
-  end if
-
 ! COHSEX approximation
 
   if(COHSEX) then 

src/MBPT/renormalization_factor_SOSEX.f90

@@ -0,0 +1,73 @@
+subroutine renormalization_factor_SOSEX(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z)
+
+! Compute renormalization factor for the SOSEX version of GW
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Omega(nS,nspin)
+  double precision,intent(in)   :: rho(nBas,nBas,nS,nspin)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: i,j,a,b,p,jb
+  double precision              :: eps
+
+! Output variables
+
+  double precision,intent(out)  :: Z(nBas)
+
+! Initialize
+
+  Z(:)  = 0d0
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nBas-nR
+          jb = jb + 1
+          do ispin=1,nspin
+            eps = e(p) - e(i) + Omega(jb,ispin) 
+            Z(p) = Z(p)  - 2d0*rho(p,i,jb,ispin)**2*(eps/(eps**2 + eta**2))**2
+          end do
+        end do
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nBas-nR
+          jb = jb + 1
+          do ispin=1,nspin
+            eps = e(p) - e(a) - Omega(jb,ispin) 
+            Z(p) = Z(p)  - 2d0*rho(p,a,jb,ispin)**2*(eps/(eps**2 + eta**2))**2
+          end do
+        end do
+      end do
+    end do
+  end do
+
+! Compute renormalization factor from derivative of SigC
+ 
+  Z(:) = 1d0/(1d0 - Z(:))
+
+end subroutine renormalization_factor_SOSEX

src/MBPT/self_energy_correlation_SOSEX_diag.f90

@@ -0,0 +1,80 @@
+subroutine self_energy_correlation_SOSEX_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,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
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: Omega(nS,nspin)
+  double precision,intent(in)   :: rho(nBas,nBas,nS,nspin)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: i,a,p,q,jb
+  double precision              :: eps
+
+! Output variables
+
+  double precision,intent(out)  :: SigC(nBas)
+  double precision,intent(out)  :: EcGM
+
+! Initialize 
+
+  SigC(:) = 0d0
+
+!-----------------------------
+! SOSEX self-energy
+!-----------------------------
+ 
+  ! Occupied part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do jb=1,nS
+        do ispin=1,nspin
+          eps = e(p) - e(i) + Omega(jb,ispin)
+          SigC(p) = SigC(p) + 2d0*rho(p,i,jb,ispin)**2*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC+1,nBas-nR
+    do a=nO+1,nBas-nR
+      do jb=1,nS
+        do ispin=1,nspin
+          eps = e(p) - e(a) - Omega(jb,ispin)
+          SigC(p) = SigC(p) + 2d0*rho(p,a,jb,ispin)**2*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+  ! GM correlation energy
+
+  EcGM = 0d0
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      do jb=1,nS
+        do ispin=1,nspin
+          eps = e(a) - e(i) + Omega(jb,ispin)
+          EcGM = EcGM - 4d0*rho(a,i,jb,ispin)*rho(a,i,jb,ispin)*eps/(eps**2 + eta**2)
+        end do
+      end do
+    end do
+  end do
+
+end subroutine self_energy_correlation_SOSEX_diag

src/QuAcK/QuAcK.f90

@@ -953,8 +953,10 @@ program QuAcK
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0)
     else
 
-      call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
+      call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
                 linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
+!     call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
+!                     eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
 
 !     call ufBSE(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
 
@@ -984,7 +986,7 @@ program QuAcK
 
     else
 
-      call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,       &
+      call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, &
                 BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
     end if
@@ -1006,14 +1008,14 @@ program QuAcK
 
     if(unrestricted) then 
     
-      call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
+      call qsUGW(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,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_AO,      &
                  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,         &
+      call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,               &
                 BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & 
                 nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
@@ -1220,7 +1222,7 @@ program QuAcK
 !   ! Long-range G0W0 calculation
 
 !   call cpu_time(start_G0W0)
-!   call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
+!   call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,       & 
 !             dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & 
 !             nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0)
 !   call cpu_time(end_G0W0)

src/eDFT/read_options_dft.f90

@@ -22,7 +22,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
 
   character(len=8),intent(out)  :: method
   integer,intent(out)           :: x_rung,c_rung
-  character(len=12),intent(out) :: x_DFA, c_DFA
+  character(len=12),intent(out) :: x_DFA,c_DFA
   integer,intent(out)           :: SGn
   integer,intent(out)           :: nEns
   logical,intent(out)           :: doNcentered
@@ -39,11 +39,11 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
 
 ! Default values
 
-  method  = 'GOK-RKS'
+  method  = 'eDFT-UKS'
   x_rung  = 1
   c_rung  = 1
-  x_DFA   = 'RS51'
+  x_DFA   = 'S51'
-  c_DFA   = 'RVWN5'
+  c_DFA   = 'VWN5'
   SGn     = 0
   wEns(:) = 0d0
 
@@ -134,15 +134,16 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
 
   if(answer == 'T') doNcentered = .true.
 
+  wEns(1) = 1d0
+  do iEns=2,nEns
+    wEns(1) = wEns(1) - wEns(iEns)
+  end do
+
   if (doNcentered) then
 
-    wEns(1) = 1d0 - wEns(2) - wEns(3)
+    do iEns=2,nEns
-    wEns(2) = (nEl(1)/nEl(2))*wEns(2)
+      wEns(iEns) = (nEl(1)/nEl(iEns))*wEns(iEns)
-    wEns(3) = (nEl(1)/nEl(3))*wEns(3)
+    end do
-
-  else
-
-    wEns(1) = 1d0 - wEns(2) - wEns(3)
 
   end if
 

src/eDFT/unrestricted_auxiliary_energy.f90

@@ -46,10 +46,6 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux)
       end do
     end do
 
-!    if (doNcentered .NE. 0) then
-!      Eaux(:,iEns) = Eaux(:,iEns)*(nEl(iEns)/nEl(1))   
-!    end if 
-
   end do
 
 end subroutine unrestricted_auxiliary_energy

src/eDFT/unrestricted_correlation_individual_energy.f90

@@ -26,7 +26,6 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
 
   double precision              :: EcLDA(nsp)
   double precision              :: EcGGA(nsp)
-  double precision              :: aC
 
 ! Output variables