Index: src/minv.f90
===================================================================
--- src/minv.f90	(revision 2)
+++ src/minv.f90	(revision 3)
@@ -4,8 +4,8 @@
       integer, parameter :: KINT = kind(1)
       integer, parameter :: KREAL = kind(1.0d0)
 
+      INTEGER(KINT), INTENT(IN) :: N
       REAL(KREAL), INTENT(INOUT) :: A(N*N)
-      INTEGER(KINT), INTENT(IN) :: N
       REAL(KREAL), INTENT(OUT) :: D
 
  
Index: src/refsor.f90
===================================================================
--- src/refsor.f90	(revision 2)
+++ src/refsor.f90	(revision 3)
@@ -39,6 +39,36 @@
 
    use VarTypes
 
+  interface
+
+     Recursive Subroutine D_subsor (XDONT, IDEB1, IFIN1)
+
+       use VarTypes 
+
+       IMPLICIT NONE
+
+       Real(KREAL), dimension (:), Intent (InOut) :: XDONT
+       Integer(KINT), Intent (In) :: IDEB1, IFIN1
+
+       Integer(KINT), Parameter :: NINS = 16 ! Max for insertion sort
+       Integer(KINT) :: ICRS, IDEB, IDCR, IFIN, IMIL
+       Real(KREAL) :: XPIV, XWRK
+     end Subroutine D_subsor
+
+     Subroutine D_inssor (XDONT)
+
+       use VarTypes 
+  
+       IMPLICIT NONE
+       Real(KREAL), dimension (:), Intent (InOut) :: XDONT
+       
+       Integer(KINT) :: ICRS, IDCR
+       Real(KREAL) :: XWRK
+
+     End Subroutine D_inssor
+
+  end interface
+
       Real (KREAL), Dimension (:), Intent (InOut) :: XDONT
 ! __________________________________________________________
 !
Index: src/Makefile
===================================================================
--- src/Makefile	(revision 2)
+++ src/Makefile	(revision 3)
@@ -5,7 +5,7 @@
 # You might also have to edit the location of
 # some libraries (like mkl for ifort)
 
-Machine=arqg
+Machine=gfortran
 
 ###########################################
 #                                         #
@@ -13,30 +13,39 @@
 #                                         #
 ###########################################
 ifeq ($(Machine),ifort)
-# Flags for Ifort
-COMP=/opt/intel/fc/9.1.036/bin/ifort
-F90=${COMP}  
-LINK=${COMP} -Vaxlib   -L/usr/lib/  -L/opt/intel/mkl/8.1/lib/32/ -lguide -lpthread -lmkl_lapack -lmkl_ia32  -lmkl
+# Flags for arq Ifort
+COMP=ifort
+#F90=${COMP}   -g -check bounds  -check format -check uninit  -traceback
+F90=${COMP}   -g  -traceback
+F77=${F90}
+#F90=${COMP}  
+LIBMATH=  -L/opt/intel/Compiler/11.1/064/mkl/lib/64  \
+    -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential \
+    -lmkl_core -Wl,--end-group 
+LINK=${COMP} -lguide -lpthread   -L/usr/lib/ 
 FreeF=-FR
 endif
 
+
 ifeq ($(Machine),g95)
 # Flags for g95 
-COMP=g95 -g  -save-temps -Wall -fbounds-check
+COMP=g95 -g   -Wall -fbounds-check
 #COMP=g95  -fbounds-check -g -ftrace=full -ftrace=frame  -save-temps
 F90=${COMP}
 F77=${COMP}
-LINK=${COMP}   -L/usr/lib/ -llapack -lblas
+LIBMATH= -llapack -lblas
+LINK=${COMP}   -L/usr/lib/
 FreeF=-ffree-form
 endif
 
 ifeq ($(Machine),gfortran)
 # Flags for gfortran
-COMP=gfortran  -fbounds-check -pedantic -std=gnu 
+COMP=gfortran  -pedantic -std=gnu 
 #COMP=g95  -fbounds-check -g -ftrace=full -ftrace=frame  -save-temps
 F90=${COMP}
 F77=${COMP}
-LINK=${COMP} -L/usr/lib/  -llapack -lblas
+LIBMATH= -llapack -lblas
+LINK=${COMP}   -L/usr/lib/
 FreeF=-ffree-form
 endif
 
@@ -46,7 +55,8 @@
 COMP=pgf90
 FLAGS= -mp  -fast -Kieee -fastsse  -tpp7 -Mipa=fast
 F90=$(COMP)  -c  $(FLAGS)
