!*************############### 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.
!(d) Before using this routine, you must confirm that the energy is CONSISTENT with its derivative. 
!                        This check of consistency is VERY IMPORTANT for energy conservation in MD.
!!!
!!  xa, ya,za present the positions of atoms. v is the energy of the system. dvdr present the force along the x,y,z directions.
!!
subroutine egrad_chamfre(Nat,V,Geom,Grad)
! WAS xa,ya,za,v,dvdr)


  use Path_module, only : Lat_a,Lat_b,Lat_c,order,orderinv,renum
  use Io_module


   IMPLICIT NONE

       INTEGER(KINT), INTENT(IN) :: Nat
       REAL(KREAL), INTENT(OUT) :: V,grad(Nat*3)
       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)

!     Bohr --> Angstr
      real(KREAL), parameter :: BOHR   = 0.52917726D+00

!! PFL 2010 Nov. 28 ->
! for now I keep this old F77 scheme. This has to be replaced by ALLOCATABLE
! things
! <- PFL 2010 Nov 28.
  REAL(KREAL) :: xa(100),ya(100),za(100),dvdr(100,3)
  REAL(KREAL) :: rshort,rlong,zshort,zlong
  REAL(KREAL), SAVE :: rnc,rnn,rhh,rch,rnh,ann(100),anc(100),anh(100),ach(100),ahh(100),a(100)
  REAL(KREAL) :: dv(100,10,3),df(100,10,3)
  REAL(KREAL) :: px(9000),py(9000),pz(9000),tmp,temp,rtemp,yf
  INTEGER(KINT) :: ka(9000),kx(9000),ky(9000),itx,ity,nf(100),mf(100,500)
  REAL(KREAL) :: xi(3),xj(3),xk(3),rij,rik,rjk,rji,rki,er,eb
  REAL(KREAL) :: fcut,tf,bodij,bodji,bod,cost,gtheta,dgdt,xt,t1,t2,t3
  REAL(KREAL) :: zi,fcut2
  INTEGER(KINT) :: np,nq,nh,nc
  INTEGER(KINT) :: i, j, nn, m, n, mid, n1, im, jm, n2,n3,n4, iat
  LOGICAL, SAVE :: First=.TRUE.
  LOGICAL :: debug, FExist

  common/zdis/zshort,zlong
  common/rdistanc/rshort,rlong


  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID
     
  END INTERFACE

!**********************************************


  debug=valid('egrad_chamfre')

 if (debug) Call Header('Entering Egrad_chamfre')

  v=0.0;dvdr=0.0
!******prepare the parameters of REBO potential
!rshort=3.5 
!rlong=3.9
zshort=3.5
zlong=4.5

rnn=3.2    !! to adjust
rnc=3.5    !! to adjust
rnh=3.0    !! to adjust
rch=3.2    !! to adjust
rhh=3.2    !! to adjust

!rph=rlong

 if (First) THEN
    First=.FALSE.
    INQUIRE(FILE="parameter.dat",EXIST=FExist)
    if (.NOT.Fexist) THEN
       WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"
       WRITE(*,*) "!!   egrad_chamfre             !!!!"
       WRITE(*,*) "The file parameter.dat does not exists."
       WRITE(*,*) "!!   egrad_chamfre             !!!!"
       WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"
       STOP
    END IF
open(1,file="parameter.dat")
do i=1,5   !Ni-Ni
read(1,*)ann(i)
enddo

do i=1,5   !Ni-C
read(1,*)anc(i)
enddo

do i=1,5   !Ni-H
read(1,*)anh(i)
enddo

do i=1,5   !H-H
read(1,*)ahh(i)
enddo
do i=1,5  !C-H
read(1,*)ach(i)
enddo
!********LONG INTERACTION **********
!do i=1,2
!read(1,*)a(i)
!enddo
close(1)
END IF
!************************************
np=45
nh=4
nc=1
nq=np+nh+nc
!************************************

 DO I=1,Nat
    Iat=Order(I)
    xa(I)=Geom(Iat,1)
    ya(I)=Geom(Iat,2)
    za(I)=Geom(Iat,3)
 END DO

px=0.0;py=0.0;pz=0.0
nn=0
 do m=-1,1
    do n=-1,1
       do i=1,nq
          if(m==0.and.n==0.and.i==1)mid=nn
            nn=nn+1
            px(nn)=xa(i)+m*lat_a(1)+n*lat_b(1)
            py(nn)=ya(i)+m*lat_a(2)+n*lat_b(2)
            pz(nn)=za(i)
            ka(nn)=i;kx(nn)=m;ky(nn)=n
       enddo
    enddo
 enddo
!***********************xyz file required for MS input
 open(123,file="111.xyz",ACCESS='APPEND')
    write(123,*)nn
    write(123,*)
   do i=1,nn
       if(ka(i)<(np+1))then
            write(123,10)'Ni',px(i),py(i),pz(i)
       elseif (ka(i)<np+5) then
            write(123,10)'H',px(i),py(i),pz(i)
       else
            write(123,10)'C',px(i),py(i),pz(i)
  10 format(a6,3f15.8)
       endif
   enddo
 close(123)
