diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90
index 3eecd1d..9640a48 100644
--- a/src/HF/cRHF.f90
+++ b/src/HF/cRHF.f90
@@ -1,5 +1,5 @@
 subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & 
-                nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P)
+                nBas,nO,S,T,V,ERI,dipole_int,X,ERHF,eHF,c,P)
 
 ! Perform complex restricted Hartree-Fock calculation
 
@@ -25,7 +25,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
   double precision,intent(in)   :: S(nBas,nBas)
   double precision,intent(in)   :: T(nBas,nBas)
   double precision,intent(in)   :: V(nBas,nBas)
-  double precision,intent(in)   :: Hc(nBas,nBas) 
   double precision,intent(in)   :: X(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
@@ -44,23 +43,25 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
   double precision              :: Conv
   double precision              :: rcond
   double precision,external     :: trace_matrix
-  double precision,allocatable  :: err(:,:)
-  double precision,allocatable  :: err_diis(:,:)
-  double precision,allocatable  :: F_diis(:,:)
-  double precision,allocatable  :: J(:,:)
-  double precision,allocatable  :: K(:,:)
-  double precision,allocatable  :: cp(:,:)
-  double precision,allocatable  :: F(:,:)
-  double precision,allocatable  :: Fp(:,:)
 
-  complex*16,allocatable        :: W(:,:)
+  double precision              :: eta
+  double precision,allocatable  :: W(:,:)
+  complex*16,allocatable        :: Hc(:,:)
+  complex*16,allocatable        :: J(:,:)
+  complex*16,allocatable        :: K(:,:)
+  complex*16,allocatable        :: cp(:,:)
+  complex*16,allocatable        :: F(:,:)
+  complex*16,allocatable        :: Fp(:,:)
+  complex*16,allocatable        :: err(:,:)
+  complex*16,allocatable        :: err_diis(:,:)
+  complex*16,allocatable        :: F_diis(:,:)
 
 ! Output variables
 
-  double precision,intent(out)  :: ERHF
-  double precision,intent(out)  :: eHF(nBas)
-  double precision,intent(inout):: c(nBas,nBas)
-  double precision,intent(out)  :: P(nBas,nBas)
+  complex*16,intent(out)        :: ERHF
+  complex*16,intent(out)        :: eHF(nBas)
+  complex*16,intent(out)        :: c(nBas,nBas)
+  complex*16,intent(out)        :: P(nBas,nBas)
 
 ! Hello world
 
@@ -73,16 +74,23 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 ! Useful quantities
 
   nBasSq = nBas*nBas
+  eta = 0.01d0
 
 ! Memory allocation
 
   allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
            Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis),     & 
-           W(nBas,nBas))
+           Hc(nBas,nBas),W(nBas,nBas))
 
 ! Read CAP integrals from file
 
-! call read_CAP_integrals()
+  call read_CAP_integrals(nBas,W)
+
+  W(:,:) = eta*W(:,:)
+
+! Define core Hamiltonian
+
+  Hc(:,:) = dcmplx(T+V,W)
 
 ! Guess coefficients and density matrix
 
@@ -205,10 +213,10 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 
   if(dotest) then
  
-    call dump_test_value('R','RHF energy',ERHF)
-    call dump_test_value('R','RHF HOMO energy',eHF(nO))
-    call dump_test_value('R','RHF LUMO energy',eHF(nO+1))
-    call dump_test_value('R','RHF dipole moment',norm2(dipole))
+!   call dump_test_value('R','RHF energy',ERHF)
+!   call dump_test_value('R','RHF HOMO energy',eHF(nO))
+!   call dump_test_value('R','RHF LUMO energy',eHF(nO+1))
+!   call dump_test_value('R','RHF dipole moment',norm2(dipole))
 
   end if
 
diff --git a/src/utils/elements.f90 b/src/utils/elements.f90
index a8a35a2..38473ef 100644
--- a/src/utils/elements.f90
+++ b/src/utils/elements.f90
@@ -175,5 +175,4 @@ function element_covalent_radius(zatom)
  element_covalent_radius = element_covalent_radius*pmtoau
 
 
-end function element_covalent_radius
-
+end function 
diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90
index 37d2b21..176862a 100644
--- a/src/utils/orthogonalization_matrix.f90
+++ b/src/utils/orthogonalization_matrix.f90
@@ -116,4 +116,4 @@ subroutine orthogonalization_matrix(nBas,S,X)
 
   end if
 