-LINK=pgf90 $(FLAGS)  -L/usr/local/pgi/linux86/6.2/lib -llapack  -lblas -lpthread 
+LIBMATH= -llapack -lblas
+LINK=${COMP} $(FLAGS)  -L/usr/local/pgi/linux86/6.2/lib  -lpthread 
 FreeF=-Mfree
 endif
 
@@ -54,7 +64,8 @@
 # Flags for PathScale
 COMP=/softs/pathscale/bin/pathf90
 F90=${COMP}
-LINK=${COMP}   -L/usr/lib/ -llapack -lblas
+LIBMATH= -llapack -lblas
+LINK=${COMP}   -L/usr/lib/ 
 FreeF=-FR
 endif
 
@@ -62,7 +73,8 @@
 # Flags for Xlf
 COMP=xlf
 F90=${COMP}
-LINK=${COMP}   -L/usr/lib/ -llapack -lblas
+LIBMATH= -llapack -lblas
+LINK=${COMP}   -L/usr/lib/ 
 FreeF=-qfree
 endif
 
@@ -100,7 +112,9 @@
 COMP=ifort
 F90=${COMP}
 F77=${COMP}
-LINK=${COMP}  -lguide -lpthread
+LIBMATH=  -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential \
+           -lmkl_core -Wl,--end-group
+LINK=${COMP}  -lguide -lpthread 
 FreeF=-FR
 endif
 
@@ -139,10 +153,10 @@
 F90=${COMP}   -g -check all  -traceback
 F77=${F90}
 #F90=${COMP}  
-LINK=${COMP} -lguide -lpthread   -L/usr/lib/ \
-  -L/opt/intel/Compiler/11.1/064/mkl/lib/64  \
+LIBMATH= -L/opt/intel/Compiler/11.1/064/mkl/lib/64  \
   -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential \
-           -lmkl_core -Wl,--end-group -lpthread 
+           -lmkl_core -Wl,--end-group
+LINK=${COMP} -lguide -lpthread   -L/usr/lib/
 FreeF=-FR
 endif
 
@@ -152,7 +166,9 @@
 F90=${COMP}  
 F77=${F90}
 #F90=${COMP}  
-LINK=${COMP}  -L/usr/lib/   -llapack -lblas
+F90=${COMP}
+LIBMATH= -llapack -lblas
+LINK=${COMP}  -L/usr/lib/  
 FreeF=-ffree-form
 endif
 
@@ -170,10 +186,11 @@
 F90=${COMP}   -g -check all -traceback
 F77=${F90}
 #F90=${COMP}  
-LINK=${COMP}    -L/softs/intel/v12.0.084/mkl/lib/intel64  \
-  -L/softs/intel/v12.0.084/lib/intel64  -L/usr/lib \
+LIBMATH=    -L/softs/intel/v12.0.084/mkl/lib/intel64  \
+ -L/softs/intel/v12.0.084/lib/intel64  -L/usr/lib \
   -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential \
            -lmkl_core -Wl,--end-group
+LINK=${COMP}
 FreeF=-FR
 endif
 
@@ -184,35 +201,13 @@
 F90=${COMP}   -g -check all -traceback
 F77=${F90}
 #F90=${COMP}  
-LINK=${COMP} -lguide -lpthread   -L/usr/lib/ \
-  -L/softs/intel/v11.1.069/mkl/lib/em64t  \
+LIBMATH=   -L/softs/intel/v11.1.069/mkl/lib/em64t  \
   -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential \
-           -lmkl_core -Wl,--end-group -lpthread
+           -lmkl_core -Wl,--end-group 
+LINK=${COMP} -lguide -lpthread   -L/usr/lib/
 FreeF=-FR
 endif
 
-
-ifeq ($(Machine),duo)
-# Flags for DUO Ifort
-COMP=/opt/intel/fc/9.1.037/bin/ifort
-#F90=${COMP}   -g -check all
-F90=${COMP}  
-LINK=${COMP} -Vaxlib   -L/usr/lib/  -L/opt/intel/mkl/8.1/lib/32/ -lguide -lpthread  -lmkl_lapack64 -lmkl
-FreeF=-FR
-endif
-
-ifeq ($(Machine),commander)
-# Flags for Commander Ifort
-COMP=/opt/intel/fc/9.0/bin/ifort
-#F90=${COMP}   -g -check all
-F90=${COMP}
-LINK=${COMP} -Vaxlib   -L/usr/lib/  -L/opt/intel/mkl/8.1/lib/32/ -lguide -lpthread -lmkl_lapack -lmkl_ia32  -lmkl
-FreeF=-FR
-endif
-
-Version=1.0
-
-
 ###############################################################
 #
 # Help Makefile 
