From 0be0d5d8c3c2a398f9fb2836bf3edab8dcc9c0a4 Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Sat, 14 Mar 2020 23:00:44 +0100
Subject: [PATCH] revive eDFT code

---
 GoXC                                          |  2 +-
 examples/molecule.F2                          |  2 +-
 examples/molecule.LiH                         |  2 +-
 input/basis                                   | 34 +++++---
 input/dft                                     |  8 +-
 input/methods                                 | 10 +--
 input/molecule                                |  5 +-
 input/molecule.xyz                            |  5 +-
 input/options                                 |  4 +-
 input/weight                                  | 34 +++++---
 scan_H2.sh                                    |  2 +-
 scan_LiH.sh                                   |  4 +-
 src/QuAcK/G0T0.f90                            | 55 +-----------
 src/QuAcK/G0W0.f90                            | 18 ++--
 src/QuAcK/QuAcK.f90                           | 13 +--
 src/QuAcK/evGW.f90                            | 15 +---
 src/QuAcK/excitation_density_Tmatrix.f90      | 30 +++----
 src/QuAcK/linear_response_pp.f90              | 85 ++++++++++---------
 src/QuAcK/renormalization_factor_Tmatrix.f90  | 48 ++++-------
 .../renormalization_factor_Tmatrix_so.f90     | 28 +++---
 src/QuAcK/self_energy_Tmatrix.f90             |  4 +-
 src/QuAcK/self_energy_Tmatrix_diag.f90        | 50 ++++-------
 src/QuAcK/self_energy_Tmatrix_diag_so.f90     | 28 +++---
 src/QuAcK/soG0T0.f90                          |  6 +-
 src/eDFT/Kohn_Sham.f90                        |  9 +-
 src/eDFT/Makefile                             | 10 +--
 src/eDFT/eDFT.f90                             | 15 +++-
 src/eDFT/exchange_energy.f90                  |  2 +-
 src/eDFT/fock_exchange_energy.f90             |  2 +-
 src/eDFT/fock_exchange_potential.f90          |  2 +-
 src/eDFT/hartree_coulomb.f90                  |  2 +-
 31 files changed, 243 insertions(+), 291 deletions(-)

