diff --git a/examples/molecule.N2 b/examples/molecule.N2
index 93b448f..6905a9c 100644
--- a/examples/molecule.N2
+++ b/examples/molecule.N2
@@ -2,4 +2,4 @@
  2   7     7   0     0
 # Znuc   x            y           z
   N     0.           0.         0.
-  N     0.           0.         3.4
+  N     0.           0.         2.
diff --git a/include/quadrature.h b/include/quadrature.h
index f4e8a13..5058ba9 100644
--- a/include/quadrature.h
+++ b/include/quadrature.h
@@ -15,3 +15,17 @@
              0.054398649583574d0 , 0.046722211728017d0, 0.03805005681419d0 , 0.028567212713429d0, 0.018476894885426d0, &
              0.0080086141288872d0  /)
 
+  integer,parameter      :: nAC = 21
+  double precision, save :: rAC(1:nAC) = &
+          (/ 0.00312391468981d0  , 0.0163865807168d0  , 0.0399503329248d0  , 0.0733183177083d0  , 0.115780018262d0   , &
+             0.166430597901d0    , 0.224190582056d0   , 0.287828939896d0   , 0.355989341599d0   , 0.42721907292d0    , &
+             0.5d0               , 0.57278092708d0    , 0.644010658401d0   , 0.712171060104d0   , 0.775809417944d0   , &
+             0.833569402099d0    , 0.884219981738d0   , 0.926681682292d0   , 0.960049667075d0   , 0.983613419283d0   , &
+             0.99687608531d0       /)
+  double precision, save :: wAc(1:nAC) = &
+          (/ 0.0080086141288872d0, 0.018476894885426d0, 0.028567212713429d0, 0.03805005681419d0 , 0.046722211728017d0, &
+             0.054398649583574d0 , 0.060915708026864d0, 0.066134469316669d0, 0.069943697395537d0, 0.072262201994985d0, &
+             0.07304056682485d0  , 0.072262201994985d0, 0.069943697395537d0, 0.066134469316669d0, 0.060915708026864d0, &
+             0.054398649583574d0 , 0.046722211728017d0, 0.03805005681419d0 , 0.028567212713429d0, 0.018476894885426d0, &
+             0.0080086141288872d0  /)
+
diff --git a/input/basis b/input/basis
index b9ca7b5..4c61da7 100644
--- a/input/basis
+++ b/input/basis
@@ -1,9 +1,16 @@
-1   3
-S   3   1.00
-     38.3600000              0.0238090        
-      5.7700000              0.1548910        
-      1.2400000              0.4699870        
+1   6
+S   4   1.00
+    234.0000000              0.0025870        
+     35.1600000              0.0195330        
+      7.9890000              0.0909980        
+      2.2120000              0.2720500        
 S   1   1.00
-      0.2976000              1.0000000        
+      0.6669000              1.0000000        
+S   1   1.00
+      0.2089000              1.0000000        
 P   1   1.00
-      1.2750000              1.0000000 
+      3.0440000              1.0000000        
+P   1   1.00
+      0.7580000              1.0000000        
+D   1   1.00
+      1.9650000              1.0000000   
diff --git a/input/methods b/input/methods
index 9ac4f33..6379c0b 100644
--- a/input/methods
+++ b/input/methods
@@ -9,8 +9,8 @@
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW 
-  F    F    F
-# G0T0 evGT qsGT 
   T    F    F
+# G0T0 evGT qsGT 
+  F    F    F
 # MCMP2
   F
diff --git a/input/options b/input/options
index 6786f5f..d04d036 100644
--- a/input/options
+++ b/input/options
@@ -1,5 +1,5 @@
 # RHF: maxSCF thresh   DIIS n_diis guess_type ortho_type
-       64    0.000000001  T     5      1          1
+       64    0.0000001  T     5      1          1
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis
@@ -7,8 +7,8 @@
 # CIS/TDHF/BSE: singlet triplet
                 T       T
 # GF: maxSCF thresh  DIIS n_diis renormalization
