scan w

input/dft

@@ -6,7 +6,7 @@
 # GGA          = 2: RB88
 # Hybrid       = 4
 # Hartree-Fock = 666
-  666 HF
+  1 S51  
 # correlation rung: 
 # Hartree      = 0
 # LDA          = 1: RVWN5,RMFL20
@@ -19,7 +19,7 @@
 # Number of states in ensemble (nEns)
   3
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.0 0.0
+  1 0.0 
 # Parameters for CC weight-dependent exchange functional
 0.000000 0.0000000 0.000000
 0.000000 0.0000000 0.0000000

input/int

@@ -1,13 +1,13 @@
 # Debuggin mode?
   F
 # Chemist notation for two-electron integral?
-  F
+  T
 # Exposant of the Slater geminal
-  0.5
+  1.0
 # One-electron integrals: Ov Kin Nuc
                           T  T   T
 # Two-electron integrals: ERI F12 Yuk Erf
-                          T   F   F   T
+                          T   F   F   F
 # Three-electron integrals: Type1 Type2 Type3
                             F     F     F
 # Four-electron integrals: Type1 Type2 Type3

input/options

@@ -9,7 +9,7 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.00367493 3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta        COHSEX SOSEX TDA_W G0W GW0 
-         256    0.00001   T    5     T  0.00367493  F      F     T     F   F  
+         256    0.00001   T    5     T  0.00367493  F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn

scripts/scan_w.sh

@@ -4,8 +4,8 @@ MOL=$1
 BASIS=$2
 
 w_start=0.0
-w_end=1.05
+w_end=1.0
-dw=0.05
+dw=0.1
 
 w2=0.0
 
@@ -15,10 +15,13 @@ CF=$4
 aw1="0.000000 0.0000000 0.000000"
 aw2="0.000000 0.0000000 0.0000000"
 
+DATA=${MOL}_${BASIS}_${XF}_${CF}.dat
+rm $DATA
+touch $DATA
+
 for w1 in $(seq $w_start $dw $w_end)
 do
 ###  w2=${w1}
-  echo "Weights = " $w1 $w2
   echo "# Restricted or unrestricted KS calculation" > input/dft
   echo "  eDFT-UKS" >> input/dft
   echo "# exchange rung:" >> input/dft
@@ -46,6 +49,15 @@ do
   echo ${aw2} >> input/dft
   echo "# GOK-DFT: maxSCF thresh   DIIS n_diis guess_type ortho_type" >> input/dft
   echo "           32    0.00001   T     5      1          1" >> input/dft
-  ./GoXC $MOL $BASIS > ${MOL}_${BASIS}_${XF}_${CF}_${w1}.out
+  OUTPUT=${MOL}_${BASIS}_${XF}_${CF}_${w1}.out
+  ./GoXC $MOL $BASIS > ${OUTPUT}
+  Ew=`grep "Ensemble energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
+  E0=`grep "Individual energy state  1:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
+  E1=`grep "Individual energy state  2:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
+  E2=`grep "Individual energy state  3:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
+  IP=`grep "Ionization Potential"  ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
+  EA=`grep "Electronic Affinity"  ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
+  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA
+  echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA >> ${DATA}
 done
 

src/QuAcK/G0W0.f90

@@ -104,6 +104,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
   call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & 
                        rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
 
+  if(print_W) call print_excitation('RPA@HF      ',ispin,nS,Omega(:,ispin))
+
 ! Compute correlation part of the self-energy 
 
   call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
@@ -143,8 +145,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
 
 ! Dump results
 
-  if(print_W) call print_excitation('RPA@G0W0    ',ispin,nS,Omega(:,ispin))
-
   call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
 
 ! Compute the RPA correlation energy

src/eDFT/UCC_lda_exchange_individual_energy.f90

@@ -86,11 +86,18 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
     r  = max(0d0,rhow(iG))
     rI = max(0d0,rho(iG))
 