diff --git a/GoXC b/GoXC
index 52f4e57..2a3b543 100755
--- a/GoXC
+++ b/GoXC
@@ -11,8 +11,8 @@ if [ $# = 2 ]
 then
   cp examples/molecule."$1"   input/molecule
   cp examples/basis."$1"."$2" input/basis
+  cp basis/"$2"               input/basis.qcaml
   cp examples/basis."$1"."$2" input/weight
-  ./bin/IntPak
   ./bin/eDFT
 fi
 
diff --git a/examples/molecule.F2 b/examples/molecule.F2
index 8f011d0..ac3ddcd 100644
--- a/examples/molecule.F2
+++ b/examples/molecule.F2
@@ -2,4 +2,4 @@
  2   9     9   0     0
 # Znuc   x            y           z
   F     0.           0.         0.
-  F     0.           0.         2.640
+  F     0.           0.         2.1
diff --git a/examples/molecule.LiH b/examples/molecule.LiH
index 83a41ac..7b1b846 100644
--- a/examples/molecule.LiH
+++ b/examples/molecule.LiH
@@ -2,4 +2,4 @@
  2   2     2   0     0
 # Znuc   x            y           z
   Li    0.           0.         0.
-  H     0.           0.         3.018
+  H     0.           0.         3.5
diff --git a/input/basis b/input/basis
index fb05e68..d92a3df 100644
--- a/input/basis
+++ b/input/basis
@@ -1,18 +1,26 @@
-1   3
+1   8
+S   6
+  1   1264.5857000              0.0019448        
+  2    189.9368100              0.0148351        
+  3     43.1590890              0.0720906        
+  4     12.0986630              0.2371542        
+  5      3.8063232              0.4691987        
+  6      1.2728903              0.3565202        
 S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+  1      3.1964631             -0.1126487
+  2      0.7478133             -0.2295064
+  3      0.2199663              1.1869167
+P   3
+  1      3.1964631              0.0559802
+  2      0.7478133              0.2615506
+  3      0.2199663              0.7939723
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0823099              1.0000000
 P   1
-  1      0.7270000              1.0000000
-2   3
-S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+  1      0.0823099              1.0000000
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0207000              1.0000000
 P   1
-  1      0.7270000              1.0000000
+  1      0.0207000              1.0000000
+D   1
+  1      0.4000000              1.0000000
diff --git a/input/dft b/input/dft
index 90787c0..4044431 100644
--- a/input/dft
+++ b/input/dft
@@ -1,12 +1,12 @@
 # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
-  666 HF
+  1 S51
 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
-  0 HF
+  0 LF19
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
-  1
+  3
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
-  0.00000 0.00000
+  0.00000 0.50000
 # eKS: maxSCF thresh   DIIS n_diis guess_type ortho_type
        64    0.0000001  T     15      1          1
diff --git a/input/methods b/input/methods
index 5b88684..8e05978 100644
--- a/input/methods
+++ b/input/methods
@@ -1,16 +1,16 @@
 # RHF UHF MOM
   T   F   F
 # MP2 MP3 MP2-F12
-  T   F   F
+  F   F   F
 # CCD CCSD CCSD(T) ringCCD ladderCCD
   F   F    F       F       F
 # CIS RPA RPAx ppRPA ADC
-  T   T   T    F     F
+  F   F   F    F     F
 # GF2 GF3
-  F   F
+  T   F
 # G0W0 evGW qsGW 
-  T    F    F
-# G0T0 evGT qsGT 
   F    F    F
+# G0T0 evGT qsGT 
+  T    F    F
 # MCMP2
   F
diff --git a/input/molecule b/input/molecule
index 8b4dc85..6a6f6d1 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,5 +1,4 @@
 # nAt nEla nElb nCore nRyd
- 2   1     1   0     0
+  1   2     2   0     0
 # Znuc   x            y           z
-  H     0.           0.         0.
-  H     0.           0.         2.3
+  Be      0.0          0.0         0.0
diff --git a/input/molecule.xyz b/input/molecule.xyz
index dbee8c5..8023e37 100644
--- a/input/molecule.xyz
+++ b/input/molecule.xyz
@@ -1,4 +1,3 @@
-  2
+  1
 
- H      0.0000000000    0.0000000000    0.0000000000
- H      0.0000000000    0.0000000000    1.2171076727
+ Be     0.0000000000    0.0000000000    0.0000000000
diff --git a/input/options b/input/options
index da3f206..8489957 100644
--- a/input/options
+++ b/input/options
@@ -7,9 +7,9 @@
 # CIS/TDHF/BSE: singlet triplet
                 T       T
 # GF: maxSCF thresh  DIIS n_diis renormalization
-      256    0.00001 T    5      3 
+       256   0.00001 T    5      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta     
-      256    0.00001    T  5     F      F     T   F   F   F   F   0.000   
+       1    0.00001    T  5     F      F     T   F   F   F   T   0.000   
 # ACFDT: AC Kx  XBS
          T  F   T
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
diff --git a/input/weight b/input/weight
index fb05e68..d92a3df 100644
--- a/input/weight
+++ b/input/weight
@@ -1,18 +1,26 @@
-1   3
+1   8
+S   6
+  1   1264.5857000              0.0019448        
+  2    189.9368100              0.0148351        
+  3     43.1590890              0.0720906        
+  4     12.0986630              0.2371542        
+  5      3.8063232              0.4691987        
+  6      1.2728903              0.3565202        
 S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+  1      3.1964631             -0.1126487
+  2      0.7478133             -0.2295064
+  3      0.2199663              1.1869167
+P   3
+  1      3.1964631              0.0559802
+  2      0.7478133              0.2615506
+  3      0.2199663              0.7939723
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0823099              1.0000000
 P   1
-  1      0.7270000              1.0000000
-2   3
-S   3
-  1     13.0100000              0.0196850        
-  2      1.9620000              0.1379770        
-  3      0.4446000              0.4781480        
+  1      0.0823099              1.0000000
 S   1
-  1      0.1220000              1.0000000        
+  1      0.0207000              1.0000000
 P   1
-  1      0.7270000              1.0000000
+  1      0.0207000              1.0000000
+D   1
+  1      0.4000000              1.0000000
diff --git a/scan_H2.sh b/scan_H2.sh
index 0d4103a..2b7a66a 100755
--- a/scan_H2.sh
+++ b/scan_H2.sh
@@ -1,7 +1,7 @@
 #! /bin/bash
 
 MOL="H2"
-BASIS="cc-pvdz"
+BASIS="cc-pvqz"
 R_START=1.0
 R_END=2.4
 DR=0.1
diff --git a/scan_LiH.sh b/scan_LiH.sh
index 08c3478..4643c11 100755
--- a/scan_LiH.sh
+++ b/scan_LiH.sh
@@ -2,9 +2,9 @@
 
 MOL="LiH"
 BASIS="cc-pvqz"
-R_START=2.5
+R_START=3.5
 R_END=3.5
-DR=0.001
+DR=0.1
 
 for R in $(seq $R_START $DR $R_END)
 do
diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90
index 159f3a4..7269a5c 100644
--- a/src/QuAcK/G0T0.f90
+++ b/src/QuAcK/G0T0.f90
@@ -1,6 +1,6 @@
 subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
-! Perform G0W0 calculation with a T-matrix self-energy (G0T0)
+! Perform one-shot calculation with a T-matrix self-energy (G0T0)
 
   implicit none
   include 'parameters.h'
@@ -57,10 +57,10 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
   allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & 
            Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & 
-           rho1s(nBas,nO-nC,nVVs),rho2s(nBas,nV-nR,nOOs), & 
+           rho1s(nBas,nO,nVVs),rho2s(nBas,nV,nOOs), & 
            Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & 
            Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & 
-           rho1t(nBas,nO-nC,nVVt),rho2t(nBas,nV-nR,nOOt), & 
+           rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), & 
            SigT(nBas),Z(nBas),eG0T0(nBas))
 
 !----------------------------------------------