-      64     0.00001 T    10      3 
+       1     0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta
-      256    0.0000001  T  5      F      F     F   F   F   F   F         0.000
+      256    0.00001    T  5      F      F     T   F   F   F   F         0.000
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/input/weight b/input/weight
index b9ca7b5..4c61da7 100644
--- a/input/weight
+++ b/input/weight
@@ -1,9 +1,16 @@
-1   3
-S   3   1.00
-     38.3600000              0.0238090        
-      5.7700000              0.1548910        
-      1.2400000              0.4699870        
+1   6
+S   4   1.00
+    234.0000000              0.0025870        
+     35.1600000              0.0195330        
+      7.9890000              0.0909980        
+      2.2120000              0.2720500        
 S   1   1.00
-      0.2976000              1.0000000        
+      0.6669000              1.0000000        
+S   1   1.00
+      0.2089000              1.0000000        
 P   1   1.00
-      1.2750000              1.0000000 
+      3.0440000              1.0000000        
+P   1   1.00
+      0.7580000              1.0000000        
+D   1   1.00
+      1.9650000              1.0000000   
diff --git a/scan_H2.sh b/scan_H2.sh
index 57b12ae..cbc6e29 100755
--- a/scan_H2.sh
+++ b/scan_H2.sh
@@ -1,7 +1,7 @@
 #! /bin/bash
 
 MOL="H2"
-BASIS="AVTZ"
+BASIS="VDZ"
 R_START=0.5
 R_END=3.0
 DR=0.01
diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix.f90
index 26f02d7..433fb64 100644
--- a/src/QuAcK/Bethe_Salpeter_A_matrix.f90
+++ b/src/QuAcK/Bethe_Salpeter_A_matrix.f90
@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A_lr)
+subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr)
 
 ! Compute the extra term for Bethe-Salpeter equation for linear response 
 
@@ -8,8 +8,10 @@ subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A_lr)
 ! Input variables
 
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: lambda
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas) 
-  double precision,intent(in)   :: Omega(nS),rho(nBas,nBas,nS)
+  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
   
 ! Local variables
 
@@ -34,8 +36,8 @@ subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A_lr)
             chi = chi + rho(i,j,kc)*rho(a,b,kc)/Omega(kc)
           enddo
 
-          A_lr(ia,jb) = A_lr(ia,jb) - ERI(i,a,j,b) + 4d0*chi
- 
+          A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi
+
         enddo
       enddo
     enddo
diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix.f90
index 903e974..006c518 100644
--- a/src/QuAcK/Bethe_Salpeter_B_matrix.f90
+++ b/src/QuAcK/Bethe_Salpeter_B_matrix.f90
@@ -1,4 +1,4 @@
-subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B_lr)
+subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr)
 
 ! Compute the extra term for Bethe-Salpeter equation for linear response 
 
@@ -8,8 +8,10 @@ subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B_lr)
 ! Input variables
 
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: lambda
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: Omega(nS),rho(nBas,nBas,nS)
+  double precision,intent(in)   :: Omega(nS)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
   
 ! Local variables
 
@@ -34,7 +36,7 @@ subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B_lr)
             chi = chi + rho(i,b,kc)*rho(a,j,kc)/Omega(kc)
           enddo
 
-          B_lr(ia,jb) = B_lr(ia,jb) - ERI(i,a,b,j) + 4d0*chi
+          B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi
 
         enddo
       enddo
diff --git a/src/QuAcK/CIS.f90 b/src/QuAcK/CIS.f90
index 66142f5..2026a36 100644
--- a/src/QuAcK/CIS.f90
+++ b/src/QuAcK/CIS.f90
@@ -18,6 +18,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   logical                       :: dump_matrix = .false.
   logical                       :: dump_trans = .false.
   integer                       :: ispin
+  double precision              :: lambda
   double precision,allocatable  :: A(:,:),Omega(:)
 
 ! Hello world
@@ -28,6 +29,10 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   write(*,*)'************************************************'
   write(*,*)
 
+! Adiabatic connection scaling
+
+  lambda = 1d0
+
 ! Switch on exchange for CIS
 
   dRPA = .false.
@@ -41,7 +46,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   if(singlet_manifold) then
 
     ispin = 1
-    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,eHF,ERI,A)
+    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (singlet state)'
@@ -63,7 +68,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
   if(triplet_manifold) then
 
     ispin = 2
-    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,eHF,ERI,A)
+    call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
  
     if(dump_matrix) then
       print*,'CIS matrix (triplet state)'
diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90
index 5c7a5b8..152fed6 100644
--- a/src/QuAcK/G0W0.f90
+++ b/src/QuAcK/G0W0.f90
@@ -5,6 +5,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
 
   implicit none
   include 'parameters.h'
+  include 'quadrature.h'
 
 ! Input variables
 
@@ -40,6 +41,12 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
   double precision,allocatable  :: rho(:,:,:,:)
   double precision,allocatable  :: rhox(:,:,:,:)
 
+  logical                       :: AC
+  integer                       :: iAC
+  double precision              :: lambda
+  double precision,allocatable  :: EcACBSE(:,:)
+  double precision,allocatable  :: EcAC(:,:)
+
 ! Output variables
 
   double precision              :: eG0W0(nBas)
@@ -75,9 +82,12 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
   allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),  & 
            rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
 
+  AC = .true.
+  allocate(EcACBSE(nAC,nspin),EcAC(nAC,nspin))
+
 ! Compute linear response
 
-  call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & 
+  call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & 
                        rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
 
 ! Compute correlation part of the self-energy 
@@ -118,15 +128,15 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
       ispin = 1
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
                            rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, &
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
-    endif
+    end if
 
    ! Triplet manifold
 
@@ -135,15 +145,15 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
       ispin = 2
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
                            rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, &
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
-    endif
+    end if
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -154,6 +164,85 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
     write(*,*)'-------------------------------------------------------------------------------'
     write(*,*)
 
-  endif
+!   Compute the BSE correlation energy via the adiabatic connection 
+
+    if(AC) then
+
+      write(*,*) '------------------------------------------------------'
+      write(*,*) 'Adiabatic connection version of BSE correlation energy'
+      write(*,*) '------------------------------------------------------'
+      write(*,*) 
+ 
+      if(singlet_manifold) then
+
+        ispin = 1
+        EcACBSE(:,ispin) = 0d0
+
+        write(*,*) '--------------'
+        write(*,*) 'Singlet states'
+        write(*,*) '--------------'
+        write(*,*) 
+
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)'
+        write(*,*) '-----------------------------------------------------------------------------------'
+
+        do iAC=1,nAC
+ 
+          lambda = rAC(iAC)
+
+          call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
+            rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+
+          call Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+
+          write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
+
+        end do
+
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,*)
+
+      end if
+ 
+      if(triplet_manifold) then
+
+        ispin = 2
+        EcACBSE(:,ispin) = 0d0
+
+        write(*,*) '--------------'
+        write(*,*) 'Triplet states'
+        write(*,*) '--------------'
+        write(*,*) 
+
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)'
+        write(*,*) '-----------------------------------------------------------------------------------'
+
+        do iAC=1,nAC
+ 
+          lambda = rAC(iAC)
+
+          call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
+            rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin))
+
+          call Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin))
+
+          write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
+
+        end do
+
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
+        write(*,*) '-----------------------------------------------------------------------------------'
+        write(*,*)
+
+      end if
+ 
+    end if
+
+  end if
 
 end subroutine G0W0
diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90
index e40e9b8..5a2b085 100644
--- a/src/QuAcK/TDHF.f90
+++ b/src/QuAcK/TDHF.f90
@@ -26,6 +26,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
   logical                       :: TDA
   logical                       :: BSE
   integer                       :: ispin
+  double precision              :: lambda
   double precision,allocatable  :: Omega(:,:)
   double precision,allocatable  :: XpY(:,:,:)
 
@@ -44,6 +45,10 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
 
   EcRPA(:) = 0d0
 
+! Adiabatic connection scaling
+
+  lambda = 1d0
+
 ! Switch on exchange for TDHF
 
   dRPA = .false.
@@ -66,7 +71,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
 
     ispin = 1
 
-    call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, &
+    call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
                          EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
     call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
 
@@ -78,7 +83,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
 
     ispin = 2
 
-    call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, &
+    call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
                          EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
     call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
 
diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90
index ad1d0de..7dbb03f 100644
--- a/src/QuAcK/evGW.f90
+++ b/src/QuAcK/evGW.f90
@@ -44,7 +44,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcGM
-  double precision              :: lambda
+  double precision              :: alpha
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: e_diis(:,:)
   double precision,allocatable  :: eGW(:)
@@ -81,7 +81,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
 ! Linear mixing
 
   linear_mixing = .false.