!******************************************
 do i=1,nq
    tmp=px(i);px(i)=px(i+mid);px(i+mid)=tmp
    tmp=py(i);py(i)=py(i+mid);py(i+mid)=tmp
    tmp=pz(i);pz(i)=pz(i+mid);pz(i+mid)=tmp
    itx=kx(i);kx(i)=kx(i+mid);kx(i+mid)=itx
    ity=ky(i);ky(i)=ky(i+mid);ky(i+mid)=ity
 enddo
 mid=0
!******************************************
  nf=0;mf=0
    do i=1,nq
        nf(i)=0
           do j=1,nn
             temp=sqrt((px(i)-px(j))**2+(py(i)-py(j))**2+(pz(i)-pz(j))**2)
!
!rtemp=rlong
    if (ka(i+mid)<np+1.and.ka(j+mid)<np+1)then     ! Ni-Ni
	     
	      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
     endif
    if (ka(i+mid)<np+1.and.ka(j+mid)>np+4) then   ! Ni-C
           
           rtemp=rnc
     endif
      if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np+4)  then ! H-C
           
            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
    if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)<np+1) then  ! H-Ni
            
             rtemp=rnh
       endif
    if (ka(i+mid)>np+4.and.ka(j+mid)<np+1) then !  C-Ni
            
             rtemp=rnc
     endif
    if (ka(i+mid)>np+4.and.ka(j+mid)>np.and.ka(j+mid)<np+5) then ! C-H
             
              rtemp=rch
     endif

      if(temp>0.1.and.temp<rtemp)then
         nf(i)=nf(i)+1
         mf(i,nf(i))=j
      endif
!
    enddo
  !write(*,*)nf(i)
 enddo
!**************************************
!stop

!*************************************Starting to calculate the sum
  yf=0.0;dvdr=0.0
   do i=1,nq
      er=0.0;eb=0.0
      dv=0.0;df=0.0
      do j=1,nf(i)

             xi(1)=px(i);      xi(2)=py(i);      xi(3)=pz(i)
             xj(1)=px(mf(i,j));xj(2)=py(mf(i,j));xj(3)=pz(mf(i,j))
             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
           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
            er=er+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut

!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

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

 !**********************************Ni-H interaction
! **********
         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
             er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut

!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

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

!**********************************Ni-C interaction
! **********
         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
            er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut

!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
enddo
!
               endif
!!**********************************H-Ni interaction
! **********
         if(n1==1.and.n3==0.and.n2==0.and.n4==0)then        ! H-Ni
           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
             er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut

!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

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
!!**********************************H-H interaction
! **********
        if(n1==1.and.n3==0.and.n2==1.and.n4==0)then   !H-H
           call fctHH(rhh,rij,fcut)
           call dfcutHH(rhh,rij,tf)
              eb=eb+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
              er=er+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut

!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
enddo
!
             endif


!!**********************************H-C interaction
! **********
        if(n1==1.and.n3==0.and.n2==1.and.n4==1)then   !H-C
            call fctCH(rch,rij,fcut)
            call dfcutCH(rch,rij,tf)
             eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
             er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut

!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,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

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
enddo
!
                endif
!!**********************************C-Ni interaction
! **********
        if(n1==1.and.n3==1.and.n2==0.and.n4==0)then   !C-Ni
            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
             er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut

!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
enddo
!

               endif
!!**********************************C-H interaction
! **********
       if(n1==1.and.n3==1.and.n2==1.and.n4==0)then   !C-H
           call fctCH(rch,rij,fcut)
           call dfcutCH(rch,rij,tf)
             eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
             er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut


   !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,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

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
enddo
!

        endif
  !!***************************ONE end****************      

enddo
!!******** CALCULATION ENERGY
!
        yf=yf+er-sqrt(eb)
!!*********
       do im=1,3
            dvdr(i,im)=dvdr(i,im)+dv(i,1,im)-0.5/(dv(i,2,im))**0.5*dv(i,3,im)
            dvdr(i,im)=dvdr(i,im)+df(i,1,im)-0.5/(df(i,2,im))**0.5*df(i,3,im)
         do j=1,nf(i+mid)
             jm=ka(mf(i+mid,j))
             dvdr(jm,im)=dvdr(jm,im)+dv(jm,1,im)-0.5/(dv(i,2,im))**0.5*dv(jm,3,im)
             dvdr(jm,im)=dvdr(jm,im)+df(jm,1,im)-0.5/(df(i,2,im))**0.5*df(jm,3,im)
         enddo
       enddo
!
enddo
!!*****************************TWO ENDS**************************
!***************************************Long-ranged interaction
!do i=nq-4,nq
!     zi=pz(i+mid) !-0.3*aa*br(3,3)
!     call fct2(zshort,zlong,zi,fcut2)
 !    call dfcut2(zshort,zlong,zi,tf)
 !    yf=yf+fcut2*(a(6)-a(7)/zi**2.0)

 !  dvdr(i,3)=dvdr(i,3)+fcut2*(a(7)*2.0/zi**3.0)+tf*(a(6)-a(7)/zi**2.0)