@@ -137,53 +137,4 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
   call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:))
 
-! Perform BSE calculation
-
-! if(BSE) then
-
-!  ! Singlet manifold
-
-!  if(singlet_manifold) then
-
-!     ispin = 1
-!     EcBSE(ispin) = 0d0
-
-!     call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,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,.false.,.false.,BSE,nBas,nC,nO,nV,nR,nS,eG0T0,ERI, &
-!                          rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
-!     call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
-
-!   endif
-
-!  ! Triplet manifold
-
-!  if(triplet_manifold) then
-
-!     ispin = 2
-!     EcBSE(ispin) = 0d0
-
-!     call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,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,eG0T0,ERI, &
-!                          rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
-!     call print_excitation('BSE  ',ispin,nS,Omega(:,ispin))
-
-!   endif
-
-!   write(*,*)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (singlet) =',EcBSE(1)
-!   write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy (triplet) =',EcBSE(2)
-!   write(*,'(2X,A40,F15.6)') 'BSE@G0T0 correlation energy           =',EcBSE(1) + EcBSE(2)
-!   write(*,'(2X,A40,F15.6)') 'BSE@G0T0 total energy                 =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
-!   write(*,*)'-------------------------------------------------------------------------------'
-!   write(*,*)
-
-! endif
-
 end subroutine G0T0
diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90
index 64134ca..0763450 100644
--- a/src/QuAcK/G0W0.f90
+++ b/src/QuAcK/G0W0.f90
@@ -1,4 +1,5 @@
-subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & 
+subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & 
+                singlet_manifold,triplet_manifold,linearize,eta,    & 
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW)
 
 ! Perform G0W0 calculation
@@ -18,6 +19,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
   logical,intent(in)            :: TDA
   logical,intent(in)            :: singlet_manifold
   logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
 
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
@@ -101,16 +103,22 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
 
   eGWlin(:) = eHF(:) + Z(:)*SigC(:)
 
+  if(linearize) then 
+  
+    eGW(:) = eGWlin(:)
+
 ! Find all the roots of the QP equation if necessary
 
 ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
 
-! Find graphical solution of the QP equation
+  else 
+  
+  ! Find graphical solution of the QP equation
 
-! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
-
-  eGW(:) = eGWlin(:)
+  call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
  
