--- forces-edip.1.10.f	Fri Apr 19 10:44:21 2002
+++ forces-edip.1.11.f	Fri Apr 19 10:44:38 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.10f
+C   Version 1.11f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -109,89 +109,84 @@
 C    u4 = 0.66;
 
 C  ******************************************
- 
-      subroutine init_EDIP()
+
+      subroutine compute_forces_EDIP(N_own, pos, L_x, L_y, L_z, 
+     &           E_potential, f)
       implicit none
 
 C     NOTE: you do not need these header files, since you will have to
 C     customize the way that atoms are passed to the force subroutine 
 C     for your particular MD program.
-      include "edip_pot_include.h"
       include "edip_neighbors_include.h"
 
-      call input_edip_params()
-      par_bg=par_a
-      par_eta = par_delta/par_Qo
-      pot_cutoff = par_a
-      delta_safe = 0.2d0
-
-      end subroutine
-
-C ******************************************
-
-      subroutine input_EDIP_params()
-      implicit none
- 
-      include "edip_pot_include.h"
- 
-      print *, "EDIP-Si parameters:"
-C   taken from Justo et al., Phys. Rev. B 58, 2539 (1998).
-      par_cap_A = 5.6714030d0
-      par_cap_B = 2.0002804d0
-      par_rh = 1.2085196d0
-      par_a = 3.1213820d0
-      par_sig = 0.5774108d0
-      par_lam = 1.4533108d0
-      par_gam = 1.1247945d0
-      par_b = 3.1213820d0
-      par_c = 2.5609104d0
-      par_delta = 78.7590539d0
-      par_mu = -0.6966326d0
-      par_Qo = 312.1341346d0
-      par_palp = 1.4074424d0
-      par_bet = -0.0070975d0
-      par_alp = 3.1083847d0
-      print *, par_cap_A,par_cap_B,par_rh,par_a,par_sig
-      print *, par_lam,par_gam,par_b,par_c,par_delta
-      print *, par_mu,par_Qo,par_palp,par_bet,par_alp
-
-      end subroutine
+c#define FIXZ
+#undef FIXZ
+c   switch on/off option using fixed Z
 
-C  ******************************************
 
-      subroutine compute_forces_EDIP(N_own, pos, L_x, L_y, L_z, 
-     &           E_potential, f)
-      implicit none
+C  ------------------------- EDIP PARAMETERS ---------------------------
 
+      double precision par_cap_A,par_cap_B,par_rh,par_a,par_sig
+      double precision par_lam,par_gam,par_b,par_c,par_delta
+      double precision par_mu,par_Qo,par_palp,par_bet,par_alp
+      double precision u1,u2,u3,u4
+      double precision asqr,Qort,muhalf,u5,bmc,cmb3inv,par_bg,par_eta
+c      double precision pot_cutoff,delta_safe
 
-c#define FIXZ
-#undef FIXZ
-c   switch on/off option using fixed Z
+C   taken from Justo et al., Phys. Rev. B 58, 2539 (1998).
+      parameter (par_cap_A = 5.6714030d0)
+      parameter (par_cap_B = 2.0002804d0)
+      parameter (par_rh = 1.2085196d0)
+      parameter (par_a = 3.1213820d0)
+      parameter (par_sig = 0.5774108d0)
+      parameter (par_lam = 1.4533108d0)
+      parameter (par_gam = 1.1247945d0)
+      parameter (par_b = 3.1213820d0)
+      parameter (par_c = 2.5609104d0)
+      parameter (par_delta = 78.7590539d0)
+      parameter (par_mu = -0.6966326d0)
+      parameter (par_Qo = 312.1341346d0)
+      parameter (par_palp = 1.4074424d0)
+      parameter (par_bet = -0.0070975d0)
+      parameter (par_alp = 3.1083847d0)
+
+      parameter (u1 = -0.165799d0)
+      parameter (u2 = 32.557d0)
+      parameter (u3 = 0.286198d0)
+      parameter (u4 = 0.66d0)
+
+C   DERIVATIVE COEFFICIENTS
+
+      parameter (asqr = par_a*par_a)
+      parameter (Qort = 17.66731826283d0) ! = sqrt(par_Qo)
+      parameter (muhalf = par_mu*0.5D0)
+      parameter (u5 = u2*u4)
+      parameter (bmc = par_b-par_c)
+      parameter (cmb3inv = 3.0D0/(par_c-par_b))
+      parameter (par_bg=par_a)
+      parameter (par_eta = par_delta/par_Qo)
+c      parameter (pot_cutoff = par_a)
+c      parameter (delta_safe = 0.2d0)
 
 
 C  ------------------------- VARIABLE DECLARATIONS -------------------------
 
 
-      integer N_own
+      integer N_own,MAX_NBRS_1
       double precision pos(3,N_own)
       double precision E_potential
       double precision f(3,N_own)
       double precision L_x, L_y, L_z
 
-      include "edip_pot_include.h"
-      include "edip_neighbors_include.h"
-
       integer i,j,k,l,n
-      double precision dx,dy,dz,r,rsqr,asqr
+      double precision dx,dy,dz,r,rsqr
       double precision rinv,rmainv,Z
       double precision dV2j,dV2ijx,dV2ijy,dV2ijz,pZ,dp
       double precision temp0,temp1,temp2
-      double precision Qort,muhalf,u5
       double precision rmbinv,winv,dwinv,tau,dtau,lcos,x,H
       double precision winv_2lam,TdW_WdT_over_W
       double precision dV3rij,dV3rik,dV3l
       double precision dEdrl,dEdrlx,dEdrly,dEdrlz
-      double precision bmc,cmb3inv
       double precision fjx,fjy,fjz,fkx,fky,fkz
 
       double precision s2_t0(MAX_NBRS_1)
@@ -265,16 +260,6 @@
       V3=0.0d0
 
 
-C   COMBINE COEFFICIENTS
-
-      asqr = par_a*par_a
-      Qort = sqrt(par_Qo)
-      muhalf = par_mu*0.5D0
-      u5 = u2*u4
-      bmc = par_b-par_c
-      cmb3inv = 3.0D0/(par_c-par_b)
-
-
 C  --- LEVEL 1: OUTER LOOP OVER ATOMS ---
 
       do i=1, N_own