graph sol for UGW

input/methods

@@ -9,11 +9,11 @@
 # CIS CID CISD
   F   F   F    
 # RPA RPAx ppRPA 
-  T   F    F    
+  F   F    F    
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0 evGW qsGW 
-  F    F    F
+  T    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -5,11 +5,11 @@
 # CC: maxSCF thresh   DIIS n_diis
       64     0.0000001  T    5
 # spin: singlet triplet spin_conserved spin_flip TDA
-        T       T       T              T          F
+        T       T       T              F          F
 # 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.001   F      F     F     F   F  
+         256    0.00001   T    5     F  0.0         F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn

src/QuAcK/SigmaC.f90

@@ -1,4 +1,4 @@
-double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
+double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -7,7 +7,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
 
 ! Input variables
 
-  integer,intent(in)            :: x
+  integer,intent(in)            :: p
   double precision,intent(in)   :: w
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas
@@ -22,7 +22,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
 
 ! Local variables
 
-  integer                       :: i,j,a,b,p,jb
+  integer                       :: i,a,jb
   double precision              :: eps
 
 ! Initialize 
@@ -32,26 +32,18 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
 ! Occupied part of the correlation self-energy
 
   do i=nC+1,nO
-    jb = 0
+    do jb=1,nS
-    do j=nC+1,nO
-      do b=nO+1,nBas-nR
-        jb = jb + 1
       eps = w - e(i) + Omega(jb)
-        SigmaC = SigmaC + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2)
+      SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2)
-      enddo
     enddo
   enddo
 
 ! Virtual part of the correlation self-energy
 
   do a=nO+1,nBas-nR
-    jb = 0
+    do jb=1,nS
-    do j=nC+1,nO
-      do b=nO+1,nBas-nR
-        jb = jb + 1
       eps = w - e(a) - Omega(jb)
-        SigmaC = SigmaC + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2)
+      SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2)
-      enddo
     enddo
   enddo
 

src/QuAcK/UG0W0.f90

@@ -45,6 +45,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 ! Local variables
 
   logical                       :: print_W = .true.
+  integer                       :: is
   integer                       :: ispin
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
@@ -126,13 +127,13 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 ! Compute self-energy !
 !---------------------!
 
-  call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,SigC)
+  call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,SigC)
 
 !--------------------------------!
 ! Compute renormalization factor !
 !--------------------------------!
 
-  call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,Z)
+  call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,Z)
 
 !-----------------------------------!
 ! Solve the quasi-particle equation !
@@ -151,10 +152,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   
   ! Find graphical solution of the QP equation
 
-! do is=1,nspin
+  do is=1,nspin
-!   call QP_graph(nBas,nC(:,is),nO(:,is),nV(:,is),nR(:,is),nS(:,is),eta,eHF(:,is),Omega(:,is), & 
+    call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),Omega_sc, & 
-!                 rho(:,:,:,ispin),eGWlin(:,is),eGW(:,is))
+                               rho_sc,eGWlin(:,is),eGW(:,is))
-! end do
+  end do
  
   end if
 

src/QuAcK/USigmaC.f90

@@ -0,0 +1,48 @@
+double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
+
+! Compute diagonal of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  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)
+  double precision,intent(in)   :: rho(nBas,nBas,nS,nspin)
+
+! Local variables
+
+  integer                       :: i,a,jb
+  double precision              :: eps
+
+! Initialize 
+
+  USigmaC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do jb=1,nS
+      eps = w - e(i) + Omega(jb)
+      USigmaC = uSigmaC + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2)
+    end do
+  end do
+
+  do a=nO+1,nBas-nR
+    do jb=1,nS
+      eps = w - e(a) - Omega(jb)
+      USigmaC = USigmaC + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2)
+    end do
+  end do
+
+end function USigmaC

src/QuAcK/dUSigmaC.f90

@@ -0,0 +1,50 @@
+double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho)
+
+! Compute the derivative of the correlation part of the self-energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: p
+  double precision,intent(in)   :: w
+  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)
+  double precision,intent(in)   :: rho(nBas,nBas,nS,nspin)
+
+! Local variables
+
+  integer                       :: i,a,jb
+  double precision              :: eps
+
+! Initialize 
+
+  dUSigmaC = 0d0
+
+! Occupied part of the correlation self-energy
+
+  do i=nC+1,nO
+    do jb=1,nS
+      eps = w - e(i) + Omega(jb)
+      dUSigmaC = dUSigmaC + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2
+    end do
+  end do
+
+! Virtual part of the correlation self-energy
+
+  do a=nO+1,nBas-nR
+    do jb=1,nS
+      eps = w - e(a) - Omega(jb)
+      dUSigmaC = dUSigmaC + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2
+    end do
+  end do
+
+end function dUSigmaC

src/QuAcK/unrestricted_Bethe_Salpeter.f90

@@ -76,6 +76,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
                                       eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin),        & 
                                       OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
 