+  end if
+
 ! Dump results
 
 ! call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 0edae02..05e6384 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -155,13 +155,13 @@ program QuAcK
 !      = nO*nV
 
   call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
-  allocate(ZNuc(nNuc),rNuc(nNuc,3))
+  allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
 
 ! Read geometry
 
   call read_geometry(nNuc,ZNuc,rNuc,ENuc)
 
-  allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), &
+  allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
            DShell(maxShell,maxK),ExpShell(maxShell,maxK))
 
 !------------------------------------------------------------------------
@@ -547,8 +547,9 @@ program QuAcK
   if(doG0W0) then
     
     call cpu_time(start_G0W0)
-    call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & 
-              eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
+    call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & 
+              singlet_manifold,triplet_manifold,linearize,eta, & 
+              nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
     call cpu_time(end_G0W0)
   
     t_G0W0 = end_G0W0 - start_G0W0
@@ -602,8 +603,8 @@ program QuAcK
   if(doG0T0) then
     
     call cpu_time(start_G0T0)
- !  call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
-    call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
+    call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
+!   call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
     call cpu_time(end_G0T0)
   
     t_G0T0 = end_G0T0 - start_G0T0
diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90
index f1e28fa..c18144c 100644
--- a/src/QuAcK/evGW.f90
+++ b/src/QuAcK/evGW.f90
@@ -1,5 +1,5 @@
 subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & 
-                singlet_manifold,triplet_manifold,linearize,eta,                   & 
+                singlet_manifold,triplet_manifold,eta,                   & 
                 nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
@@ -25,7 +25,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   logical,intent(in)            :: GW0
   logical,intent(in)            :: singlet_manifold
   logical,intent(in)            :: triplet_manifold
-  logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
   double precision,intent(in)   :: eHF(nBas)
@@ -144,17 +143,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
 
     endif
 
-    ! Solve the quasi-particle equation (linearized or not)
+    ! Solve the quasi-particle equation
 
-    if(linearize) then 
-
-      eGW(:) = eHF(:) + Z(:)*SigC(:)
-
-    else
-
-      eGW(:) = eHF(:) + SigC(:)
-
-    endif
+    eGW(:) = eHF(:) + SigC(:)
 
     ! Convergence criteria
 
diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90
index c2ad49b..df0f510 100644
--- a/src/QuAcK/excitation_density_Tmatrix.f90
+++ b/src/QuAcK/excitation_density_Tmatrix.f90
@@ -24,8 +24,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
 ! Output variables
 
-  double precision,intent(out)  :: rho1(nBas,nBas,nVV)
-  double precision,intent(out)  :: rho2(nBas,nBas,nOO)
+  double precision,intent(out)  :: rho1(nBas,nO,nVV)
+  double precision,intent(out)  :: rho2(nBas,nV,nOO)
 
   rho1(:,:,:) = 0d0   
   rho2(:,:,:) = 0d0   
@@ -66,7 +66,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
         end do
       end do
 
-      do a=nO+1,nBas-nR
+      do a=1,nV-nR
         do ij=1,nOO
 
           cd = 0
@@ -74,9 +74,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do d=nO+1,c
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-!                          + ERI(p,a,c,d)*X2(cd,ij)
-                           + (ERI(p,a,c,d) + ERI(p,a,d,c))*X2(cd,ij) &
-                           /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d)))
+!                          + ERI(p,nO+a,c,d)*X2(cd,ij)
+                           + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
+                           /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
 
             end do
           end do
@@ -86,9 +86,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=nC+1,k
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-!                          + ERI(p,a,k,l)*Y2(kl,ij)  
-                           + (ERI(p,a,k,l) + ERI(p,a,l,k))*Y2(kl,ij) &
-                           /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l)))
+!                          + ERI(p,nO+a,k,l)*Y2(kl,ij)  
+                           + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
+                           /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
 
             end do
           end do
@@ -131,7 +131,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
         end do
       end do
 
-      do a=nO+1,nBas-nR
+      do a=1,nV-nR
         do ij=1,nOO
 
           cd = 0
@@ -139,7 +139,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do d=nO+1,c-1
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
+                           + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
             end do
           end do
 
