diff --git a/README b/README
new file mode 100644
index 0000000000000000000000000000000000000000..39b44e59ea30a7693e97501e38ff90347c0b68c5
--- /dev/null
+++ b/README
@@ -0,0 +1,3 @@
+This is a project based on DFTATOM (Certik 2013) and edited by Levämäki
+(2015) for orbital-free calculations. The idea of the project here is to
+add the necessary gamma-dependent terms.
\ No newline at end of file
diff --git a/main_Be b/main_Be
index d4c3dbe03c863711423bcdfd73edf875ffb0a950..a4a298f7e03a9f7289e349a1fdbd1ebbb4ee4db5 100755
Binary files a/main_Be and b/main_Be differ
diff --git a/solverhovp.f90 b/solverhovp.f90
index e6eeb204eba3e50005501f627530b181ccb4df34..cfe7b383d21e693ba8cddb965420f7197938dc74 100644
--- a/solverhovp.f90
+++ b/solverhovp.f90
@@ -336,6 +336,7 @@ contains
     end do
     do i=1,AdamsOutwardOrder
        call DiffBeta(Pdot(i,:),Qdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
+       call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
        !print *,Qdot(:,i)
     end do
     endPoint=ctp-1
@@ -369,6 +370,7 @@ contains
           Q(i+1,j)=M21*Ptmp+M22*Qtmp
        end do
        call DiffBeta(Pdot(i+1,:),Qdot(i+1,:),orbs(:,i+1),P(i+1,:),Q(i+1,:),e1s,mu)
+       call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
        do PECEiter=1,PECEiterMax
           do j=1,NbgGrid
              fvPOutward(:,j)=Dr(i-AdamsOutwardOrder+2:i+1)*4*((mu-vP(i-AdamsOutwardOrder+2:i+1,j))*&
@@ -398,7 +400,7 @@ contains
              Q(i+1,j)=M21*Ptmp+M22*Qtmp
           end do
           call DiffBeta(Pdot(i+1,:),Qdot(i+1,:),orbs(:,i+1),P(i+1,:),Q(i+1,:),e1s,mu)
-
+		  call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
           !print *,'PECEiter, i+1, P(i+1,1) = ',PECEiter,i,P(i+1,1)
 
        end do
@@ -462,6 +464,7 @@ contains
     end do
     do i=1,AdamsOutwardOrder
        call DiffBeta(Pdot(i,:),Qdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
+	   call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
        !print *,Qdot(:,i)
     end do
     endPoint=ctp-1
@@ -494,6 +497,7 @@ contains
           Q(i+1,j)=M21*Ptmp+M22*Qtmp
        end do
        call DiffBeta(Pdot(i+1,:),Qdot(i+1,:),orbs(:,i+1),P(i+1,:),Q(i+1,:),e1s,mu)
+       call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
        do PECEiter=1,PECEiterMax
           do j=1,NbgGrid
              fvPOutward(:,j)=Dr(i-AdamsOutwardOrder+2:i+1)*4*((mu-vP(i-AdamsOutwardOrder+2:i+1,j))*&
@@ -523,6 +527,7 @@ contains
              Q(i+1,j)=M21*Ptmp+M22*Qtmp
           end do
           call DiffBeta(Pdot(i+1,:),Qdot(i+1,:),orbs(:,i+1),P(i+1,:),Q(i+1,:),e1s,mu)
+          call DiffGamma(Pdot(i,:),Pgdot(i,:),Qdot(i,:),Qgdot(i,:),orbs(:,i),P(i,:),Q(i,:),e1s,mu)
        end do
        !print *,orbs(:,i+1)
     end do
@@ -720,5 +725,18 @@ contains
     !
     !print *,Qdot
   end subroutine DiffBeta
+  subroutine DiffGamma(Pdot, Pgdot, Qdot, Qgdot,orbs,P,Q, e1s,mu)
+  	use data, only: betaGrid, gammaGrid,NbgGrid, Norbs,leigens,lambdas
+  	real(dp), intent(in) :: P(:),Q(:), Pdot(:), Qdot(:), e1s, mu
+    real(dp), intent(inout) :: Pgdot(:),Qgdot(:),orbs(:)
+  	real(dp) :: orbsPrime(Norbs)
+  	integer :: N,NRHS,LDA,LDB,INFO
+  	real(dp) :: A(Norbs,Norbs),Awork(Norbs,Norbs),B(Norbs,2) 
+  	integer :: IPIV(Norbs)
+  	real(dp) :: rho(Norbs),rhoPrime(Norbs),energies(Norbs)
+  	integer :: i
 
+	Pgdot = 0.0_dp
+	Qgdot = 0.0_dp
+  end subroutine DiffGamma
 end module solverhovp