diff --git a/include/parameters.h b/include/parameters.h
index cf69187..2109908 100644
--- a/include/parameters.h
+++ b/include/parameters.h
@@ -1,4 +1,7 @@
+  integer,parameter              :: ncart = 3
   integer,parameter              :: nspin = 2
+  integer,parameter              :: nsp   = 3
+  integer,parameter              :: maxEns = 10
   integer,parameter              :: maxShell = 50
   integer,parameter              :: n1eInt = 3
   integer,parameter              :: n2eInt = 4
@@ -6,6 +9,7 @@
   integer,parameter              :: n4eInt = 3
   integer,parameter              :: maxK = 20
 
+  double precision,parameter    :: threshold = 1d-15
   double precision,parameter    :: pi = acos(-1d0)
   double precision,parameter    :: HaToeV = 27.21138602d0
   double precision,parameter    :: pmtoau = 0.0188973d0
diff --git a/input/int b/input/int
index ac39e96..86aa017 100644
--- a/input/int
+++ b/input/int
@@ -2,6 +2,8 @@
   F
 # Chemist notation for two-electron integral?
   T
+# Exposant of the Slater geminal
+  1.0
 # One-electron integrals: Ov Kin Nuc
                           T  T   T
 # Two-electron integrals: ERI F12 Yuk Erf
diff --git a/input/molecule b/input/molecule
index e9384c4..375e1fb 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,4 @@
-# nAt nEl nCore nRyd
-  1   2   0     0
+# nAt nEl nEla nElb
+  1   2   1     1
 # Znuc   x            y           z
-  2.0      0.0          0.0         0.0
+  He      0.0          0.0         0.0
diff --git a/src/IntPak/CalcNBasis.f90 b/src/IntPak/CalcNBasis.f90
deleted file mode 100644
index 49fb2a9..0000000
--- a/src/IntPak/CalcNBasis.f90
+++ /dev/null
@@ -1,28 +0,0 @@
-subroutine CalcNBasis(nShell,atot,NBasis)
-
-  implicit none
-
-! Input variables
-
-  integer,intent(in)   :: nShell
-  integer,intent(in)   :: atot(nShell)
-
-! Local variables
-
-  integer              :: iShell
-
-! Output variables
-
-  integer,intent(out)  :: NBasis
-
-  NBasis = 0
-  do iShell=1,nShell
-    NBasis = NBasis + (atot(iShell)*atot(iShell) + 3*atot(iShell) + 2)/2
-  enddo
-
-  write(*,'(A28)') '------------------'
-  write(*,'(A28,1X,I16)') 'Number of basis functions',NBasis
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-end subroutine CalcNBasis
diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90
index 2f6d5eb..1f51110 100644
--- a/src/IntPak/IntPak.f90
+++ b/src/IntPak/IntPak.f90
@@ -18,11 +18,13 @@ program IntPak
   logical                       :: do3eInt(n3eInt)
   logical                       :: do4eInt(n4eInt)
 
-  integer                       :: NAtoms,NBasis,iType
+  integer                       :: nNuc,nBas,iType
+  integer                       :: nEl,nO,nV
   double precision              :: ExpS
+  double precision              :: ENuc
   integer                       :: KG
   double precision,allocatable  :: DG(:),ExpG(:)
-  double precision,allocatable  :: ZNuc(:),XYZAtoms(:,:)
+  double precision,allocatable  :: ZNuc(:),rNuc(:,:)
 
   integer                       :: nShell
   integer,allocatable           :: TotAngMomShell(:),KShell(:)
@@ -52,7 +54,7 @@ program IntPak
 
 ! Read options for integral calculations
 
-  call read_options(debug,chemist_notation,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt)
+  call read_options(debug,chemist_notation,ExpS,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt)
 
 ! Which integrals do you want?
 
@@ -61,26 +63,29 @@ program IntPak
 ! Read input information
 !------------------------------------------------------------------------
 
-  call ReadNAtoms(NAtoms)
+! Read number of atoms, number of electrons of the system
+! nO   = number of occupied orbitals
+! nV   = number of virtual orbitals (see below)
+! nBas = number of basis functions (see below)
+!      = nO + nV
 
-  allocate(ZNuc(1:NAtoms),XYZAtoms(1:NAtoms,1:3))
+  call read_molecule(nNuc,nEl,nO)
 
