Index: egrad_chamfre.f90
===================================================================
--- egrad_chamfre.f90	(revision 3)
+++ egrad_chamfre.f90	(revision 4)
@@ -1,4 +1,4 @@
-!*************############### Attention	###############*************
+!*************############### Attention  ###############*************
 !(a) Two parameters aa and br(3,3) must be provided in main program
 !(b) At first, z=0 for the outmost surface 
 !(c) The units of all parameters are not atomic units, but angstrom, eV etc.
@@ -173,9 +173,9 @@
 !
 !rtemp=rlong
     if (ka(i+mid)<np+1.and.ka(j+mid)<np+1)then     ! Ni-Ni
-	     
-	      rtemp=rnn
-	 endif
+       
+        rtemp=rnn
+   endif
     if (ka(i+mid)<np+1.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then  ! Ni-H
           
            rtemp=rnh
@@ -189,9 +189,9 @@
             rtemp=rch
      endif
     if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then  ! H-H
-	       
-	        rtemp=rhh
-	 endif
+         
+          rtemp=rhh
+   endif
     if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)<np+1) then  ! H-Ni
             
              rtemp=rnh
@@ -228,10 +228,10 @@
              rij=sqrt((xj(1)-xi(1))**2.0+(xj(2)-xi(2))**2.0+(xj(3)-xi(3))**2.0)
              n1=ka(i)/(np+1);n2=ka(mf(i,j))/(np+1)    !!!  Ni 46  H  50  C
              n3=ka(i)/(np+5);n4=ka(mf(i,j))/(np+5)
-	
+  
  !**********************************Ni-Ni interaction
  ! **********
-         if(n1==0.and.n3==0.and.n2==0.and.n4==0)then		 !Ni-Ni
+         if(n1==0.and.n3==0.and.n2==0.and.n4==0)then     !Ni-Ni
            call fctNiNi(rnn,rij,fcut)
            call dfcutNiNi(rnn,rij,tf)
             eb=eb+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
