diff --git a/GoDuck b/GoDuck
index dcf207b..97684a7 100755
--- a/GoDuck
+++ b/GoDuck
@@ -11,8 +11,9 @@ 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/IntPak
   ./bin/QuAcK
 fi
 
diff --git a/GoQCaml b/GoQCaml
new file mode 100755
index 0000000..881460c
--- /dev/null
+++ b/GoQCaml
@@ -0,0 +1,4 @@
+#! /bin/bash
+
+cd int
+../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz
diff --git a/examples/molecule.N2 b/examples/molecule.N2
index 93b448f..35e31a1 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.0
diff --git a/input/molecule b/input/molecule
index 93b448f..35e31a1 100644
--- a/input/molecule
+++ b/input/molecule
@@ -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.0
diff --git a/src/QuAcK/ACFDT.f90 b/src/QuAcK/ACFDT.f90
new file mode 100644
index 0000000..ff31901
--- /dev/null
+++ b/src/QuAcK/ACFDT.f90
@@ -0,0 +1,148 @@
+subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold, & 
+                 nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
+
+! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
+
+  implicit none
+  include 'parameters.h'
+  include 'quadrature.h'
+
+! Input variables
+
+  logical,intent(in)            :: doXBS
+  logical,intent(in)            :: exchange_kernel
+  logical,intent(in)            :: dRPA
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: BSE
+  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: triplet_manifold
+
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: e(nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+  double precision              :: Omega(nS,nspin)
+  double precision              :: XpY(nS,nS,nspin)
+  double precision              :: XmY(nS,nS,nspin)
+  double precision              :: rho(nBas,nBas,nS,nspin)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: iAC
+  double precision              :: lambda
+  double precision,allocatable  :: Ec(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcAC(nspin)
+
+! Memory allocation
+
+  allocate(Ec(nAC,nspin))
+
+! Antisymmetrized kernel version
+
+  if(exchange_kernel) then
+
+    write(*,*)
+    write(*,*) '*** Exchange kernel version ***'
+    write(*,*)
+
+  end if
+
+! Singlet manifold
+
+  if(singlet_manifold) then
+
+    ispin = 1
+    EcAC(ispin) = 0d0
+    Ec(:,ispin) = 0d0
+
+    write(*,*) '--------------'
+    write(*,*) 'Singlet states'
+    write(*,*) '--------------'
+    write(*,*) 
+
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
+    write(*,*) '-----------------------------------------------------------------------------------'
+
+    do iAC=1,nAC
+ 
+      lambda = rAC(iAC)
+
+      if(doXBS) then
+
+        call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
+                             rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
+
+      end if
+
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
+                       rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+
+      call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
+
+      write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
+
+    end do
+
+    EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
+
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,*)
+
+  end if
+ 
+! Triplet manifold
+
+  if(triplet_manifold) then
+
+    ispin = 2
+    EcAC(ispin) = 0d0
+    Ec(:,ispin) = 0d0
+
+    write(*,*) '--------------'
+    write(*,*) 'Triplet states'
+    write(*,*) '--------------'
+    write(*,*) 
+
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
+    write(*,*) '-----------------------------------------------------------------------------------'
+
+    do iAC=1,nAC
+ 
+      lambda = rAC(iAC)
+
+      if(doXBS) then
+
+        call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
+                             rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
+
+      end if  
+
+      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
+                           rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+
+      call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
+
+      write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
+
+    end do
+
+    EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
+
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
+    write(*,*) '-----------------------------------------------------------------------------------'
+    write(*,*)
+
+  end if
+
+end subroutine ACFDT
diff --git a/src/QuAcK/ACFDT_correlation_energy.f90 b/src/QuAcK/ACFDT_correlation_energy.f90
new file mode 100644
index 0000000..f43b87d
--- /dev/null
+++ b/src/QuAcK/ACFDT_correlation_energy.f90
@@ -0,0 +1,93 @@
+subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC)
+
+! Compute the correlation energy via the adiabatic connection formula
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: ispin
+  logical,intent(in)            :: exchange_kernel
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: XpY(nS,nS)
+  double precision,intent(in)   :: XmY(nS,nS)
+
+! Local variables
+
+  integer                       :: i,j,a,b
+  integer                       :: ia,jb,kc
+  double precision              :: delta_spin
+  double precision              :: delta_Kx
+  double precision,allocatable  :: P(:,:)
+  double precision,allocatable  :: Ap(:,:)
+  double precision,allocatable  :: Bp(:,:)
+  double precision,allocatable  :: X(:,:)
+  double precision,allocatable  :: Y(:,:)
+  double precision,external     :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)  :: EcAC
+
+! Singlet or triplet manifold?
+
+  delta_spin = 0d0
+  if(ispin == 1) delta_spin = +1d0
+  if(ispin == 2) delta_spin = -1d0
+
+! Exchange kernel
+
+  delta_Kx = 0d0
+  if(exchange_kernel) delta_Kx = 1d0
+  
+! Memory allocation
+
+  allocate(P(nS,nS),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS))
+
+! Compute P = (X+Y)^T.(X+Y) - 1
+
+  P(:,:) = matmul(transpose(XpY),XpY)
+
+  do ia=1,nS
+    P(ia,ia) = P(ia,ia) - 1d0
+  enddo
+
+! Compute Aiajb = (ia|bj) and Biajb = (ia|jb)
+
+  ia = 0
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      ia = ia + 1
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nBas-nR
+          jb = jb + 1
+
+            Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) & 
+                      - delta_Kx*ERI(i,b,j,a)
+
+            Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) &
+                      - delta_Kx*ERI(i,j,b,a)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+! Compute Tr(K x P_lambda)
+
+! EcAC = trace_matrix(nS,matmul(Ap,P))
+
+  X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
+  Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
+
+  EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
+       + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
+       - trace_matrix(nS,Ap)
+
+! call matout(nS,nS,matmul(transpose(X),X) - matmul(transpose(Y),Y))
+
+end subroutine ACFDT_correlation_energy
+
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index dc61dec..7703f34 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -205,9 +205,8 @@ program QuAcK
 
   else
 
-!   call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
-    call 
-!   call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
+    call system('./GoQCaml')
+    call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
 
   end if