!enddo
!*******************************************************************
  v=yf 
   do im=1,3
      do jm=1,nq
          dvdr(jm,im)=dvdr(jm,im)  
      enddo
  enddo

!PFL: We convert dvdr into Grad
      do jm=1,nq
          Iat=Jm
          if (renum) Iat=order(jm)
          Grad(3*Iat-2:3*Iat)=dvdr(jm,1:3)
      enddo
  
 if (debug) WRITE(*,*) "EGRAD_CHAMFRE, E=",V
  
 if (debug) call Header('Egrad_chamfre OVER')

!********************************************ENDENDENDENDENDENDENDEND
return
end subroutine egrad_chamfre

!******************************************
!**********PARAMETERS**********************
subroutine fctNiNi(r2,r,fcut)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)

 REAL(KREAL) :: r2,r,fcut,r1,pi
pi=4.0*atan(1.0)
r1=r2-0.5

if(r<=r1)then
fcut=1.0
endif
if(r>r1.and.r<=r2)then
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
fcut=0.0
endif
return
end subroutine fctNiNi
!!************
subroutine dfcutNiNi(r2,r,df)
   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)
  REAL(KREAL) :: r, r2,df,pi,r1
pi=4.0*atan(1.0)
r1=r2-0.5

if(r<=r1)then
df=0.0
endif
if(r>r1.and.r<=r2)then
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
df=0.0
endif
return
end subroutine dfcutNiNi


!!************
subroutine fctNiH(r2,r,fcut)
   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)

 REAL(KREAL) ::r2,r,fcut
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=2.7
 rlong=3.0
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 fcut=1.0
endif
if(r>r1.and.r<=r2)then
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
 fcut=0.0
endif
return
end subroutine fctNiH
!!**********

subroutine dfcutNiH(r2,r,df)


   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::r2,r,df
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=2.7
 rlong=3.0
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 df=0.0
endif
if(r>r1.and.r<=r2)then
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
 df=0.0
endif
return
end subroutine dfcutNiH

!!************
subroutine fctNiC(r2,r,fcut)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::r2,r,fcut
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=3.2
 rlong=3.5
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 fcut=1.0
endif
if(r>r1.and.r<=r2)then
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
 fcut=0.0
endif
return
end subroutine fctNiC

!!***********
subroutine dfcutNiC(r2,r,df)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)

     
REAL(KREAL) ::r2,r,df
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=3.2
 rlong=3.5
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 df=0.0
endif
if(r>r1.and.r<=r2)then
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
 df=0.0
endif
return
end subroutine dfcutNiC

!!************
subroutine fctCH(r2,r,fcut)


   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::r2,r,fcut
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=2.1
 rlong=3.2
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 fcut=1.0
endif
if(r>r1.and.r<=r2)then
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
 fcut=0.0
endif
return
end subroutine fctCH
!!*****************

subroutine dfcutCH(r2,r,df)


   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)

     
REAL(KREAL) ::r2,r,df
REAL(KREAL) ::rshort,rlong,r1,pi

 rshort=2.1
 rlong=3.2
 pi=4.0*atan(1.0)
 r1=rshort
 r2=rlong
if(r<=r1)then
 df=0.0
endif
if(r>r1.and.r<=r2)then
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
 df=0.0
endif
return
end subroutine dfcutCH
!!************
subroutine fctHH(r2,r,fcut)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::r2,r,fcut,r1,pi
pi=4.0*atan(1.0)
r1=r2-1.1
if(r<=r1)then
 fcut=1.0
endif
if(r>r1.and.r<=r2)then
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
 fcut=0.0
endif
return
end subroutine fctHH
!!************
subroutine dfcutHH(r2,r,df)


   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::r2,r,df,r1,pi
pi=4.0*atan(1.0)
r1=r2-1.1
if(r<=r1)then
 df=0.0
endif
if(r>r1.and.r<=r2)then
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
 df=0.0
endif
return
end subroutine dfcutHH

!!*********************LONG INTERACTION

subroutine dfcut2(zshort,zlong,r,df)


   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::zshort,zlong,r,df,pi, r1, r2

pi=4.0*atan(1.0)

r1=zshort
r2=zlong

if(r<=r1)then
df=0.0
endif
if(r>r1.and.r<=r2)then
df=pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
endif
if(r>r2)then
df=0.0
endif

return
end subroutine dfcut2

subroutine fct2(zshort,zlong,r,fcut)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


REAL(KREAL) ::zshort,zlong,r,fcut,pi,r1,r2

pi=4.0*atan(1.0)

r1=zshort
r2=zlong

if(r<=r1)then
fcut=0.0
endif
if(r>r1.and.r<=r2)then
fcut=1.0-(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
endif
if(r>r2)then
fcut=1.0
endif

return
end subroutine fct2

