Chimie Théorique » OpenPath

!*************############### 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