-  call ReadGeometry(NAtoms,ZNuc,XYZAtoms)
+  allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3))
+
+! Read geometry
+
+  call read_geometry(nNuc,ZNuc,rNuc,ENuc)
 
   allocate(CenterShell(1:maxShell,1:3),TotAngMomShell(1:maxShell),KShell(1:maxShell), &
            DShell(1:maxShell,1:maxK),ExpShell(1:maxShell,1:maxK))
 
-  call ReadBasis(NAtoms,XYZAtoms,nShell,CenterShell, &
-                     TotAngMomShell,KShell,DShell,ExpShell)
-
-  call CalcNBasis(nShell,TotAngMomShell,NBasis)
-
-  call ReadGeminal(ExpS)
+  call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
 
 !------------------------------------------------------------------------
 ! Memory allocation
 !------------------------------------------------------------------------
-  allocate(S(1:NBasis,1:NBasis))
+  allocate(S(1:nBas,1:nBas))
 
 !------------------------------------------------------------------------
 ! Compute one-electron overlap integrals
@@ -90,7 +95,7 @@ program IntPak
     iType = 1
 
     call cpu_time(start_1eInt(iType))
-    call ComputeOv(debug,NBasis,nShell,                  &
+    call ComputeOv(debug,nBas,nShell,                  &
                     CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
                     np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType),S)
     call cpu_time(end_1eInt(iType))
@@ -148,7 +153,7 @@ program IntPak
     call cpu_time(start_1eInt(iType))
     call ComputeNuc(debug,nShell,                        &
                     CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
-                    NAtoms,ZNuc,XYZAtoms,                              &
+                    nNuc,ZNuc,rNuc,                              &
                     np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType))
     call cpu_time(end_1eInt(iType))
 