@@ -148,7 +148,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=nC+1,k-1
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
+                           + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
             end do
           end do
 
@@ -190,7 +190,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
         end do
       end do
 
-      do a=nO+1,nBas-nR
+      do a=1,nV-nR
         do ij=1,nOO
 
           cd = 0
@@ -198,7 +198,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
             do d=c+1,nBas-nR
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij)
+                           + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
             end do
           end do
 
@@ -207,7 +207,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k+1,nO
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-                           + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij)
+                           + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
             end do
           end do
 
diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90
index 8321cda..4a42f94 100644
--- a/src/QuAcK/linear_response_pp.f90
+++ b/src/QuAcK/linear_response_pp.f90
@@ -1,6 +1,7 @@
-subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
+subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & 
+                              e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
 
-! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013)
+! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
 
   implicit none
   include 'parameters.h'
@@ -16,6 +17,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   
 ! Local variables
 
+  logical                       :: dump_matrices = .false.
   integer                       :: ab,cd,ij,kl
   integer                       :: p,q,r,s
   double precision              :: trace_matrix
@@ -65,47 +67,52 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   M(    1:nVV    ,nVV+1:nOO+nVV) = -           B(1:nVV,1:nOO)
   M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(B(1:nVV,1:nOO))
 
-  open(unit=42,file='B.dat')
-  open(unit=43,file='C.dat')
-  open(unit=44,file='D.dat')
-  open(unit=45,file='ERI.dat')
-  open(unit=46,file='eps.dat')
- 
-  do ab=1,nVV
-    do ij=1,nOO
-      if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij)
-    end do
-  end do
- 
-  do ab=1,nVV
-    do cd=1,nVV
-      if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd)
-    end do
-  end do
- 
-  do ij=1,nOO
-    do kl=1,nOO
-      if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl)
-    end do
-  end do
- 
-  do p=1,nBas
-    write(46,*) p,e(p)
-    do q=1,nBas
-      do r=1,nBas
-        do s=1,nBas
-      if(abs(ERI(p,q,r,s)) > 1d-15) write(45,*) p,q,r,s,ERI(p,q,r,s)
+! Dump ppRPA matrices 
+
+  if(dump_matrices) then 
+
+    open(unit=42,file='B.dat')
+    open(unit=43,file='C.dat')
+    open(unit=44,file='D.dat')
+    open(unit=45,file='ERI.dat')
+    open(unit=46,file='eps.dat')
+  
+    do ab=1,nVV
+      do ij=1,nOO
+        if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij)
       end do
     end do
+  
+    do ab=1,nVV
+      do cd=1,nVV
+        if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd)
+      end do
     end do
-  end do
-
-  close(42)
-  close(43)
-  close(44)
-  close(45)
-  close(46)
+  
+    do ij=1,nOO
+      do kl=1,nOO
+        if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl)
+      end do
+    end do
+  
+    do p=1,nBas
+      write(46,*) p,e(p)
+      do q=1,nBas
+        do r=1,nBas
+          do s=1,nBas
+        if(abs(ERI(p,q,r,s)) > 1d-15) write(45,*) p,q,r,s,ERI(p,q,r,s)
+        end do
+      end do
+      end do
+    end do
+ 
+    close(42)
+    close(43)
+    close(44)
+    close(45)
+    close(46)
 
+  end if
 
 ! print*, 'pp-RPA matrix'
 ! call matout(nOO+nVV,nOO+nVV,M(:,:))
diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90
index 76a0c5a..14ff12a 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix.f90
@@ -15,9 +15,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
   integer,intent(in)            :: nVVs,nVVt
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1s(nVVs),Omega1t(nVVt)
-  double precision,intent(in)   :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
+  double precision,intent(in)   :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
   double precision,intent(in)   :: Omega2s(nOOs),Omega2t(nOOt)
-  double precision,intent(in)   :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt)
+  double precision,intent(in)   :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
 
 ! Local variables
 
@@ -40,13 +40,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=nO+1,c
-          cd = cd + 1
-          eps = e(p) + e(i) - Omega1s(cd)
-          Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
-        enddo
+      do cd=1,nVVs
+        eps = e(p) + e(i) - Omega1s(cd)
+        Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
       enddo
     enddo
   enddo