-    if(r > threshold .and. rI > threshold) then
+    if(r > threshold) then
 
       e_p  =         Cx*r**(1d0/3d0)
       dedr = 1d0/3d0*Cx*r**(-2d0/3d0)
-      Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
+
+      Ex = Ex - weight(iG)*dedr*r*r
+
+      if(rI > threshold) then
+
+      Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI)
+
+      endif
 
     endif
 

src/eDFT/UVWN5_lda_correlation_individual_energy.f90

@@ -61,7 +61,40 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
    
 !   spin-up contribution
    
-    if(ra > threshold .and. raI > threshold) then
+    if(ra > threshold .or. raI > threshold) then
+   
+      r  = ra
+      rI = raI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+   
+      x_f = x*x + b_f*x + c_f
+   
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+   
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_f = 2d0*x + b_f
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdra_f = drsdra*dxdrs*decdx_f
+
+      Ec(1) = Ec(1) + weight(iG)*(ec_z + decdra_f*r)*rI
+   
+    end if
+
+!   up-down contribution
+   
+    if(ra > threshold .or. raI > threshold) then
    
       r = ra + rb
       rI = raI + rbI
@@ -131,74 +164,39 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
 
 !   spin-down contribution
    
-    if(rb > threshold .and. rbI > threshold) then
+    if(rb > threshold .or. rbI > threshold) then
  
-      r = ra + rb
+      r = rb
-      rI = raI + rbI
+      rI = rbI
 
       rs = (4d0*pi*r/3d0)**(-1d0/3d0)
-      z = (ra - rb)/r
       x = sqrt(rs)
 
-      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
-      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
-   
-      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
-   
-      x_p = x*x + b_p*x + c_p
       x_f = x*x + b_f*x + c_f
-      x_a = x*x + b_a*x + c_a
 
-      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
       xx0_f = x0_f*x0_f + b_f*x0_f + c_f
-      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
 
-      q_p = sqrt(4d0*c_p - b_p*b_p)
       q_f = sqrt(4d0*c_f - b_f*b_f)
-      q_a = sqrt(4d0*c_a - b_a*b_a)
-   
-      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
-                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
 
       ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
                  - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
 
-      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
-                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
-
-      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
-   
-      dzdrb = -(1d0 + z)/r
-      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
-      dfzdrb = dzdrb*dfzdz
-
-      drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
       dxdrs  = 0.5d0/sqrt(rs)
 
-      dxdx_p = 2d0*x + b_p
       dxdx_f = 2d0*x + b_f
-      dxdx_a = 2d0*x + b_a
-
-      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
-                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
 
       decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
                       - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
 
-      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+      decdra_f = drsdra*dxdrs*decdx_f
-                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
 
-      decdrb_p = drsdrb*dxdrs*decdx_p
+      Ec(3) = Ec(3) + weight(iG)*(ec_z + decdra_f*r)*rI  
-      decdrb_f = drsdrb*dxdrs*decdx_f
-      decdrb_a = drsdrb*dxdrs*decdx_a
-
-      decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
-             + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
-
-      Ec(2) = Ec(2) + weight(iG)*(ec_z + decdrb*r)*rI
 
     end if
 
   end do
 
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
 end subroutine UVWN5_lda_correlation_individual_energy

src/eDFT/print_UKS.f90

@@ -17,18 +17,31 @@ subroutine print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
   double precision,intent(in)        :: Ec(nsp)
   double precision,intent(in)        :: Ew
 
+  integer                            :: ispin
   integer                            :: HOMO(nspin)
   integer                            :: LUMO(nspin)
   double precision                   :: Gap(nspin)
 
 ! HOMO and LUMO
 
-  HOMO(:) = nO(:)
 
-  LUMO(:) = HOMO(:) + 1
+  do ispin=1,nspin
 