diff --git a/src/IntPak/ReadBasis.f90 b/src/IntPak/ReadBasis.f90
deleted file mode 100644
index 9a1db2a..0000000
--- a/src/IntPak/ReadBasis.f90
+++ /dev/null
@@ -1,176 +0,0 @@
-subroutine ReadBasis(NAtoms,XYZAtoms,nShell,CenterShell, &
-                     TotAngMomShell,KShell,DShell,ExpShell)
-
-! Read basis set information
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: NAtoms
-  double precision,intent(in)   :: XYZAtoms(NAtoms,3)
-
-! Local variables
-
-  integer                       :: nShAt,iAt
-  integer                       :: i,j,k
-  character                     :: shelltype
-
-! Output variables
-
-  integer,intent(out)           :: nShell
-  double precision,intent(out)  :: CenterShell(maxShell,3)
-  integer,intent(out)           :: TotAngMomShell(maxShell),KShell(maxShell)
-  double precision,intent(out)  :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
-
-!------------------------------------------------------------------------
-! Primary basis set information
-!------------------------------------------------------------------------
-
-! Open file with basis set specification
-
-  open(unit=2,file='input/basis')
-
-! Read basis information
-
-  write(*,'(A28)') 'Gaussian basis set'
-  write(*,'(A28)') '------------------'
-
-  nShell = 0
-  do i=1,NAtoms
-    read(2,*) iAt,nShAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
-    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
-    write(*,'(A28)') '------------------'
-
-!   Basis function centers
-
-    do j=1,nShAt
-      nShell = nShell + 1
-      do k=1,3
-        CenterShell(nShell,k) = XYZAtoms(iAt,k)
-      enddo
-
-! Shell type and contraction degree
-
-      read(2,*) shelltype,KShell(nShell)
-        if(shelltype == "S") then
-          TotAngMomShell(nShell) = 0
-          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "P") then
-          TotAngMomShell(nShell) = 1
-          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "D") then
-          TotAngMomShell(nShell) = 2
-          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "F") then
-          TotAngMomShell(nShell) = 3
-          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "G") then
-          TotAngMomShell(nShell) = 4
-          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "H") then
-          TotAngMomShell(nShell) = 5
-          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "I") then
-          TotAngMomShell(nShell) = 6
-          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
-        endif
-
-! Read exponents and contraction coefficients
-
-        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
-        do k=1,Kshell(nShell)
-          read(2,*) ExpShell(nShell,k),DShell(nShell,k)
-          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
-        enddo
-    enddo
-    write(*,'(A28)') '------------------'
-  enddo
-
-! Total number of shells
-
-  write(*,'(A28,1X,I16)') 'Number of shells in OBS',nShell
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with basis set specification
-
-  close(unit=2)
-
-!------------------------------------------------------------------------
-! Auxiliary basis set information
-!------------------------------------------------------------------------
-
-! Open file with auxilairy basis specification
-
-  open(unit=3,file='input/auxbasis')
-
-! Read basis information
-
-  write(*,'(A28)') 'Auxiliary basis set'
-  write(*,'(A28)') '-------------------'
-
-  do i=1,NAtoms
-    read(3,*) iAt,nShAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
-    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
-    write(*,'(A28)') '------------------'
-
-!   Basis function centers
-
-    do j=1,nShAt
-      nShell = nShell + 1
-      do k=1,3
-        CenterShell(nShell,k) = XYZAtoms(iAt,k)
-      enddo
-
-! Shell type and contraction degree
-
-      read(3,*) shelltype,KShell(nShell)
-        if(shelltype == "S") then
-          TotAngMomShell(nShell) = 0
-          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "P") then
-          TotAngMomShell(nShell) = 1
-          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "D") then
-          TotAngMomShell(nShell) = 2
-          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "F") then
-          TotAngMomShell(nShell) = 3
-          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "G") then
-          TotAngMomShell(nShell) = 4
-          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "H") then
-          TotAngMomShell(nShell) = 5
-          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "I") then
-          TotAngMomShell(nShell) = 6
-          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
-        endif
-
-! Read exponents and contraction coefficients
-
-        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
-        do k=1,Kshell(nShell)
-          read(3,*) ExpShell(nShell,k),DShell(nShell,k)
-          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
-        enddo
-    enddo
-    write(*,'(A28)') '------------------'
-  enddo
-
-! Total number of shells
-
-  write(*,'(A28,1X,I16)') 'Number of shells in ABS',nShell
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with basis set specification
-
-  close(unit=3)
-
-end subroutine ReadBasis
diff --git a/src/IntPak/ReadNAtoms.f90 b/src/IntPak/ReadNAtoms.f90
deleted file mode 100644
index 9f78096..0000000
--- a/src/IntPak/ReadNAtoms.f90
+++ /dev/null
@@ -1,20 +0,0 @@
-subroutine ReadNAtoms(NAtoms)
-
-! Read number of atoms
-
-  implicit none
-
-! Input variables
-  integer,intent(out)           :: NAtoms
-
-! Open file with geometry specification
-  open(unit=1,file='input/molecule')
-
-! Read number of atoms
-  read(1,*)
-  read(1,*) NAtoms
-
-! Close file with geometry specification
-  close(unit=1)
-
-end subroutine ReadNAtoms
diff --git a/src/MCQC/MCQC.f90 b/src/MCQC/MCQC.f90
index 6d94385..42c8f3a 100644
--- a/src/MCQC/MCQC.f90
+++ b/src/MCQC/MCQC.f90
@@ -65,13 +65,13 @@ program MCQC
 ! Hello World
 
   write(*,*)
-  write(*,*) '********************************'
-  write(*,*) '*            Quack             *'
-  write(*,*) '*   __        __        __     *'
-  write(*,*) '* <(o )___  <(o )___  <(o )___ *'
-  write(*,*) '* ( ._> /   ( ._> /   ( ._> /  *'
-  write(*,*) '*|----------------------------|*'
-  write(*,*) '********************************'
+  write(*,*) '******************************************************************************************'
+  write(*,*) '*            Quack                       Quack                         Quack             *'
+  write(*,*) '*   __        __        __       __        __        __       __        __        __     *'
+  write(*,*) '* <(o )___  <(o )___  <(o )___ <(o )___  <(o )___  <(o )___ <(o )___  <(o )___  <(o )___ *'
+  write(*,*) '* ( ._> /   ( ._> /   ( ._> /  ( ._> /   ( ._> /   ( ._> /  ( ._> /   ( ._> /   ( ._> /  *'
+  write(*,*) '*|--------------------------------------------------------------------------------------|*'
+  write(*,*) '******************************************************************************************'
   write(*,*)
 
 ! Which calculations do you want to do?