@@ -54,14 +50,10 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
 ! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=nC+1,k
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2s(kl)
-          Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
-        enddo
+    do a=1,nV-nR
+      do kl=1,nOOs
+        eps = e(p) + e(nO+a) - Omega2s(kl)
+        Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
       enddo
     enddo
   enddo
@@ -74,13 +66,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=nO+1,c-1
-          cd = cd + 1
-          eps = e(p) + e(i) - Omega1t(cd)
-          Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2
-        enddo
+      do cd=1,nVVt
+        eps = e(p) + e(i) - Omega1t(cd)
+        Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2
       enddo
     enddo
   enddo
@@ -88,14 +76,10 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
 ! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=nC+1,k-1
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2t(kl)
-          Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2
-        enddo
+    do a=1,nV-nR
+      do kl=1,nOOt
+        eps = e(p) + e(nO+a) - Omega2t(kl)
+        Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/renormalization_factor_Tmatrix_so.f90 b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
index e06ec78..aa3692a 100644
--- a/src/QuAcK/renormalization_factor_Tmatrix_so.f90
+++ b/src/QuAcK/renormalization_factor_Tmatrix_so.f90
@@ -13,9 +13,9 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1(nVV)
-  double precision,intent(in)   :: rho1(nBas,nBas,nVV)
+  double precision,intent(in)   :: rho1(nBas,nO,nVV)
   double precision,intent(in)   :: Omega2(nOO)
-  double precision,intent(in)   :: rho2(nBas,nBas,nOO)
+  double precision,intent(in)   :: rho2(nBas,nV,nOO)
 
 ! Local variables
 
@@ -38,13 +38,13 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=c+1,nBas-nR
-          cd = cd + 1
+      do cd=1,nVV
+!     do c=nO+1,nBas-nR
+!       do d=c+1,nBas-nR
+!         cd = cd + 1
           eps = e(p) + e(i) - Omega1(cd)
           Z(p) = Z(p) - rho1(p,i,cd)**2/eps**2
-        enddo
+!       enddo
       enddo
     enddo
   enddo
@@ -52,14 +52,14 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
 ! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=k+1,nO
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2(kl)
+    do a=1,nV-nR
+      do kl=1,nOO
+!     do k=nC+1,nO
+!       do l=k+1,nO
+!         kl = kl + 1
+          eps = e(p) + e(nO+a) - Omega2(kl)
           Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
-        enddo
+!       enddo
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/self_energy_Tmatrix.f90 b/src/QuAcK/self_energy_Tmatrix.f90
index 25cc8e6..361b17e 100644
--- a/src/QuAcK/self_energy_Tmatrix.f90
+++ b/src/QuAcK/self_energy_Tmatrix.f90
@@ -17,9 +17,9 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,Omega2,rho1
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1(nVV)
-  double precision,intent(in)   :: rho1(nBas,nBas,nVV)
+  double precision,intent(in)   :: rho1(nBas,nO,nVV)
   double precision,intent(in)   :: Omega2(nOO)
-  double precision,intent(in)   :: rho2(nBas,nBas,nOO)
+  double precision,intent(in)   :: rho2(nBas,nV,nOO)
 
 ! Local variables
 
diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90
index c01784b..b793985 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag.f90
@@ -19,9 +19,9 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
   integer,intent(in)            :: nVVs,nVVt
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1s(nVVs),Omega1t(nVVt)
-  double precision,intent(in)   :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt)
+  double precision,intent(in)   :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
   double precision,intent(in)   :: Omega2s(nOOs),Omega2t(nOOt)
-  double precision,intent(in)   :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt)
+  double precision,intent(in)   :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
 
 ! Local variables
 
@@ -44,13 +44,9 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=nO+1,c
-          cd = cd + 1
-          eps = e(p) + e(i) - Omega1s(cd)
-          SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps
-        enddo
+      do cd=1,nVVs
+        eps = e(p) + e(i) - Omega1s(cd)
+        SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps
       enddo
     enddo
   enddo
@@ -58,14 +54,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
   ! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=nC+1,k
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2s(kl)
-          SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps
-        enddo
+    do a=1,nV-nR
+      do kl=1,nOOs
+        eps = e(p) + e(nO+a) - Omega2s(kl)
+        SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps
       enddo
     enddo
   enddo