@@ -239,26 +239,36 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut*(-ann(2)/ann(3))*(xi(im)-xj(im))/rij
-dv(i,2,im)=dv(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut*(-2.0*ann(1)/ann(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut*(-ann(2)/ann(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut*(-2.0*ann(1)/ann(3))*(xj(im)-xi(im))/rij
+   dv(i,1,im)=dv(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut &
+        *(-ann(2)/ann(3))*(xi(im)-xj(im))/rij
+   dv(i,2,im)=dv(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
+   dv(i,3,im)=dv(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut &
+        *(-2.0*ann(1)/ann(3))*(xi(im)-xj(im))/rij
+   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) &
+     *(rij/ann(3)-1.0))*fcut*(-ann(2)/ann(3))*(xj(im)-xi(im))/rij
+   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) &
+     *(rij/ann(3)-1.0))*fcut
+   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) &
+        *(rij/ann(3)-1.0))*fcut*(-2.0*ann(1)/ann(3))*(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(i,2,im)=df(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(i,1,im)=df(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(i,2,im)=df(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
+   df(i,3,im)=df(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) &
+     *(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) &
+        *(rij/ann(3)-1.0))*fcut
+   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) &
+        *(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
-             endif
+endif
 
  !**********************************Ni-H interaction
 ! **********
-         if(n1==0.and.n3==0.and.n2==1.and.n4==0)then		 !Ni-H
+         if(n1==0.and.n3==0.and.n2==1.and.n4==0)then     !Ni-H
            call fctNiH(rnh,rij,fcut)
            call dfcutNiH(rnh,rij,tf)
              eb=eb+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
@@ -266,19 +276,31 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
-dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xj(im)-xi(im))/rij
+   dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut &
+        *(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
+   dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
+   dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut &
+        *(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
+   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)* &
+        exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3)) &
+        *(xj(im)-xi(im))/rij
+   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)* &
+        exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
+   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)* &
+        exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3)) &
+        *(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
+   df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) &
+        *(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1) &
+        *(rij/anh(3)-1.0))*fcut
+   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1) &
+        *(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
 !
@@ -286,7 +308,7 @@
 
 !**********************************Ni-C interaction
 ! **********
-         if(n1==0.and.n3==0.and.n2==1.and.n4==1)then		 !Ni-C
+         if(n1==0.and.n3==0.and.n2==1.and.n4==1)then     !Ni-C
            call fctNiC(rnc,rij,fcut)
            call dfcutNiC(rnc,rij,tf)
             eb=eb+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
@@ -294,19 +316,29 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
-dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
-
-df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut &
+        *(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
+   dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
+   dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut &
+        *(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
+   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
+        anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
+   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut
+   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
+   
+   df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
+   df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
+        anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut
+   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
                endif
@@ -320,19 +352,29 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
-dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xj(im)-xi(im))/rij
+   dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut &
+        *(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
+   dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
+   dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut &
+        *(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
+   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) &
+     *(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xj(im)-xi(im))/rij
+   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* & 
+        anh(1)*(rij/anh(3)-1.0))*fcut
+   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* &
+        anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
+   df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0* &
+        anh(2)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* &
+        anh(1)*(rij/anh(3)-1.0))*fcut
+   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* &
+        anh(1)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
               endif
@@ -346,19 +388,29 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut*(-ahh(2)/ahh(3))*(xi(im)-xj(im))/rij
-dv(i,2,im)=dv(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut*(-2.0*ahh(1)/ahh(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut*(-ahh(2)/ahh(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut*(-2.0*ahh(1)/ahh(3))*(xj(im)-xi(im))/rij
-
-df(i,1,im)=df(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(i,2,im)=df(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   dv(i,1,im)=dv(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut &
+        *(-ahh(2)/ahh(3))*(xi(im)-xj(im))/rij
+   dv(i,2,im)=dv(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
+   dv(i,3,im)=dv(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut &
+       *(-2.0*ahh(1)/ahh(3))*(xi(im)-xj(im))/rij
+   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* &
+       ahh(2)*(rij/ahh(3)-1.0))*fcut*(-ahh(2)/ahh(3))*(xj(im)-xi(im))/rij
+   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* &
+       ahh(1)*(rij/ahh(3)-1.0))*fcut
+   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* &
+       ahh(1)*(rij/ahh(3)-1.0))*fcut*(-2.0*ahh(1)/ahh(3))*(xj(im)-xi(im))/rij
+   
+   df(i,1,im)=df(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(i,2,im)=df(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
+   df(i,3,im)=df(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* &
+        ahh(2)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
+   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* &
+        ahh(1)*(rij/ahh(3)-1.0))*fcut
+   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* &
+        ahh(1)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
              endif
@@ -374,19 +426,29 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
+dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut &
+        *(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
 dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
+dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut &
+        *(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
+dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
+        ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
+dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
+        ach(1)*(rij/ach(3)-1.0))*fcut
+dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
+        ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf*(xi(im)-xj(im))/rij
+df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
 df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
+        ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
+        ach(1)*(rij/ach(3)-1.0))*fcut
+df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
+        ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
                 endif
@@ -400,19 +462,29 @@
 
 !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
+dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut &
+        *(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
 dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
+dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut &
+        *(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
+dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
+        anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
+dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut
+dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf*(xi(im)-xj(im))/rij
+df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
 df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf &
+        *(xi(im)-xj(im))/rij
+df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
+        anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*fcut
+df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
+        anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
 
@@ -428,19 +500,29 @@
 
    !Calculation of FORCE
 do im=1,3
-dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
+dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut &
+       *(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
 dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
-dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
-dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
+dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut &
+       *(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
+dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
+       ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
+dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
+       ach(1)*(rij/ach(3)-1.0))*fcut
+dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
+       ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
 
-df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf*(xi(im)-xj(im))/rij
+df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf &
+       *(xi(im)-xj(im))/rij
 df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf*(xi(im)-xj(im))/rij
-df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
-df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
-df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf &
+       *(xi(im)-xj(im))/rij
+df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
+       ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
+df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
+       ach(1)*(rij/ach(3)-1.0))*fcut
+df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
+       ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
 enddo
 !
 