diff --git a/src/MCQC/NormCoeff.f90 b/src/MCQC/NormCoeff.f90
deleted file mode 100644
index 9e6cabf..0000000
--- a/src/MCQC/NormCoeff.f90
+++ /dev/null
@@ -1,29 +0,0 @@
-function NormCoeff(alpha,a)
-
-  implicit none
-
-! Input variables
-
-  double precision,intent(in)   :: alpha
-  integer,intent(in)            :: a(3)
-
-! local variable
-  double precision              :: pi,dfa(3),dfac
-  integer                       :: atot
-
-! Output variable
-  double precision NormCoeff
-
-  pi = 4d0*atan(1d0)
-  atot = a(1) + a(2) + a(3)
-
-  dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1)))
-  dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2)))
-  dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3)))
-
-
-  NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot
-  NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3))
-  NormCoeff = sqrt(NormCoeff)
-
-end function NormCoeff
diff --git a/src/MCQC/qsGW_PT.f90 b/src/MCQC/qsGW_PT.f90
index 6b23129..adb2239 100644
--- a/src/MCQC/qsGW_PT.f90
+++ b/src/MCQC/qsGW_PT.f90
@@ -16,7 +16,6 @@ subroutine qsGW_PT(nBas,nC,nO,nV,nR,nS,e0,SigCm)
   integer                       :: x,y,z,t
   double precision              :: eps
   double precision,allocatable  :: e1(:),e2(:),e3(:),e4(:)
-  double precision,parameter    :: threshold = 1d-15
 
 ! Allocation
 