+!   call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc)
+
     call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
                                          XpY_RPA_sc,rho_RPA_sc)
 

src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90

@@ -1,6 +1,5 @@
 subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
-                                                ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,     & 
+                                                ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr)
-                                                Omega,rho,A_lr)
 
 ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
 
@@ -40,7 +39,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
 ! Build part A of the BSE matrix !
 !--------------------------------!
 
-  ! alpha-alpha block
+  ! aaaa block
 
   ia = 0
   do i=nC(1)+1,nO(1)
@@ -55,7 +54,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
             chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps &
-                      + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
+                      + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
           enddo
 
           A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi
@@ -65,7 +64,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
     enddo
   enddo
 
-  ! alpha-beta block
+  ! aabb block
 
   ia = 0
   do i=nC(1)+1,nO(1)
@@ -90,7 +89,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
     enddo
   enddo
 
-  ! beta-alpha block
+  ! bbaa block
 
   ia = 0
   do i=nC(2)+1,nO(2)
@@ -104,8 +103,8 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           chi = 0d0
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps &
+            chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps &
-                      + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
+                      + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps
           enddo
 
           A_lr(nSa+ia,jb) = A_lr(nSa+ia,jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi
@@ -115,7 +114,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
     enddo
   enddo
 
-  ! beta-beta block
+  ! bbbb block
 
   ia = 0
   do i=nC(2)+1,nO(2)
@@ -129,7 +128,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           chi = 0d0
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps &
+            chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps &
                       + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps
           enddo
 

src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90

@@ -51,7 +51,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
             chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & 
-                      + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
+                      + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
           enddo
 
           B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi
@@ -100,8 +100,8 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           chi = 0d0
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & 
+            chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & 
-                      + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
+                      + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps
           enddo
 
           B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi
@@ -125,7 +125,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
           chi = 0d0
           do kc=1,nSt
             eps = Omega(kc)**2 + eta**2
-            chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & 
+            chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & 
                       + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps
           enddo
 

src/QuAcK/unrestricted_QP_graph.f90

@@ -0,0 +1,83 @@
+subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
+
+! Compute the graphical solution of the QP equation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  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)   :: eta
+  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS,nspin)
+
+  double precision,intent(in)   :: eGWlin(nBas)
+
+! Local variables
+
+  integer                       :: p
+  integer                       :: nIt
+  integer,parameter             :: maxIt = 10
+  double precision,parameter    :: thresh = 1d-6
+  double precision,external     :: USigmaC,dUSigmaC
+  double precision              :: sig,dsig
+  double precision              :: f,df
+  double precision              :: w
+
+! Output variables
+
+  double precision,intent(out)  :: eGW(nBas)
+
+! Run Newton's algorithm to find the root
+ 
+  do p=nC+1,nBas-nR
+
+    write(*,*) '-----------------'
+    write(*,'(A10,I3)') 'Orbital ',p
+    write(*,*) '-----------------'
+
+    w = eGWlin(p)
+    nIt = 0
+    f = 1d0
+    write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f
+    
+    do while (abs(f) > thresh .and. nIt < maxIt)
+    
+      nIt = nIt + 1
+
+      sig  =  USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
+      dsig = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
+      f  = w - eHF(p) - sig
+      df = 1d0 - dsig
+    
+      w = w - f/df
+
+      write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig
+    
+    
+    end do
+ 
+    if(nIt == maxIt) then 
+
+      write(*,*) 'Newton root search has not converged!'
+      eGW(p) = eGWlin(p)
+
+    else
+
+      eGW(p) = w
+
+      write(*,'(A32,F16.10)')   'Quasiparticle energy (eV)   ',eGW(p)*HaToeV
+      write(*,*)
+
+    end if
+          
+  end do
+
+end subroutine unrestricted_QP_graph

src/QuAcK/unrestricted_renormalization_factor.f90

@@ -1,4 +1,4 @@
-subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,Z)
+subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z)
 
 ! Compute the renormalization factor in the unrestricted formalism
 
@@ -13,8 +13,6 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,
   integer,intent(in)            :: nO(nspin)
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
-  integer,intent(in)            :: nSa
-  integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
   double precision,intent(in)   :: e(nBas,nspin)
   double precision,intent(in)   :: Omega(nSt)
@@ -89,5 +87,4 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,
 
   Z(:,:) = 1d0/(1d0 + Z(:,:))
 
-
 end subroutine unrestricted_renormalization_factor

src/QuAcK/unrestricted_self_energy_correlation_diag.f90

@@ -1,4 +1,4 @@
-subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,SigC)
+subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -13,8 +13,6 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS
   integer,intent(in)            :: nO(nspin)
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
-  integer,intent(in)            :: nSa
-  integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
   double precision,intent(in)   :: e(nBas,nspin)
   double precision,intent(in)   :: Omega(nSt)