@@ -248,6 +243,7 @@
       egrad_vasp.f90 \
       egrad_turbomole.f90 \
       egrad_chamfre.f90 \
+      egrad_LEPS.f90 \
       CalcHess.f90 \
       Hupdate_all.f90 \
       Hupdate_MS.f90 \
@@ -341,8 +337,8 @@
 SRCLAPACKU=$(wildcard ./lapack/util/*.f )
 SRCBLAS=$(wildcard ./blas/*.f )
 
-OBJLAPACK=$(SCRLAPACKD:.f=.o) \
-        $(SCRLAPACKU:.f=.o)
+OBJLAPACK=$(SRCLAPACKD:.f=.o) \
+        $(SRCLAPACKU:.f=.o)
 
 OBJBLAS=$(SRCBLAS:.f=.o)
 
@@ -355,7 +351,7 @@
 	ln -s ${PWD}/../utils/xyz2path ../exe/.
 
 path: ${MODo} ${OBJ} 
-	${LINK} -o Path.exe ${OBJ} ${MODo} 
+	${LINK} -o Path.exe ${OBJ} ${MODo} ${LIBMATH} 
 	@echo ""
 	@echo "Path.exe has been created."
 
@@ -389,7 +385,7 @@
 	ar rcs ./blas/blas.a ${OBJBLAS}
 
 clean:
-	rm -f *~ "#"* *.s ${OBJ} ${MOD}  ${MODo} ${OBJLAPACK} ${OBJBLAS}  
+	rm -f *~ "#"* *.s ${OBJ} *.mod  ${MODo} ${OBJLAPACK} ${OBJBLAS}  
 
 veryclean: clean
 	rm -f Path.exe ./lapack/lapack.a  ./blas/blas.a ../utils/xyz2scan ../utils/xyz2path
Index: src/Path.f90
===================================================================
--- src/Path.f90	(revision 2)
+++ src/Path.f90	(revision 3)
@@ -590,6 +590,7 @@
 
   if (Prog.EQ.'TM') Prog="TURBOMOLE"
   if (Prog.EQ.'TEST2') Prog="CHAMFRE"
+  if (Prog.EQ.'TEST3') Prog="LEPS"
 
 
   Select Case (Prog)
Index: src/egrad_LEPS.f90
===================================================================
--- src/egrad_LEPS.f90	(revision 0)
+++ src/egrad_LEPS.f90	(revision 3)
@@ -0,0 +1,146 @@
+  SUBROUTINE egrad_LEPS(nat,e,geom,grad) 
+! This program computes the ernergy and gradien in cartesian coordinates
+! for the cartesian geometry Geom
+
+   IMPLICIT NONE
+     integer, parameter :: KINT = kind(1)
+     integer, parameter :: KREAL = kind(1.0d0)
+
+! Number of atoms
+       INTEGER(KINT), INTENT(IN) :: Nat
+! Input geometry
+       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
+! output energy and gradient
+       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)
+
+
+!     Bohr --> Angstr
+      real(KREAL), parameter :: BOHR   = 0.52917726D+00
+!
+! Parameters to define the surface
+      INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)
+      INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)
+      REAL(KREAL), PARAMETER :: hh=0.001d0
+      
+! Variables
+       INTEGER(KINT) :: i,j,iat,jat
+       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
+       REAL(KREAL) :: xp,yp
+
+       CHARACTER(132) :: Line
+       LOGICAL :: Debug
+  
+      real(KREAL), parameter :: zero = 0.0_KREAL,   one = 1.0_KREAL
+      REAL(KREAL), external :: ELEPS_xyz
+
+
+  INTERFACE
+     function valid(string) result (isValid)
+       CHARACTER(*), intent(in) :: string
+       logical                  :: isValid
+     END function VALID
+     
+  END INTERFACE
+
+  debug=valid('egrad_leps')
+
+ if (debug) WRITE(*,*) '================ Entering Egrad_leps ==================='
+ if (debug) THEN
+    WRITE(*,*) "Cartesian Geometry"
+    DO I=1,Nat
+       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3)
+    END DO
+ END IF
+      
+      ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3))
+
+      e=ELEPS_xyz(nat,Geom)
+
+! We now calculate the gradients using numerical derivatives
+      Grad=0.d0
+      GradTmp=0.d0
+      do iat=1,3
+         do jat=1,nat
+            xyztmp=geom
+            do i=1,6
+               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh
+               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*ELEPS_xyz(nat,xyztmp)
+            end do
+         end do
+      end do
+      gradTmp=gradTmp/(60.*hh)
+      do iat=1,nat
+       grad(3*iat-2:3*iat)=gradTmp(iat,1:3)
+      end do
+
+ if (debug) THEN
+    WRITE(*,*) "Cartesian gradient "
+    DO I=1,Nat
+       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)
+    END DO
+ END IF
+
+
+      deallocate(xyztmp,gradTmp)
+
+
+ if (debug) WRITE(*,*) '================ Exiting Egrad_leps ==================='
+
+! ======================================================================
+      end 
+
+
+
+      function ELEPS_xyz(natoms,Xyz)
+
+        use Path_module, only : order
+
+   IMPLICIT NONE
+     integer, parameter :: KINT = kind(1)
+     integer, parameter :: KREAL = kind(1.0d0)
+
+
+      INTEGER(KINT) ,INTENT(IN) :: natoms
+      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
+      REAL(KREAL) :: ELEPS_xyz
+
+      INTEGER(KINT) :: i
+      REAL(KREAL) :: rAB, rBC, rAC
+
+      REAL(KREAL), PARAMETER :: a=0.05, b=0.30, c=0.05
+      REAL(KREAL), PARAMETER :: dAB=4.746, dBC=4.747, dAC=3.445
+      REAL(KREAL), PARAMETER :: r0=0.742, alpha=1.942
+
+      REAL(KREAL) :: Qbc,QAc,Qab,Jab,Jbc,Jac
+
+      Real(KREAL), Parameter ::  autoA=0.52917715d0,autoeV=27.21183d0
+
+      rAB=0.
+      rAC=0.
+      rBC=0.
+      do i=1,3
+         rBC=rBC+(xyz(Order(3),i)-xyz(order(2),i))**2
+         rAB=rAB+(xyz(order(1),i)-xyz(order(2),i))**2
+         rAC=rAC+(xyz(order(1),i)-xyz(order(3),i))**2
+      end do
+      rAB=sqrt(rAB)
+      rBC=sqrt(rBC)
+      rAC=sqrt(rAC)
+
+! For consistency with previous articles : rAC is constrained to be 3.742
+      rAC=3.742
+
+      Qab = dab*(1.5*exp(-2*alpha*(rAB-r0))-exp(-alpha*(rAB-r0)))/(2+2*a)
+      Qbc = dbc*(1.5*exp(-2*alpha*(rBC-r0))-exp(-alpha*(rBC-r0)))/(2+2*b)
+      Qac = dac*(1.5*exp (-2*alpha*(rAC-r0))-exp(-alpha*(rAC-r0)))/(2+2*c)
+      Jab = dab*(exp(-2*alpha*(rAB-r0))-6*exp(-alpha*(rAB-r0)))/(4+4*a)
+      Jbc = dbc*(exp(-2*alpha*(rBC-r0))-6*exp(-alpha*(rBC-r0)))/(4+4*b)
+      Jac = dac*(exp(-2*alpha*(rAC-r0))-6*exp(-alpha*(rAC-r0)))/(4+4*c)
+
+! V(x,y) = Qab(x)+Qbc(y)+Qac(rac)-sqrt(Jab(x)**2+Jbc(y)**2+Jac(x)**2-Jab(x)*Jbc(y)-Jbc(y)*Jac(rac)-Jab(x)*Jac(rac))
+
+
+      ELEPS_xyz=(Qab+Qbc+Qac-sqrt(Jab**2+Jbc**2+Jac**2-Jab*Jbc-Jbc*Jac-Jab*Jac))/autoeV
+
+      return
+    end function ELEPS_xyz
Index: src/Space_all.f90
===================================================================
--- src/Space_all.f90	(revision 2)
+++ src/Space_all.f90	(revision 3)
@@ -14,7 +14,8 @@
 !      INCLUDE 'SIZES'
 
  !MSET=mth iteration, MRESET=Maximum no. of iterations. 
-      INTEGER(KINT) :: NCoord,MRESET,MSET(NGeomF),NGeomF,IGeom
+      INTEGER(KINT) :: NCoord,MRESET,NGeomF,IGeom
+      INTEGER(KINT) :: MSET(NGeomF)
  !Geom=geometry with NCoord coordinates. GRAD=gradients.
       REAL(KREAL) :: Geom(NCoord),GRAD(NCoord)
       REAL(KREAL) :: Heat