diff --git a/src/MCQC/read_auxiliary_basis.f90 b/src/MCQC/read_auxiliary_basis.f90
deleted file mode 100644
index cf4f6ed..0000000
--- a/src/MCQC/read_auxiliary_basis.f90
+++ /dev/null
@@ -1,176 +0,0 @@
-subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, &
-                        TotAngMomShell,KShell,DShell,ExpShell)
-
-! Read auxiliary basis set information
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: NAtoms
-  double precision,intent(in)   :: XYZAtoms(NAtoms,3)
-
-! Local variables
-
-  integer                       :: nShAt,iAt
-  integer                       :: i,j,k
-  character                     :: shelltype
-
-! Output variables
-
-  integer,intent(out)           :: nShell
-  double precision,intent(out)  :: CenterShell(maxShell,3)
-  integer,intent(out)           :: TotAngMomShell(maxShell),KShell(maxShell)
-  double precision,intent(out)  :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
-
-!------------------------------------------------------------------------
-! Primary basis set information
-!------------------------------------------------------------------------
-
-! Open file with basis set specification
-
-  open(unit=2,file='input/basis')
-
-! Read basis information
-
-  write(*,'(A28)') 'Gaussian basis set'
-  write(*,'(A28)') '------------------'
-
-  nShell = 0
-  do i=1,NAtoms
-    read(2,*) iAt,nShAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
-    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
-    write(*,'(A28)') '------------------'
-
-!   Basis function centers
-
-    do j=1,nShAt
-      nShell = nShell + 1
-      do k=1,3
-        CenterShell(nShell,k) = XYZAtoms(iAt,k)
-      enddo
-
-! Shell type and contraction degree
-
-      read(2,*) shelltype,KShell(nShell)
-        if(shelltype == "S") then
-          TotAngMomShell(nShell) = 0
-          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "P") then
-          TotAngMomShell(nShell) = 1
-          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "D") then
-          TotAngMomShell(nShell) = 2
-          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "F") then
-          TotAngMomShell(nShell) = 3
-          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "G") then
-          TotAngMomShell(nShell) = 4
-          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "H") then
-          TotAngMomShell(nShell) = 5
-          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "I") then
-          TotAngMomShell(nShell) = 6
-          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
-        endif
-
-! Read exponents and contraction coefficients
-
-        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
-        do k=1,Kshell(nShell)
-          read(2,*) ExpShell(nShell,k),DShell(nShell,k)
-          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
-        enddo
-    enddo
-    write(*,'(A28)') '------------------'
-  enddo
-
-! Total number of shells
-
-  write(*,'(A28,1X,I16)') 'Number of shells in OBS',nShell
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with basis set specification
-
-  close(unit=2)
-
-!------------------------------------------------------------------------
-! Auxiliary basis set information
-!------------------------------------------------------------------------
-
-! Open file with auxilairy basis specification
-
-  open(unit=3,file='input/auxbasis')
-
-! Read basis information
-
-  write(*,'(A28)') 'Auxiliary basis set'
-  write(*,'(A28)') '-------------------'
-
-  do i=1,NAtoms
-    read(3,*) iAt,nShAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
-    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
-    write(*,'(A28)') '------------------'
-
-!   Basis function centers
-
-    do j=1,nShAt
-      nShell = nShell + 1
-      do k=1,3
-        CenterShell(nShell,k) = XYZAtoms(iAt,k)
-      enddo
-
-! Shell type and contraction degree
-
-      read(3,*) shelltype,KShell(nShell)
-        if(shelltype == "S") then
-          TotAngMomShell(nShell) = 0
-          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "P") then
-          TotAngMomShell(nShell) = 1
-          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "D") then
-          TotAngMomShell(nShell) = 2
-          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "F") then
-          TotAngMomShell(nShell) = 3
-          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "G") then
-          TotAngMomShell(nShell) = 4
-          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "H") then
-          TotAngMomShell(nShell) = 5
-          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "I") then
-          TotAngMomShell(nShell) = 6
-          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
-        endif
-
-! Read exponents and contraction coefficients
-
-        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
-        do k=1,Kshell(nShell)
-          read(3,*) ExpShell(nShell,k),DShell(nShell,k)
-          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
-        enddo
-    enddo
-    write(*,'(A28)') '------------------'
-  enddo
-
-! Total number of shells
-
-  write(*,'(A28,1X,I16)') 'Number of shells in ABS',nShell
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with basis set specification
-
-  close(unit=3)
-
-end subroutine read_auxiliary_basis
diff --git a/src/MCQC/read_basis.f90 b/src/MCQC/read_basis.f90
deleted file mode 100644
index cc82700..0000000
--- a/src/MCQC/read_basis.f90
+++ /dev/null
@@ -1,128 +0,0 @@
-subroutine read_basis(nAt,rAt,nBas,nC,nO,nV,nR,nS, &
-                      nShell,atot,CenterShell,TotAngMomShell,KShell,DShell,ExpShell)
-
-! Read basis set information
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: nAt,nC,nO,nR,atot(maxShell)
-  double precision,intent(in)   :: rAt(nAt,3)
-
-! Local variables
-
-  integer                       :: nShAt,iAt,iShell
-  integer                       :: i,j,k
-  character                     :: shelltype
-
-! Output variables
-
-  integer,intent(out)           :: nShell,nBas,nV,nS
-  double precision,intent(out)  :: CenterShell(maxShell,3)
-  integer,intent(out)           :: TotAngMomShell(maxShell),KShell(maxShell)
-  double precision,intent(out)  :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
-
-!------------------------------------------------------------------------
-! Primary basis set information
-!------------------------------------------------------------------------
-
-! Open file with basis set specification
-
-  open(unit=2,file='input/basis')
-
-! Read basis information
-
-  write(*,'(A28)') 'Gaussian basis set'
-  write(*,'(A28)') '------------------'
-
-  nShell = 0
-  do i=1,nAt
-    read(2,*) iAt,nShAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
-    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
-    write(*,'(A28)') '------------------'
-
-!   Basis function centers
-
-    do j=1,nShAt
-      nShell = nShell + 1
-      do k=1,3
-        CenterShell(nShell,k) = rAt(iAt,k)
-      enddo
-
-! Shell type and contraction degree
-
-      read(2,*) shelltype,KShell(nShell)
-        if(shelltype == "S") then
-          TotAngMomShell(nShell) = 0
-          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "P") then
-          TotAngMomShell(nShell) = 1
-          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "D") then
-          TotAngMomShell(nShell) = 2
-          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "F") then
-          TotAngMomShell(nShell) = 3
-          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "G") then
-          TotAngMomShell(nShell) = 4
-          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "H") then
-          TotAngMomShell(nShell) = 5
-          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
-        elseif(shelltype == "I") then
-          TotAngMomShell(nShell) = 6
-          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
-        endif
-
-! Read exponents and contraction coefficients
-
-        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
-        do k=1,Kshell(nShell)
-          read(2,*) ExpShell(nShell,k),DShell(nShell,k)
-          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
-        enddo
-    enddo
-    write(*,'(A28)') '------------------'
-  enddo
-
-! Total number of shells
-
-  write(*,'(A28,1X,I16)') 'Number of shells in OBS',nShell
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with basis set specification
-
-  close(unit=2)
-
-! Calculate number of basis functions
-
-  nBas = 0
-  do iShell=1,nShell
-    nBas = nBas + (atot(iShell)*atot(iShell) + 3*atot(iShell) + 2)/2
-  enddo
-
-  write(*,'(A28)') '------------------'
-  write(*,'(A28,1X,I16)') 'Number of basis functions',NBas
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Number of virtual orbitals
-
-  nV = nBas - nO 
-
-  if(nR > nV) then
-    write(*,*) 'Number of Rydberg orbitals greater than number of virtual orbitals!'
-    stop
-  endif   
-
-! Number of single excitation
-
-  nS = (nO - nC)*(nV - nR)
-
-
-end subroutine read_basis
diff --git a/src/MCQC/read_geometry.f90 b/src/MCQC/read_geometry.f90
deleted file mode 100644
index 8f0fc56..0000000
--- a/src/MCQC/read_geometry.f90
+++ /dev/null
@@ -1,58 +0,0 @@
-subroutine read_geometry(nAt,ZNuc,rA,ENuc)
-
-! Read molecular geometry
-
-  implicit none
-
-! Ouput variables
-  integer,intent(in)            :: nAt
-
-! Local variables
-  integer                       :: i,j
-  double precision              :: RAB
-
-! Ouput variables
-  double precision,intent(out)  :: ZNuc(NAt),rA(nAt,3),ENuc
-  
-
-! Open file with geometry specification
-  open(unit=1,file='input/molecule')
-
-! Read number of atoms
-  read(1,*) 
-  read(1,*) 
-  read(1,*) 
-
-  do i=1,nAt
-    read(1,*) ZNuc(i),rA(i,1),rA(i,2),rA(i,3)
-  enddo
-
-! Compute nuclear repulsion energy
-  ENuc = 0
-
-  do i=1,nAt-1
-    do j=i+1,nAt
-      RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2
-      ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
-    enddo
-  enddo
-
-! Close file with geometry specification
-  close(unit=1)
-
-! Print geometry
-  write(*,'(A28)') '------------------'
-  write(*,'(A28)') 'Molecular geometry'
-  write(*,'(A28)') '------------------'
-  do i=1,NAt
-    write(*,'(A28,1X,I16)') 'Atom n. ',i
-    write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
-    write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,3)
-  enddo
-  write(*,*)
-  write(*,'(A28)') '------------------'
-  write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-end subroutine read_geometry
diff --git a/src/MCQC/read_integrals.f90 b/src/MCQC/read_integrals.f90
index 0226865..fa435d1 100644
--- a/src/MCQC/read_integrals.f90
+++ b/src/MCQC/read_integrals.f90
@@ -13,7 +13,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
   logical                       :: debug
   integer                       :: mu,nu,la,si
   double precision              :: Ov,Kin,Nuc,ERI