-  lambda = 0.2d0
+  alpha = 0.2d0
 
 ! Memory allocation
 
@@ -165,7 +165,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
 
     if(linear_mixing) then
  
-      eGW(:) = lambda*eGW(:) + (1d0 - lambda)*eOld(:)
+      eGW(:) = alpha*eGW(:) + (1d0 - alpha)*eOld(:)
  
     else
 
@@ -220,11 +220,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
       ispin = 1
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
                            rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
@@ -237,11 +237,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
       ispin = 2
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
                            rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
 
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90
index 5672d55..cda7170 100644
--- a/src/QuAcK/linear_response.f90
+++ b/src/QuAcK/linear_response.f90
@@ -1,4 +1,4 @@
-subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRPA,Omega,XpY)
+subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY)
 
 ! Compute linear response
 
@@ -9,7 +9,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRP
 
   logical,intent(in)            :: dRPA,TDA,BSE
   integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas),rho(nBas,nBas,nS)
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   
 ! Local variables
 
@@ -28,16 +31,16 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRP
 
 ! Build A and B matrices 
 
-  call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A)
-  if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,A)
+  call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
+  if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A)
 
 ! Tamm-Dancoff approximation
 
   B = 0d0
   if(.not. TDA) then
 
-    call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B)
-    if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,ERI,Omega,rho,B)
+    call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
+    if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B)
 
   endif
 
diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/QuAcK/linear_response_A_matrix.f90
index c61ffd9..1a1d391 100644
--- a/src/QuAcK/linear_response_A_matrix.f90
+++ b/src/QuAcK/linear_response_A_matrix.f90
@@ -1,4 +1,4 @@
-subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
+subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_lr)
 
 ! Compute linear response
 
@@ -9,7 +9,9 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
 
   logical,intent(in)            :: dRPA
   integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
@@ -45,8 +47,8 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,e,ERI,A_lr)
           jb = jb + 1
 
           A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                      + (1d0 + delta_spin)*ERI(i,b,a,j) &
-                      - (1d0 - delta_dRPA)*ERI(i,b,j,a)
+                      + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
+                      - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
 
         enddo
       enddo
diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/QuAcK/linear_response_B_matrix.f90
index c91dc0a..5226059 100644
--- a/src/QuAcK/linear_response_B_matrix.f90
+++ b/src/QuAcK/linear_response_B_matrix.f90
@@ -1,4 +1,4 @@
-subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
+subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_lr)
 
 ! Compute linear response
 
@@ -9,6 +9,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
 
   logical,intent(in)            :: dRPA
   integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: lambda
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   
 ! Local variables
@@ -43,8 +44,8 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
         do b=nO+1,nBas-nR
           jb = jb + 1
 
-          B_lr(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) &
-                      - (1d0 - delta_dRPA)*ERI(i,j,b,a)
+          B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
+                      - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
 
         enddo
       enddo
diff --git a/src/QuAcK/mo_guess.f90 b/src/QuAcK/mo_guess.f90
index 9321d13..f7b416b 100644
--- a/src/QuAcK/mo_guess.f90
+++ b/src/QuAcK/mo_guess.f90
@@ -43,6 +43,11 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
 
     call random_number(c)
 
+  else
+
+    print*,'Wrong guess option'
+    stop
+
   endif
 
   P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90
index fca9912..2f76bff 100644
--- a/src/QuAcK/print_excitation.f90
+++ b/src/QuAcK/print_excitation.f90
@@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
 
 ! Input variables
 
-  character*12,intent(in)             :: method
+  character*4,intent(in)             :: method
   integer,intent(in)                 :: ispin,nS
   double precision,intent(in)        :: Omega(nS)
 
diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90
index b12c070..7c34a49 100644
--- a/src/QuAcK/qsGW.f90
+++ b/src/QuAcK/qsGW.f90
@@ -140,7 +140,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
 
     if(.not. GW0 .or. nSCF == 0) then
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
                            rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
 
     endif
@@ -253,11 +253,11 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
       ispin = 1
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
                              rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
      
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & 
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, & 
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
 
@@ -269,11 +269,11 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
       ispin = 2
       EcBSE(ispin) = 0d0
 
-      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
+      call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
                              rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
      
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, &
                            rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
       call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))