-  Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1)
+    if(nO(ispin) > 0) then 
-  Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2)
+
+      HOMO(ispin) = nO(ispin)
+      LUMO(ispin) = HOMO(ispin) + 1
+      Gap(ispin)  = eps(LUMO(ispin),ispin) - eps(HOMO(ispin),ispin)
+
+    else
+
+      HOMO(ispin) = 0
+      LUMO(ispin) = 0
+      Gap(ispin)  = 0d0
+
+    end if
+
+  end do
 
 ! Dump results
 

src/eDFT/print_unrestricted_individual_energy.f90

@@ -137,7 +137,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
 !------------------------------------------------------------------------
 
   write(*,'(A60)') '-------------------------------------------------'
-  write(*,'(A60)') ' IP and EA FROM AUXILIARY ENERGIES '
+  write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES '
   write(*,'(A60)') '-------------------------------------------------'
 
     write(*,'(A43,F16.10,A4)') ' Ionization Potential  1 -> 2:',Omaux(2)+OmxcDD(2),' au'

src/eDFT/restricted_lda_correlation_derivative_discontinuity.f90

@@ -28,11 +28,11 @@ subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGr
 
       Ec(:) = 0d0
 
-    case ('RVWN5')
+    case ('VWN5')
 
       Ec(:) = 0d0
 
-    case ('RMFL20')
+    case ('MFL20')
 
       call RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),Ec(:))
 

src/eDFT/restricted_lda_correlation_energy.f90

@@ -29,11 +29,11 @@ subroutine restricted_lda_correlation_energy(DFA,LDA_centered,nEns,wEns,nGrid,we
 
       Ec = 0d0
 
-    case ('RVWN5')
+    case ('VWN5')
 
       call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),Ec)
 
-    case ('RMFL20')
+    case ('MFL20')
 
       call RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns(:),nGrid,weight(:),rho(:),Ec)
 

src/eDFT/restricted_lda_correlation_individual_energy.f90

@@ -26,13 +26,13 @@ subroutine restricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wE
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
-    case ('RVWN5')
+    case ('VWN5')
 
       call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ec)
 
 !   Marut-Fromager-Loos weight-dependent correlation functional
 
-    case ('RMFL20')
+    case ('MFL20')
 
       call RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ec)
 

src/eDFT/restricted_lda_correlation_potential.f90

@@ -31,11 +31,11 @@ subroutine restricted_lda_correlation_potential(DFA,LDA_centered,nEns,wEns,nGrid
 
       Fc(:,:) = 0d0
 
-    case ('RVWN5')
+    case ('VWN5')
 
       call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:))
 
-    case ('RMFL20')
+    case ('MFL20')
 
       call RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:))
 

src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90

@@ -28,13 +28,13 @@ subroutine unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,n
 
 !   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
 
-    case ('UW38')
+    case ('W38')
 
       Ec(:,:) = 0d0
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
-    case ('UVWN5')
+    case ('VWN5')
 
       Ec(:,:) = 0d0
 

src/eDFT/unrestricted_lda_correlation_energy.f90

@@ -28,19 +28,19 @@ subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec
 
       Ec(:) = 0d0
 
-    case ('UW38')
+    case ('W38')
 
       call UW38_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
-    case ('UVWN5')
+    case ('VWN5')
 
       call UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
 
 !   Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101
 
-    case ('UC16')
+    case ('C16')
 
       call UC16_lda_correlation_energy(nGrid,weight,rho,Ec)
 

src/eDFT/unrestricted_lda_correlation_individual_energy.f90

@@ -26,7 +26,7 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
-    case ('UVWN5')
+    case ('VWN5')
 
       call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
 

src/eDFT/unrestricted_lda_correlation_potential.f90

@@ -32,19 +32,19 @@ include 'parameters.h'
 
 !   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
 
-    case ('UW38')
+    case ('W38')
 
       call UW38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
 
 !   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
 
-    case ('UVWN5')
+    case ('VWN5')
 
       call UVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
 
 !   Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101
 
-    case ('UC16')
+    case ('C16')
 
       call UC16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))