-  double precision              :: scale
 
 ! Output variables
 
@@ -23,8 +22,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
 
   debug = .false.
 
-  scale = 1d0
-
   open(unit=8 ,file='int/Ov.dat')
   open(unit=9 ,file='int/Kin.dat')
   open(unit=10,file='int/Nuc.dat')
@@ -44,7 +41,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
   T = 0d0
   do 
     read(9,*,end=9) mu,nu,Kin
-    T(mu,nu) = Kin/scale**2
+    T(mu,nu) = Kin
   enddo
   9 close(unit=9)
 
@@ -65,25 +62,22 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
 
   G = 0d0
   do 
-    read(11,*,end=11) mu,la,nu,si,ERI
-!   read(11,*,end=11) ERI,mu,nu,la,si
-
-    ERI = ERI/scale
-!   <12|34>
+    read(11,*,end=11) mu,nu,la,si,ERI
+!   (12|34)
     G(mu,nu,la,si) = ERI
-!   <32|14>
-    G(la,nu,mu,si) = ERI
-!   <14|32>
-    G(mu,si,la,nu) = ERI
-!   <34|12>
-    G(la,si,mu,nu) = ERI
-!   <41|23>
-    G(si,mu,nu,la) = ERI
-!   <23|41>
-    G(nu,la,si,mu) = ERI
-!   <21|43>
+!   (21|34)
+    G(nu,mu,la,si) = ERI
+!   (12|43)
+    G(mu,nu,si,la) = ERI
+!   (21|43)
     G(nu,mu,si,la) = ERI