-end subroutine orthogonalization_matrix
+end subroutine 
diff --git a/src/utils/read_CAP_integrals.f90 b/src/utils/read_CAP_integrals.f90
new file mode 100644
index 0000000..388b064
--- /dev/null
+++ b/src/utils/read_CAP_integrals.f90
@@ -0,0 +1,47 @@
+subroutine read_CAP_integrals(nBas,W)
+
+! Read one- and two-electron integrals from files
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+
+! Local variables
+
+  logical                       :: debug
+  integer                       :: mu,nu
+  double precision              :: wx,wy,wz
+
+! Output variables
+
+  double precision,intent(out)  :: W(nBas,nBas)
+
+! Open file with integrals
+
+  debug = .false.
+
+  open(unit=31,file='int/CAP.dat')
+
+! Read CAP integrals
+
+  W(:,:) = 0d0
+  do 
+    read(31,*,end=31) mu,nu,wx,wy,wz
+    W(mu,nu) = wx + wy + wz
+    W(nu,mu) = wx + wy + wz
+  end do
+  31 close(unit=31)
+
+! Print results
+  if(debug) then
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'CAP integrals'
+    write(*,'(A28)') '----------------------'
+    call matout(nBas,nBas,W)
+    write(*,*)
+  end if
+
+end subroutine 
diff --git a/src/utils/read_auxiliary_basis.f90 b/src/utils/read_auxiliary_basis.f90
index 2819348..e44cadb 100644
--- a/src/utils/read_auxiliary_basis.f90
+++ b/src/utils/read_auxiliary_basis.f90
@@ -173,4 +173,4 @@ subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, &
 
   close(unit=3)
 
-end subroutine read_auxiliary_basis
+end subroutine 
diff --git a/src/utils/read_basis.f90 b/src/utils/read_basis.f90
index c2efad3..43cdc36 100644
--- a/src/utils/read_basis.f90
+++ b/src/utils/read_basis.f90
@@ -182,4 +182,4 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh
 
   nV(:) = nBas - nO(:)
 
-end subroutine read_basis
+end subroutine 
diff --git a/src/utils/read_dipole_integrals.f90 b/src/utils/read_dipole_integrals.f90
index 553b46c..095bb18 100644
--- a/src/utils/read_dipole_integrals.f90
+++ b/src/utils/read_dipole_integrals.f90
@@ -68,4 +68,4 @@ subroutine read_dipole_integrals(nBas,R)
     call matout(nBas,nBas,R(:,:,3))
   end if
 
-end subroutine read_dipole_integrals
+end subroutine 
diff --git a/src/utils/read_geometry.f90 b/src/utils/read_geometry.f90
index 5ee3916..f1803cf 100644
--- a/src/utils/read_geometry.f90
+++ b/src/utils/read_geometry.f90
@@ -73,4 +73,4 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
   write(*,'(A28)') '------------------'
   write(*,*)
 
-end subroutine read_geometry
+end subroutine 
diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90
index 3b2baeb..9e215e4 100644
--- a/src/utils/read_integrals.f90
+++ b/src/utils/read_integrals.f90
@@ -130,4 +130,4 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
     write(*,*)
   end if
 
-end subroutine read_integrals
+end subroutine 
diff --git a/src/utils/read_molecule.f90 b/src/utils/read_molecule.f90
index 4c018f3..8bcec0b 100644
--- a/src/utils/read_molecule.f90
+++ b/src/utils/read_molecule.f90
@@ -59,4 +59,4 @@ subroutine read_molecule(nNuc,nO,nC,nR)
 
   close(unit=1)
 
-end subroutine read_molecule
+end subroutine 
diff --git a/src/utils/spatial_to_spin_ERI.f90 b/src/utils/spatial_to_spin_ERI.f90
index 8b82d85..7ccff58 100644
--- a/src/utils/spatial_to_spin_ERI.f90
+++ b/src/utils/spatial_to_spin_ERI.f90
@@ -32,4 +32,4 @@ subroutine spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
       end do
     end do
 
-end subroutine spatial_to_spin_ERI
+end subroutine 
diff --git a/src/utils/spatial_to_spin_MO_energy.f90 b/src/utils/spatial_to_spin_MO_energy.f90
index 158d6ad..2cfd4e8 100644
--- a/src/utils/spatial_to_spin_MO_energy.f90
+++ b/src/utils/spatial_to_spin_MO_energy.f90
@@ -26,4 +26,4 @@ subroutine spatial_to_spin_MO_energy(nBas,e,nBas2,se)
 ! print*,'MO energies in spinorbital basis'
 ! call matout(nBas2,1,se)
 
-end subroutine spatial_to_spin_MO_energy
+end subroutine 
diff --git a/src/utils/spatial_to_spin_fock.f90 b/src/utils/spatial_to_spin_fock.f90
index db694ab..c37759c 100644
--- a/src/utils/spatial_to_spin_fock.f90
+++ b/src/utils/spatial_to_spin_fock.f90
@@ -26,4 +26,4 @@ subroutine spatial_to_spin_fock(nBas,F,nBas2,sF)
       end do
     end do
 
-end subroutine spatial_to_spin_fock
+end subroutine