@@ -78,28 +70,20 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=nO+1,c-1
-          cd = cd + 1
-          eps = e(p) + e(i) - Omega1t(cd)
-          SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps
-        enddo
+      do cd=1,nVVt
+        eps = e(p) + e(i) - Omega1t(cd)
+        SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps
       enddo
     enddo
   enddo
 
-  ! Virtual part of the T-matrix self-energy
+! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=nC+1,k-1
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2t(kl)
-          SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps
-        enddo
+    do a=1,nV-nR
+      do kl=1,nOOt
+        eps = e(p) + e(nO+a) - Omega2t(kl)
+        SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/self_energy_Tmatrix_diag_so.f90 b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
index 3f52f24..ab96502 100644
--- a/src/QuAcK/self_energy_Tmatrix_diag_so.f90
+++ b/src/QuAcK/self_energy_Tmatrix_diag_so.f90
@@ -17,9 +17,9 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: e(nBas)
   double precision,intent(in)   :: Omega1(nVV)
-  double precision,intent(in)   :: rho1(nBas,nBas,nVV)
+  double precision,intent(in)   :: rho1(nBas,nO,nVV)
   double precision,intent(in)   :: Omega2(nOO)
-  double precision,intent(in)   :: rho2(nBas,nBas,nOO)
+  double precision,intent(in)   :: rho2(nBas,nV,nOO)
 
 ! Local variables
 
@@ -42,13 +42,13 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
 
   do p=nC+1,nBas-nR
     do i=nC+1,nO
-      cd = 0
-      do c=nO+1,nBas-nR
-        do d=c+1,nBas-nR
-          cd = cd + 1
+      do cd=1,nVV
+!     do c=nO+1,nBas-nR
+!       do d=c+1,nBas-nR
+!         cd = cd + 1
           eps = e(p) + e(i) - Omega1(cd)
           SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
-        enddo
+!       enddo
       enddo
     enddo
   enddo
@@ -56,14 +56,14 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
 ! Virtual part of the T-matrix self-energy
 
   do p=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      kl = 0
-      do k=nC+1,nO
-        do l=k+1,nO
-          kl = kl + 1
-          eps = e(p) + e(a) - Omega2(kl)
+    do a=1,nV-nR
+      do kl=1,nOO
+!     do k=nC+1,nO
+!       do l=k+1,nO
+!         kl = kl + 1
+          eps = e(p) + e(nO+a) - Omega2(kl)
           SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
-        enddo
+!       enddo
       enddo
     enddo
   enddo
diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90
index 01a1602..4e829af 100644
--- a/src/QuAcK/soG0T0.f90
+++ b/src/QuAcK/soG0T0.f90
@@ -66,9 +66,9 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Memory allocation
 
-  allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV),         & 
-           Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO),         & 
-           rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), & 
+  allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV),     & 
+           Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO),     & 
+           rho1(nBas2,nO2,nVV),rho2(nBas2,nV2,nOO), & 
            eG0T0(nBas2),SigT(nBas2),Z(nBas2))
 
 !----------------------------------------------
diff --git a/src/eDFT/Kohn_Sham.f90 b/src/eDFT/Kohn_Sham.f90
index 4010a5c..bed604f 100644
--- a/src/eDFT/Kohn_Sham.f90
+++ b/src/eDFT/Kohn_Sham.f90
@@ -1,5 +1,5 @@
-subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis, &
-                     guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
+subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis,guess_type, &
+                     nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew)
 
 ! Perform unrestricted Kohn-Sham calculation for ensembles
 
@@ -224,7 +224,7 @@ subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thr
 
       do ispin=1,nspin
         do iEns=1,nEns
-          call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
+!         call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
         end do
       end do
 
@@ -271,7 +271,8 @@ subroutine Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thr
 
     n_diis = min(n_diis+1,max_diis)
     do ispin=1,nspin
-      call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
+      call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, & 
+                              err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
     end do
 
 !   Reset DIIS if required
diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile
index 2c3316c..448a2e1 100644
--- a/src/eDFT/Makefile
+++ b/src/eDFT/Makefile
@@ -14,19 +14,19 @@ ifeq ($(PROFIL),1)
   FC += -p -fno-inline
 endif
 
-LIBS = ~/Dropbox/quack/lib/*.a
-#LIBS = -lblas -llapack
+LIBS = ../../lib/*.a
+LIBS += -lblas -llapack
 
 SRCF90 = $(wildcard *.f90)
 
-SRC = $(wildcard *.f)
+SRC = $(wildcard *.F)
 
-OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
+OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.F,$(ODIR)/%.o,$(SRC))
 
 $(ODIR)/%.o: %.f90 
 	$(FC) -c -o $@ $< $(FFLAGS) 
 
-$(ODIR)/%.o: %.f
+$(ODIR)/%.o: %.F
 	$(FC) -c -o $@ $< $(FFLAGS) 
 
 $(BDIR)/eDFT: $(OBJ)
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index 9ebae91..510bd40 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -31,6 +31,7 @@ program eDFT
   double precision,allocatable  :: dAO(:,:,:)
 
   double precision              :: start_KS,end_KS,t_KS
+  double precision              :: start_int,end_int,t_int
 
   integer                       :: nEns
   double precision,allocatable  :: wEns(:)
@@ -96,7 +97,17 @@ program eDFT
 
 ! Read integrals
 
-  call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI)  
+  call cpu_time(start_int)
+
+  call system('./GoQCaml')
+  call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI)
+
+  call cpu_time(end_int)
+
+  t_int = end_int - start_int
+  write(*,*)
+  write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds'
+  write(*,*)
 
 ! Orthogonalization X = S^(-1/2)
 
@@ -122,7 +133,7 @@ program eDFT
 !------------------------------------------------------------------------
 
     call cpu_time(start_KS)
-    call Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(1:nEns),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
+    call Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
                    nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS)
     call cpu_time(end_KS)
 
diff --git a/src/eDFT/exchange_energy.f90 b/src/eDFT/exchange_energy.f90
index 8e63441..34b2497 100644
--- a/src/eDFT/exchange_energy.f90
+++ b/src/eDFT/exchange_energy.f90
@@ -17,7 +17,7 @@ subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,
   double precision,intent(in)   :: P(nBas,nBas)
   double precision,intent(in)   :: FxHF(nBas,nBas)
   double precision,intent(in)   :: rho(nGrid)
-  double precision,intent(in)   :: drho(3,nGrid)
+  double precision,intent(in)   :: drho(ncart,nGrid)
 
 ! Local variables
 
diff --git a/src/eDFT/fock_exchange_energy.f90 b/src/eDFT/fock_exchange_energy.f90
index f019c04..17e7d49 100644
--- a/src/eDFT/fock_exchange_energy.f90
+++ b/src/eDFT/fock_exchange_energy.f90
@@ -20,6 +20,6 @@ subroutine fock_exchange_energy(nBas,P,Fx,Ex)
 
 ! Compute HF exchange energy
 
-  Ex = trace_matrix(nBas,matmul(P,Fx))
+  Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx))
 
 end subroutine fock_exchange_energy
diff --git a/src/eDFT/fock_exchange_potential.f90 b/src/eDFT/fock_exchange_potential.f90
index 59af3b1..483f08a 100644
--- a/src/eDFT/fock_exchange_potential.f90
+++ b/src/eDFT/fock_exchange_potential.f90
@@ -25,7 +25,7 @@ subroutine fock_exchange_potential(nBas,P,ERI,Fx)
     do si=1,nBas
       do la=1,nBas
         do mu=1,nBas
-          Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,si,la,nu)
+          Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
         enddo
       enddo
     enddo
diff --git a/src/eDFT/hartree_coulomb.f90 b/src/eDFT/hartree_coulomb.f90
index 468d086..42f83d7 100644
--- a/src/eDFT/hartree_coulomb.f90
+++ b/src/eDFT/hartree_coulomb.f90
@@ -23,7 +23,7 @@ subroutine hartree_coulomb(nBas,P,ERI,J)
     do nu=1,nBas
       do la=1,nBas
         do si=1,nBas
-          J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,nu,la,si)
+          J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
         enddo
       enddo
     enddo