-!   <43|21>
+!   (34|12)
+    G(la,si,mu,nu) = ERI
+!   (43|12)
+    G(si,la,mu,nu) = ERI
+!   (34|21)
+    G(la,si,nu,mu) = ERI
+!   (43|21)
     G(si,la,nu,mu) = ERI
   enddo
   11 close(unit=11)
diff --git a/src/MCQC/read_molecule.f90 b/src/MCQC/read_molecule.f90
deleted file mode 100644
index b7f14b9..0000000
--- a/src/MCQC/read_molecule.f90
+++ /dev/null
@@ -1,67 +0,0 @@
-subroutine read_molecule(nAt,nEl,nC,nO,nR)
-
-! Read number of atoms nAt, 
-!      number of electrons nEl, 
-!      number of core electrons nC,
-!      number of Rydberg orbitals nR
-
-  implicit none
-
-! Input variables
-  integer,intent(out)           :: nAt,nEl,nC,nO,nR
-
-! Open file with geometry specification
-
-  open(unit=1,file='input/molecule')
-
-! Read number of atoms and number of electrons
-
-  read(1,*) 
-  read(1,*) nAt,nEl,nC,nR
-
-! Number of occupied orbitals
-
-  if(mod(nEl,2) /= 0) then
-    write(*,*) 'closed-shell system required!'
-!   stop
-  endif
-  nO = nEl/2
-
-! Number of core orbitals
-
-  if(mod(nC,2) /= 0) then
-    write(*,*) 'Number of core electrons not even!'
-    stop
-  endif
-  nC = nC/2
-
-  if(nC > nO) then
-    write(*,*) 'Number of core electrons greater than number of electrons!'
-    stop
-  endif
-
-! Print results
-
-  write(*,'(A28)') '----------------------'
-  write(*,'(A28,1X,I16)') 'Number of atoms',nAt
-  write(*,'(A28)') '----------------------'
-  write(*,*)
-  write(*,'(A28)') '----------------------'
-  write(*,'(A28,1X,I16)') 'Number of electrons',nEl
-  write(*,'(A28)') '----------------------'
-  write(*,*)
-  write(*,'(A28)') '----------------------'
-  write(*,'(A28,1X,I16)') 'Number of core electrons',2*nC
-  write(*,'(A28)') '----------------------'
-  write(*,*)
-  write(*,*)
-  write(*,'(A28)') '----------------------'
-  write(*,'(A28,1X,I16)') 'Number of Rydberg orbitals',nR
-  write(*,'(A28)') '----------------------'
-  write(*,*)
-
-! Close file with geometry specification
-
-  close(unit=1)
-
-end subroutine read_molecule
diff --git a/src/xcDFT/eDFT.f90 b/src/xcDFT/eDFT.f90
index 6ae8a60..9db226b 100644
--- a/src/xcDFT/eDFT.f90
+++ b/src/xcDFT/eDFT.f90
@@ -4,10 +4,10 @@ program eDFT
 
   include 'parameters.h'
 
-  integer                       :: nAt,nBas,nEl(nspin),nO(nspin),nV(nspin)
+  integer                       :: nNuc,nBas,nEl(nspin),nO(nspin),nV(nspin)
   double precision              :: ENuc,EKS
 
-  double precision,allocatable  :: ZNuc(:),rAt(:,:)
+  double precision,allocatable  :: ZNuc(:),rNuc(:,:)
 
   integer                       :: nShell
   integer,allocatable           :: TotAngMomShell(:)
@@ -55,12 +55,12 @@ program eDFT
 ! nBas = number of basis functions (see below)
 !      = nO + nV
 
-  call read_molecule(nAt,nEl,nO)
-  allocate(ZNuc(nAt),rAt(nAt,ncart))
+  call read_molecule(nNuc,nEl,nO)
+  allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
 
 ! Read geometry
 
-  call read_geometry(nAt,ZNuc,rAt,ENuc)
+  call read_geometry(nNuc,ZNuc,rNuc,ENuc)
 
   allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
            DShell(maxShell,maxK),ExpShell(maxShell,maxK))
@@ -69,7 +69,7 @@ program eDFT
 ! Read basis set information
 !------------------------------------------------------------------------
 
-  call read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
+  call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
 
 !------------------------------------------------------------------------
 ! Read one- and two-electron integrals
diff --git a/src/xcDFT/read_integrals.f90 b/src/xcDFT/read_integrals.f90
deleted file mode 100644
index fa435d1..0000000
--- a/src/xcDFT/read_integrals.f90
+++ /dev/null
@@ -1,114 +0,0 @@
-subroutine read_integrals(nBas,S,T,V,Hc,G)
-
-! Read one- and two-electron integrals from files
-
-  implicit none
-
-! Input variables
-
-  integer,intent(in)            :: nBas
-
-! Local variables
-
-  logical                       :: debug
-  integer                       :: mu,nu,la,si
-  double precision              :: Ov,Kin,Nuc,ERI
-
-! Output variables
-
-  double precision,intent(out)  :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
-
-! Open file with integrals
-
-  debug = .false.
-
-  open(unit=8 ,file='int/Ov.dat')
-  open(unit=9 ,file='int/Kin.dat')
-  open(unit=10,file='int/Nuc.dat')
-  open(unit=11,file='int/ERI.dat')
-
-! Read overlap integrals
-
-  S = 0d0
-  do 
-    read(8,*,end=8) mu,nu,Ov
-    S(mu,nu) = Ov
-  enddo
-  8 close(unit=8)
-
-! Read kinetic integrals
-
-  T = 0d0
-  do 
-    read(9,*,end=9) mu,nu,Kin
-    T(mu,nu) = Kin
-  enddo
-  9 close(unit=9)
-
-! Read nuclear integrals
-
-  V = 0d0
-  do 
-    read(10,*,end=10) mu,nu,Nuc
-    V(mu,nu) = Nuc
-  enddo
-  10 close(unit=10)
-
-! Define core Hamiltonian
-
-  Hc = T + V
-
-! Read nuclear integrals
-
-  G = 0d0
-  do 
-    read(11,*,end=11) mu,nu,la,si,ERI
-!   (12|34)
-    G(mu,nu,la,si) = ERI
-!   (21|34)
-    G(nu,mu,la,si) = ERI
-!   (12|43)
-    G(mu,nu,si,la) = ERI
-!   (21|43)
-    G(nu,mu,si,la) = ERI
-!   (34|12)
-    G(la,si,mu,nu) = ERI
-!   (43|12)
-    G(si,la,mu,nu) = ERI
-!   (34|21)
-    G(la,si,nu,mu) = ERI
-!   (43|21)
-    G(si,la,nu,mu) = ERI
-  enddo
-  11 close(unit=11)
-
-
-! Print results
-  if(debug) then
-    write(*,'(A28)') '----------------------'
-    write(*,'(A28)') 'Overlap integrals'
-    write(*,'(A28)') '----------------------'
-    call matout(nBas,nBas,S)
-    write(*,*)
-    write(*,'(A28)') '----------------------'
-    write(*,'(A28)') 'Kinetic integrals'
-    write(*,'(A28)') '----------------------'
-    call matout(nBas,nBas,T)
-    write(*,*)
-    write(*,'(A28)') '----------------------'
-    write(*,'(A28)') 'Nuclear integrals'
-    write(*,'(A28)') '----------------------'
-    call matout(nBas,nBas,V)
-    write(*,*)
-    write(*,'(A28)') '----------------------'
-    write(*,'(A28)') 'Electron repulsion integrals'
-    write(*,'(A28)') '----------------------'
-    do la=1,nBas
-      do si=1,nBas
-        call matout(nBas,nBas,G(1,1,la,si))
-      enddo
-    enddo
-    write(*,*)
-  endif
-
-end subroutine read_integrals
diff --git a/src/xcDFT/read_options.f90 b/src/xcDFT/read_options.f90
index ff96627..38574f8 100644
--- a/src/xcDFT/read_options.f90
+++ b/src/xcDFT/read_options.f90
@@ -31,7 +31,7 @@ subroutine read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DI
 
 ! Open file with method specification
 
-  open(unit=1,file='input/options')
+  open(unit=1,file='input/dft')
 
 ! Default values