C--------------------------------------------------------------------
C       PROGRAM DEFECT:
C       VERSION 2.0 (JUNE.6,1996):
C       1. USE TERSOFF POTENTIAL
C       2. USE CONJUGATE GRADIENT METHOD TO MINIMIZE THE ENERGY
C --------------------------------------------------------------------
C       USE THE FOLLOWING UNITS:                                 
C       ENERGY UNIT = ev  = 1.602E-19 JOULE
C       MASS   UNIT = AMU = 1.66054E-27 KG                      
C       LENGTH UNIT = ANGSTROM                                 
C       TIME   UNIT = 1.018043E-14 SECOND = 10.18043 FS
C       STRESS UNIT = 1.602192 MBAR
C --------------------------------------------------------------------
 
      PROGRAM DEFECT
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=50,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),LAST,
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM),N,MODE,VIRILSUM(3,3),NEW,
     A                HSUM(3,3),STRESSUM(3,3),STRESS(3,3),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2,R_SI_SI,S_SI_SI,
     A                R_C_C,S_C_C,R_SI_C,S_SI_C
      COMMON SEED

      LOGICAL WRITE_LP,UPDATE,NEW
      INTEGER PERTURBATION

      REAL XX(3*NM)
      DOUBLE PRECISION SX(NM),SY(NM),SZ(NM),S(NM),H(3,3),HH(3,3),DELTA
      DATA LP,LP_CONFIG  /25,26/
      CHARACTER NAME *20

      UENERGY_IN_JOULE = 1.602192E-19
      ULENGTH_IN_METER = 1.E-10
      UMASS_IN_KG = 1.66053E-27
      UTIME_IN_SECOND = 
     A   ULENGTH_IN_METER*SQRT(UMASS_IN_KG/UENERGY_IN_JOULE)
      UTIME_IN_PS = UTIME_IN_SECOND/1.E-12 
      UTIME_IN_FS = UTIME_IN_SECOND/1.E-15
      USTRESS_IN_PASCAL = UENERGY_IN_JOULE/ULENGTH_IN_METER**3
      USTRESS_IN_BAR  = USTRESS_IN_PASCAL/1.E5
      USTRESS_IN_MBAR = USTRESS_IN_PASCAL/1.E11
      AVO   = 6.0221367E+23
      BOLZ  = 1.380658E-23
      HBAR  = 1.054589E-34
      PM(1) = 28.086
      PM(2) = 12.011 

      PA(1)   = 1.8308E3
      PB(1)   = 4.7118E2
      PL1(1)  = 2.4799
      PL2(1)  = 1.7322
      PN(1)   = 7.8734E-1
      PBT(1)  = (1.1E-6)**PN(1)
      PC(1)   = (1.0039E5)**2.E0
      PD(1)   = (1.6217E1)**2.E0
      PH(1)   = -5.9825E-1
      PR(1)   = 2.70
      PS(1)   = 3.00
      C2D2(1) = 1.E0+PC(1)/PD(1)
      PWR(1)  = -1.E0/(2.E0*PN(1))
      PWR1(1) = PWR(1)-1.E0
      PWR2(1) = PN(1)-1.E0

C     SET OF PARAMTERS FOR C. 
      PA(2)   = 1.5448E3
      PB(2)   = 3.8963E2
      PL1(2)  = 3.4653
      PL2(2)  = 2.3064
      PN(2)   = 9.9054E-1
      PBT(2)  = (4.1612E-6)**PN(2)
      PC(2)   = (1.9981E4)**2.E0
      PD(2)   = (7.0340)**2.E0
      PH(2)   = -3.9953E-1
      PR(2)   = 1.80
      PS(2)   = 2.10
      C2D2(2) = 1.E0+PC(2)/PD(2)
      PWR(2)  = -1.E0/(2.E0*PN(2))
      PWR1(2) = PWR(2)-1.E0
      PWR2(2) = PN(2)-1.E0

      RCUT = 3.00
      RLIST = 1.1*RCUT
      SEED  = 0.2573

      OPEN (UNIT=LP_CONFIG, FORM='UNFORMATTED',FILE="config")
      UPDATE = .TRUE.
      NEW = .TRUE.

      READ *,LAST,CUBE,TOL,N,NPOINTER,IW,PERTURBATION
      READ *,R_SI_SI,S_SI_SI,R_C_C,S_C_C,R_SI_C,S_SI_C,CHI_SI_C
      NC = NINT((N/8.)**(1.E0/3.E0))
      NOLD = N

      CUBE = CUBE*NC
      CUBE2 = CUBE*0.5
      FCUBE2 = -CUBE2
      IF (LAST.EQ.1) THEN 
      IW  = 1
      WRITE (6,468) R_SI_SI,S_SI_SI,R_C_C,S_C_C,R_SI_C,S_SI_C, 
     A              CHI_SI_C
 468  FORMAT ('R(Si-Si) = ',F4.2,'  S(Si-Si) = ',F4.2,
     A        '  R(C-C) = ',F4.2,'  S(C-C) = ',F4.2,
     A        '  R(Si-C) = ',F4.2,'  S(Si-C) = ',F4.2,
     A        '  CHI(Si-C) = ',F6.4)
      ENDIF

 154  IF (LAST.EQ.2) THEN 
       READ (LP_CONFIG) N,H,HH,HH,HH,HH,HH
       READ (LP_CONFIG) SX,SY,SZ
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) NTYPE,DELTA
       CLOSE (UNIT = LP_CONFIG)
       OPEN (UNIT=LP_CONFIG, FORM='UNFORMATTED',FILE="config")
       epp = -6.386389
      ELSE
C     --------------------------------------------------
C     DIAMOND CUBIC STRUCTURE
       N = NOLD
       A = 1.E0/FLOAT(NC)
       X(1) = 0.E0
       Y(1) = 0.E0
       Z(1) = 0.E0
       X(3) = 1.E0/2.E0
       Y(3) = 1.E0/2.E0
       Z(3) = 0.E0 
       X(5) = 0.E0
       Y(5) = 1.E0/2.E0
       Z(5) = 1.E0/2.E0 
       X(7) = 1.E0/2.E0
       Y(7) = 0.E0
       Z(7) = 1.E0/2.E0

       IP = 0
       DO 217 I=1,NC
       DO 217 J=1,NC
       DO 217 K=1,NC
       DO 218 L=1,7,2
       SX(IP+L) = (X(L)+I-7.E0/8.E0)*A
       SY(IP+L) = (Y(L)+J-7.E0/8.E0)*A
       SZ(IP+L) = (Z(L)+K-7.E0/8.E0)*A
       NTYPE (IP+L) = 1
       SX(IP+L+1) = (X(L)+I-5.E0/8.E0)*A
       SY(IP+L+1) = (Y(L)+J-5.E0/8.E0)*A
       SZ(IP+L+1) = (Z(L)+K-5.E0/8.E0)*A
       NTYPE (IP+L+1) = 2
 218   CONTINUE
       IP = IP+8
 217   CONTINUE


       if (.not.((last.eq.1).and.(iw.ne.1))) then 
       IP = 1
       DO 52 I=1,N
       X(I) = CUBE*SX(I)
       Y(I) = CUBE*SY(I)
       Z(I) = CUBE*SZ(I)
       XOLD(I) = X(I)
       YOLD(I) = Y(I)
       ZOLD(I) = Z(I)
       DX(I) = 0.
       DY(I) = 0.
       DZ(I) = 0.
       XX(IP)   = X(I)
       XX(IP+1) = Y(I)
       XX(IP+2) = Z(I)
       IP = IP+3
 52    CONTINUE
       UPDATE = .TRUE.
       NEW = .TRUE.
       call tersoff_value
       epp = pote/n
       WRITE (6,408) epp,-6.4325
 408   FORMAT(' epp           ',F7.3,' eV',F9.3,' eV')
       endif

       IF (IW.EQ.1) THEN
c      carbon vacancy 
       sx(2) = sx(n)
       sy(2) = sy(n)
       sz(2) = sz(n)
       ntype(2) = ntype(n)
       n = n-1
       name = " Vc "
       eab  = 5.9
c      ab initio result

       ELSEIF (IW.EQ.2) THEN 
c      silicon vacancy
       sx(1) = sx(n)
       sy(1) = sy(n)
       sz(1) = sz(n)
       ntype(1) = ntype(n)
       n = n-1
       name = " VSi "
       eab  = 6.8

       ELSEIF (IW.EQ.3) THEN 
c      carbon antisite
       ntype(1) = 2
       name = " C_si "
       eab  = 1.1

       ELSEIF (IW.EQ.4) THEN 
c      silicon antisite
       ntype(2) = 1
       name = " Si_c "
       eab  = 7.3

       ELSEIF (IW.EQ.5) THEN 
c      interstitial C at Tc site
       n = n + 1
       ntype(n) = 2
       sx(n) = sx(1)+A/2.
       sy(n) = sy(1)+A/2.
       sz(n) = sz(1)+A/2.
       name = " CTc "
       eab  = 11.0

       ELSEIF (IW.EQ.6) THEN 
c      interstitial C at TSi site
       n = n + 1
       ntype(n) = 2
       sx(n) = sx(2)+A/2.
       sy(n) = sy(2)+A/2.
       sz(n) = sz(2)+A/2.
       name = " CTsi "
       eab  = 8.6

       ELSEIF (IW.EQ.7) THEN 
c      interstitial Si at TSi site
       n = n + 1
       ntype(n) = 1
       sx(n) = sx(2)+A/2.
       sy(n) = sy(2)+A/2.
       sz(n) = sz(2)+A/2.
       name = " SiTsi "
       eab  = 15.0

       ELSEIF (IW.EQ.8) THEN 
c      interstitial Si at Tc site
       n = n + 1
       ntype(n) = 1
       sx(n) = sx(1)+A/2.
       sy(n) = sy(1)+A/2.
       sz(n) = sz(1)+A/2.
       name = " SiTc "
       eab  = 14.7

       ELSEIF (IW.EQ.9) THEN 
c      nearest neighbour antisite pair
       ntype(1) = 2
       ntype(2) = 1
       name = " antisite pair "
       eab = 5.9

       ELSEIF (IW.EQ.10) THEN 
c      di-vacancy 
       sx(1) = sx(n-1)
       sy(1) = sy(n-1)
       sz(1) = sz(n-1)
       ntype(1) = ntype(n-1)
       sx(2) = sx(n)
       sy(2) = sy(n)
       sz(2) = sz(n)
       ntype(2) = ntype(n)
       n = n-2
       name = " divacancy     "
       eab = 8.1
       ENDIF
C     --------------------------------------------------
       ENDIF


       IF (PERTURBATION.EQ.1) THEN 
C       test: a little perturbation
        do 543 i=1,n
        call random
        sz(i) = sz(i)+A*0.001*(seed-0.5)
 543    continue
       ENDIF

       IP = 1
       DO 57 I=1,N
       X(I) = CUBE*SX(I)
       Y(I) = CUBE*SY(I)
       Z(I) = CUBE*SZ(I)
       XOLD(I) = X(I)
       YOLD(I) = Y(I)
       ZOLD(I) = Z(I)
       DX(I) = 0.
       DY(I) = 0.
       DZ(I) = 0.
       XX(IP)   = X(I)
       XX(IP+1) = Y(I)
       XX(IP+2) = Z(I)
       IP = IP+3
 57    CONTINUE
       UPDATE = .TRUE.
       NEW = .TRUE.
       call tersoff_value

       IF (LAST.NE.1) THEN 
       print *, pote, pote/n,n
       print *,' '

      if (npointer.eq.0) goto 903
C     print out the defect coordinates
      ip = npointer
      if (ntype(ip).eq.1) print *,'  The indicated atom is a Si'
      if (ntype(ip).eq.2) print *,'  The indicated atom is a C'
      do 54 i=1,index(ip)
      jpt = list(ip,i)
      rx = x(jpt)-x(ip)
      ry = y(jpt)-y(ip)
      rz = z(jpt)-z(ip)
      if (rx.lt.fcube2) rx = rx+cube
      if (rx.gt.cube2)  rx = rx-cube
      if (ry.lt.fcube2) ry = ry+cube
      if (ry.gt.cube2)  ry = ry-cube
      if (rz.lt.fcube2) rz = rz+cube
      if (rz.gt.cube2)  rz = rz-cube
      r = sqrt(rx*rx+ry*ry+rz*rz)
      if (ntype(jpt).eq.1) print *,jpt,'   Si ',r
      if (ntype(jpt).eq.2) print *,jpt,'    C ',r
 54   continue
      print *, ' '
      ENDIF

 903  CALL FRPRMN(XX,3*N,TOL,ITER,FMIN,VALUE,FORCE)
      
      IF (LAST.NE.1) THEN 
      PRINT *,'THE MININUM CONVERGE TO', FMIN,' eV'
      PRINT *,'AFTER ',ITER,' ITERATIONS.'
      PRINT *,'THE DEFECT FORMATION ENERGY IS',
     A        FMIN-N*EPP,' eV for',name

      if (npointer.ne.0) then 
C     print out the defect coordinates
      ip = npointer
      print *,' '
      if (ntype(ip).eq.1) print *,'  The indicated atom is a Si'
      if (ntype(ip).eq.2) print *,'  The indicated atom is a C'
      do 558 i=1,index(ip)
      jpt = list(ip,i)
      rx = x(jpt)-x(ip)
      ry = y(jpt)-y(ip)
      rz = z(jpt)-z(ip)
      if (rx.lt.fcube2) rx = rx+cube
      if (rx.gt.cube2)  rx = rx-cube
      if (ry.lt.fcube2) ry = ry+cube
      if (ry.gt.cube2)  ry = ry-cube
      if (rz.lt.fcube2) rz = rz+cube
      if (rz.gt.cube2)  rz = rz-cube
      r = sqrt(rx*rx+ry*ry+rz*rz)
      if (ntype(jpt).eq.1) print *,jpt,'   Si ',r
      if (ntype(jpt).eq.2) print *,jpt,'    C ',r
 558  continue
      print *, ' '
      endif

      ELSE
      WRITE (6,438) name,FMIN-N*epp,eab
 438  FORMAT(A15,F7.3,' eV',F9.3,' eV')
      ENDIF

      IF (LAST.EQ.1.AND.IW.NE.10) THEN 
      IW = IW+1
      GOTO 154
      ENDIF

C     WRITE CONFIG FILE
      DO 987 I=1,3
      DO 987 J=1,3
      IF (I.EQ.J) THEN
      H(I,J) = CUBE
      ELSE 
      H(I,J) =0.
      ENDIF
      HH(I,J) = 0.
 987  CONTINUE
      DO 46 I=1,NM
 46   S(I) = 0.
      DO 47 I=1,N
      SX(I) = X(I)/CUBE
      SY(I) = Y(I)/CUBE
      SZ(I) = Z(I)/CUBE
      IF(SX(I).LT.0.D0) SX(I)=SX(I)+1.D0
      IF(SX(I).GT.1.D0) SX(I)=SX(I)-1.D0
      IF(SY(I).LT.0.D0) SY(I)=SY(I)+1.D0
      IF(SY(I).GT.1.D0) SY(I)=SY(I)-1.D0
      IF(SZ(I).LT.0.D0) SZ(I)=SZ(I)+1.D0
      IF(SZ(I).GT.1.D0) SZ(I)=SZ(I)-1.D0
 47   CONTINUE
      DELTA = 0.049
      WRITE (LP_CONFIG) N,H,HH,HH,HH,HH,HH
      WRITE (LP_CONFIG) SX,SY,SZ
      WRITE (LP_CONFIG) S,S,S
      WRITE (LP_CONFIG) S,S,S
      WRITE (LP_CONFIG) S,S,S
      WRITE (LP_CONFIG) S,S,S
      WRITE (LP_CONFIG) S,S,S
      WRITE (LP_CONFIG) NTYPE,DELTA

C     print out the defect coordinates
c      do 54 i=1,index(1)
c      jpt = list(105,i)
c 54   print *,jpt,x(jpt)-x(105),y(jpt)-y(105),
c     a        z(jpt)-z(105)
c      print *,' '
c      do 55 i=1,index(n)
c      jpt = list(n,i)
c 55   print *,jpt,x(jpt)-x(105),y(jpt)-y(105),
c     a        z(jpt)-z(105)

      STOP
      END


      FUNCTION VALUE(XX)
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=50,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),LAST,
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM),N,MODE,VIRILSUM(3,3),NEW,
     A                HSUM(3,3),STRESSUM(3,3),STRESS(3,3),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      LOGICAL WRITE_LP,UPDATE,NEW
      REAL VALUE,XX(3*NM)
      EXTERNAL TERSOFF_VALUE
      
      IP = 1
      DO 465 I=1,N
      X(I) = XX(IP) 
      Y(I) = XX(IP+1)
      Z(I) = XX(IP+2)
      IP = IP+3
 465  CONTINUE

      CALL TERSOFF_VALUE
      VALUE = POTE
      RETURN
      END


      FUNCTION FORCE(XX,F)
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=50,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),LAST,
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM),N,MODE,VIRILSUM(3,3),NEW,
     A                HSUM(3,3),STRESSUM(3,3),STRESS(3,3),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      LOGICAL WRITE_LP
      REAL XX(3*NM),F(3*NM)

      IP = 1
      DO 465 I=1,N
      X(I) = XX(IP) 
      Y(I) = XX(IP+1)
      Z(I) = XX(IP+2)
      IP = IP+3
 465  CONTINUE

      CALL TERSOFF_FORCE

      IP = 1
      DO 462 I=1,N
      F(IP)   = FX(I) 
      F(IP+1) = FY(I)
      F(IP+2) = FZ(I)
      IP = IP+3
 462  CONTINUE
    
      RETURN
      END



      SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func)
      REAL ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY
      EXTERNAL func
      PARAMETER (GOLD=1.618034, GLIMIT=1000., TINY=1.e-20)
      REAL dum,fu,q,r,u,ulim
      fa=func(ax)
      fb=func(bx)
      if(fb.gt.fa)then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif
      cx=bx+GOLD*(bx-ax)
      fc=func(cx)
1     if(fb.ge.fc)then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r))
        ulim=bx+GLIMIT*(cx-bx)
        if((bx-u)*(u-cx).gt.0.)then
          fu=func(u)
          if(fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            return
          else if(fu.gt.fb)then
            cx=u
            fc=fu
            return
          endif
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        else if((cx-u)*(u-ulim).gt.0.)then
          fu=func(u)
          if(fu.lt.fc)then
            bx=cx
            cx=u
            u=cx+GOLD*(cx-bx)
            fb=fc
            fc=fu
            fu=func(u)
          endif
        else if((u-ulim)*(ulim-cx).ge.0.)then
          u=ulim
          fu=func(u)
        else
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        goto 1
      endif
      return
      END

      FUNCTION brent(ax,bx,cx,f,tol,xmin)
      INTEGER ITMAX
      REAL brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS
      EXTERNAL f
      PARAMETER (ITMAX=1000,CGOLD=.3819660,ZEPS=1.0e-10)
      INTEGER iter
      REAL a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
      a=min(ax,cx)
      b=max(ax,cx)
      v=bx
      w=v
      x=v
      e=0.
      fx=f(x)
      fv=fx
      fw=fx
      do 11 iter=1,ITMAX
        xm=0.5*(a+b)
        tol1=tol*abs(x)+ZEPS
        tol2=2.*tol1
        if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3
        if(abs(e).gt.tol1) then
          r=(x-w)*(fx-fv)
          q=(x-v)*(fx-fw)
          p=(x-v)*q-(x-w)*r
          q=2.*(q-r)
          if(q.gt.0.) p=-p
          q=abs(q)
          etemp=e
          e=d
          if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) 
     *goto 1
          d=p/q
          u=x+d
          if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x)
          goto 2
        endif
1       if(x.ge.xm) then
          e=a-x
        else
          e=b-x
        endif
        d=CGOLD*e
2       if(abs(d).ge.tol1) then
          u=x+d
        else
          u=x+sign(tol1,d)
        endif
        fu=f(u)
        if(fu.le.fx) then
          if(u.ge.x) then
            a=x
          else
            b=x
          endif
          v=w
          fv=fw
          w=x
          fw=fx
          x=u
          fx=fu
        else
          if(u.lt.x) then
            a=u
          else
            b=u
          endif
          if(fu.le.fw .or. w.eq.x) then
            v=w
            fv=fw
            w=u
            fw=fu
          else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
            v=u
            fv=fu
          endif
        endif
11    continue
      pause 'brent exceed maximum iterations'
3     xmin=x
      brent=fx
      return
      END

      FUNCTION f1dim(x)
      INTEGER NM,NMAX
      REAL f1dim,VALUE,x
      PARAMETER (NM=70,NMAX=3*NM)
CU    USES VALUE
      INTEGER j,ncom
      REAL pcom(NMAX),xicom(NMAX),xt(NMAX)
      COMMON /f1com/ pcom,xicom,ncom
      EXTERNAL VALUE
      do 11 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
11    continue
      f1dim=VALUE(xt)
      return
      END

      SUBROUTINE linmin(p,xi,n,fret)
      PARAMETER (NM=70,NMAX=3*NM,TOL=1.e-5)
      INTEGER n,NM,NMAX
      REAL fret,p(NMAX),xi(NMAX),TOL
CU    USES brent,f1dim,mnbrak
      INTEGER j,ncom
      REAL ax,bx,fa,fb,fx,xmin,xx,pcom(NMAX),xicom(NMAX),brent
      COMMON /f1com/ pcom,xicom,ncom
      EXTERNAL f1dim
      ncom=n
      do 11 j=1,n
        pcom(j)=p(j)
        xicom(j)=xi(j)
11    continue
      ax=0.
      xx=1.
      call mnbrak(ax,xx,bx,fa,fx,fb,f1dim)
      fret=brent(ax,xx,bx,f1dim,TOL,xmin)
      do 12 j=1,n
        xi(j)=xmin*xi(j)
        p(j)=p(j)+xi(j)
12    continue
      return
      END

      SUBROUTINE frprmn(p,n,ftol,iter,fret)
      PARAMETER (NM=70,NMAX=3*NM,ITMAX=10000,EPS=1.e-10)
      INTEGER iter,n,NM,NMAX,ITMAX
      REAL fret,ftol,p(NMAX),EPS,VALUE
      EXTERNAL VALUE
CU    USES FORCE,VALUE,linmin
      INTEGER its,j
      REAL dgg,fp,gam,gg,g(NMAX),h(NMAX),xi(NMAX)
      fp=VALUE(p)
      call FORCE(p,xi)
      do 11 j=1,n
        g(j)=-xi(j)
        h(j)=g(j)
        xi(j)=h(j)
11    continue
      do 14 its=1,ITMAX
        iter=its
        call linmin(p,xi,n,fret)
        if(2.*abs(fret-fp).le.ftol*(abs(fret)+abs(fp)+EPS))return
        fp=VALUE(p)
        call FORCE(p,xi)
        gg=0.
        dgg=0.
        do 12 j=1,n
          gg=gg+g(j)**2
C         dgg=dgg+xi(j)**2
          dgg=dgg+(xi(j)+g(j))*xi(j)
12      continue
        if(gg.eq.0.)return
        gam=dgg/gg
        do 13 j=1,n
          g(j)=-xi(j)
          h(j)=g(j)+gam*h(j)
          xi(j)=h(j)
13      continue


14    continue
      pause 'frprmn maximum iterations exceeded'
      return
      END



      SUBROUTINE TERSOFF_FORCE
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=50,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),LAST,
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM),N,MODE,VIRILSUM(3,3),NEW,
     A                HSUM(3,3),STRESSUM(3,3),STRESS(3,3),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2,R_SI_SI,S_SI_SI,
     A                R_C_C,S_C_C,R_SI_C,S_SI_C
      LOGICAL WRITE_LP

      DIMENSION ZET2(NX,NX),ZET3(NX,NX),FC(NM),FCP(NM),
     A  RX(NX),RY(NX),RZ(NX),RXT(NX),RYT(NX),RZT(NX),
     A  R(NX), RJK2(NX,NX),RJK(NX,NX),RJKX(NX,NX),RJKY(NX,NX),
     A  RJKZ(NX,NX),PRJKX(NX,NX),PRJKY(NX,NX),PRJKZ(NX,NX),
     A  ZET(NX),SZET1(NX),SZET11(NX),SFTM2X(NX),SFTM2Y(NX),
     A  SFTM2Z(NX),SFTM3X(NX),SFTM3Y(NX),SFTM3Z(NX),PF(NX),
     A  FRST(NX),FF2(NX),FF3(NX),R2(NX)

      DO 10 I=1,N
      FX(I) = 0.E0
      FY(I) = 0.E0
      FZ(I) = 0.E0
      EP(I) = 0.E0
 10   CONTINUE

      DO 500 IPT = 1,N
      N1 = NTYPE(IPT)
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)
      NEG = 0

      DO 600 J = 1,INDEX(IPT)
      JPT = LIST(IPT,J)
      N2  = NTYPE(JPT)
      RSMAX = PS(N1)*PS(N2)
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      RR   = RXX*RXX+RYY*RYY+RZZ*RZZ
      IF(RR.GT.RSMAX) GOTO 600

      NEG = NEG + 1
      R1  = SQRT(RR)
      RX(NEG)  = RXX
      RY(NEG)  = RYY
      RZ(NEG)  = RZZ
      RXT(NEG) = RXX/R1
      RYT(NEG) = RYY/R1
      RZT(NEG) = RZZ/R1
      R2(NEG)  = RR
      R(NEG)   = R1
      LIST2(NEG) = JPT

      CA(NEG) = SQRT(PA(N1)*PA(N2))
      CB(NEG) = SQRT(PB(N1)*PB(N2))
      CR(NEG) = SQRT(PR(N1)*PR(N2))
      CS(NEG) = SQRT(PS(N1)*PS(N2))

      if (n1.eq.1.and.n2.eq.1) then 
      cr(neg) = r_si_si
      cs(neg) = s_si_si
      elseif (n1.eq.2.and.n2.eq.2) then 
      cr(neg) = r_c_c
      cs(neg) = s_c_c
      else
      cr(neg) = r_si_c
      cs(neg) = s_si_c
      endif

      RLM1(NEG) = (PL1(N1)+PL1(N2))*0.5
      RLM2(NEG) = (PL2(N1)+PL2(N2))*0.5
      IF(N1.EQ.N2) THEN
      CHIJ(NEG) = 1.E0
      ELSE
      CHIJ(NEG) = CHI_SI_C
      ENDIF

      IF(R1.LT.CR(NEG)) THEN
       FC(NEG)   = 1.E0
       FCP(NEG)  = 0.E0
      ELSE
       FC(NEG)=0.5*(1.E0+COS(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG))))
       FCP(NEG)=-0.5*PI*SIN(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG)))/(CS(NEG)-CR(NEG))
      ENDIF
 600  CONTINUE


      DO 700 J = 1,NEG
      JPT = LIST2(J)
      N22=NTYPE(JPT)
      ZET(J) = 0.E0
      SZET1(J)  = 0.E0
      SZET11(J) = 0.E0
      SFTM2X(J) = 0.E0
      SFTM2Y(J) = 0.E0
      SFTM2Z(J) = 0.E0
      SFTM3X(J) = 0.E0
      SFTM3Y(J) = 0.E0
      SFTM3Z(J) = 0.E0

      DO 800 K = 1,NEG
      IF(K.EQ.J) GOTO 800
      KP = LIST2(K)

      RMU = (RX(J)*RX(K)+RY(J)*RY(K)+RZ(J)*RZ(K))/R(J)/R(K)
      RJK2(J,K) = R2(J)+R2(K)-2.E0*RMU*R(J)*R(K)
      RJK(J,K)  = SQRT(RJK2(J,K))
      RJKX(J,K) = RX(J)-RX(K)
      RJKY(J,K) = RY(J)-RY(K)
      RJKZ(J,K) = RZ(J)-RZ(K)
      PRJKX(J,K) = RJKX(J,K)/RJK(J,K)
      PRJKY(J,K) = RJKY(J,K)/RJK(J,K)
      PRJKZ(J,K) = RJKZ(J,K)/RJK(J,K)

      RMU1 = RMU-PH(N1)
      RMU2 = RMU1*RMU1
      RMU3 = PD(N1)+RMU2
      GMU  = C2D2(N1)-PC(N1)/RMU3
      GMUP = 2.E0*RMU1*PC(N1)/RMU3**2

      COM1 = FCP(K)*GMU
      COM2 = FC(K)*GMUP
      ZET(J) = ZET(J)+FC(K)*GMU
      
      TIP1 = RMU/R(J)-1.E0/R(K)
      TIP2 = RMU/R(K)-1.E0/R(J)
      TIP3 = RJK(J,K)/R(J)/R(K)
      
      ZET1 = COM2*(-TIP1)
      SZET1(J)  = SZET1(J)+ZET1
      ZET2(J,K) = COM1+COM2*(-TIP2)
      ZET3(J,K) = COM2*(-TIP3)

      SFTM2X(J) = SFTM2X(J) + ZET2(J,K)*RXT(K)
      SFTM2Y(J) = SFTM2Y(J) + ZET2(J,K)*RYT(K)
      SFTM2Z(J) = SFTM2Z(J) + ZET2(J,K)*RZT(K)
      SFTM3X(J) = SFTM3X(J) + ZET3(J,K)*PRJKX(J,K)
      SFTM3Y(J) = SFTM3Y(J) + ZET3(J,K)*PRJKY(J,K)
      SFTM3Z(J) = SFTM3Z(J) + ZET3(J,K)*PRJKZ(J,K)

 800  CONTINUE
 700  CONTINUE


      DO 704 J=1,NEG
      JPT = LIST2(J)
      N22 = NTYPE(JPT)

      ELM2  = EXP(-RLM2(J)*R(J))
      ELM2B = ELM2*CB(J)
      ELM1  = EXP(-RLM1(J)*R(J))
      ELM1A = ELM1*CA(J)
      FCD1  = FCP(J)-RLM1(J)*FC(J)
      FCD2  = FCP(J)-RLM2(J)*FC(J)
 
      IF(ZET(J).NE.0.E0) THEN
       BEZET  = 1.E0+ PBT(N1)*(ZET(J)**PN(N1))
       BIJ    = (BEZET**PWR(N1))*CHIJ(J)
       BIJ1   = 0.5*PBT(N1)*(ZET(J)**PWR2(N1))
       BIJP   = BIJ1*(BEZET**PWR1(N1))*CHIJ(J)
       PF(J)  = FC(J)*ELM2B*BIJP
      ELSE
       BIJ = 1.E0
       PF(J)  = 0.E0
      ENDIF

      FRST(J)  = FCD1*ELM1A - BIJ*FCD2*ELM2B

C     CALCULATE POTENTIAL ENERGY PER PARTICLE
      VIJR = FC(J)*ELM1A
      VIJA = -FC(J)*BIJ*ELM2B
      VIJ  = VIJR+VIJA
      EP(IPT) = EP(IPT) + VIJ/4.E0
      EP(JPT) = EP(JPT) + VIJ/4.E0

C     THREE BODY FORCES 
      DO 840 K=1,NEG
      KPT = LIST2(K)
      IF (K.EQ.J) GOTO 840

      FF2(J) = PF(J)*ZET2(J,K)
      FF3(J) = PF(J)*ZET3(J,K)

C     CALC FORCE ON PARTICLE K
      FKX = -(FF2(J)*RX(K)/R(K)+FF3(J)*(-PRJKX(J,K)))*0.5
      FKY = -(FF2(J)*RY(K)/R(K)+FF3(J)*(-PRJKY(J,K)))*0.5
      FKZ = -(FF2(J)*RZ(K)/R(K)+FF3(J)*(-PRJKZ(J,K)))*0.5
      FX(KPT) = FX(KPT) + FKX
      FY(KPT) = FY(KPT) + FKY
      FZ(KPT) = FZ(KPT) + FKZ
 840  CONTINUE
C     --------------------------------------------------------------

C     CALC FORCE ON PARTICLE I
      FIX = ((FRST(J)+PF(J)*SZET1(J))*RXT(J)+PF(J)*SFTM2X(J))*0.5
      FIY = ((FRST(J)+PF(J)*SZET1(J))*RYT(J)+PF(J)*SFTM2Y(J))*0.5
      FIZ = ((FRST(J)+PF(J)*SZET1(J))*RZT(J)+PF(J)*SFTM2Z(J))*0.5
      FX(IPT) = FX(IPT)+FIX
      FY(IPT) = FY(IPT)+FIY
      FZ(IPT) = FZ(IPT)+FIZ

C     CALC FORCE ON PARTICLE J
      FJX=-((FRST(J)+PF(J)*SZET1(J))*RXT(J)+SFTM3X(J)*PF(J))*0.5
      FJY=-((FRST(J)+PF(J)*SZET1(J))*RYT(J)+SFTM3Y(J)*PF(J))*0.5
      FJZ=-((FRST(J)+PF(J)*SZET1(J))*RZT(J)+SFTM3Z(J)*PF(J))*0.5
      FX(JPT) = FX(JPT)+FJX
      FY(JPT) = FY(JPT)+FJY
      FZ(JPT) = FZ(JPT)+FJZ
     
 704  CONTINUE
 500  CONTINUE

      POTE = 0.E0
      DO 400 I =1,N
 400  POTE  = POTE  + EP(I)
 
      if (last.ne.1) print *,pote

      RETURN
      END


      SUBROUTINE TERSOFF_VALUE
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=50,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),LAST,
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM),N,MODE,VIRILSUM(3,3),NEW,
     A                HSUM(3,3),STRESSUM(3,3),STRESS(3,3),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2,R_SI_SI,S_SI_SI,
     A                R_C_C,S_C_C,R_SI_C,S_SI_C

      LOGICAL WRITE_LP,UPDATE,NEW

      DIMENSION FC(NM),RX(NX),RY(NX),RZ(NX),RXT(NX),RYT(NX),RZT(NX),
     A          R(NX), ZET(NX),R2(NX)

      DO 65 I=1,N
 65   EP(I) = 0.

      IF(.NOT.NEW) THEN
C     ------  CHECK FOR NLIST UPDATE ------
      DISPMX = 0.E0
      DO 30 I = 1,N
      DX(I) = X(I)-XOLD(I)
      DY(I) = Y(I)-YOLD(I)
      DZ(I) = Z(I)-ZOLD(I)
      DR = SQRT((DX(I)-DXOLD(I))**2+(DY(I)-DYOLD(I))**2+
     A          (DZ(I)-DZOLD(I))**2)
      IF (DR.GT.DISPMX) DISPMX = DR
30    CONTINUE
      UPDATE = (2.5E0*DISPMX).GT.(RLIST-RCUT)
      ENDIF

      IF (UPDATE) THEN
C    ---------------------------------------------------------  
C    UPDATE THE VERLET LIST

      IF (NEW) NEW=.FALSE.
      DO 160 IPT = 1,N
 160  INDEX(IPT) = 0
      
      DO 170 IPT = 1,N
      N1 = NTYPE(IPT)
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)

      DO 270 JPT = IPT+1,N
      N2=NTYPE(JPT)
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      R1  = SQRT(RXX*RXX+RYY*RYY+RZZ*RZZ)

      IF (R1.LT.RLIST) THEN
      INDEX(IPT)=INDEX(IPT)+1
      LIST(IPT,INDEX(IPT))=JPT
      INDEX(JPT)=INDEX(JPT)+1
      LIST(JPT,INDEX(JPT))=IPT
      ENDIF
 270  CONTINUE
 170  CONTINUE
      
      DO 450 I = 1, N
      DXOLD(I) = DX(I)
      DYOLD(I) = DY(I)
      DZOLD(I) = DZ(I)
 450  CONTINUE
C     ---------------------------------------------------------
      ENDIF

      DO 500 IPT = 1,N
      N1  = NTYPE(IPT)
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)
      NEG = 0

      DO 600 J = 1,INDEX(IPT)
      JPT = LIST(IPT,J)
      N2  = NTYPE(JPT)

      RSMAX = PS(N1)*PS(N2)
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      RR = RXX*RXX+RYY*RYY+RZZ*RZZ
      IF(RR.GT.RSMAX) GOTO 600

      NEG = NEG + 1
      R1  = SQRT(RR)
      RX(NEG)  = RXX
      RY(NEG)  = RYY
      RZ(NEG)  = RZZ
      RXT(NEG) = RXX/R1
      RYT(NEG) = RYY/R1
      RZT(NEG) = RZZ/R1
      R2(NEG)  = RR
      R(NEG)   = R1
      LIST2(NEG) = JPT

      CA(NEG) = SQRT(PA(N1)*PA(N2))
      CB(NEG) = SQRT(PB(N1)*PB(N2))
      CR(NEG) = SQRT(PR(N1)*PR(N2))
      CS(NEG) = SQRT(PS(N1)*PS(N2))
      RLM1(NEG) = (PL1(N1)+PL1(N2))*0.5
      RLM2(NEG) = (PL2(N1)+PL2(N2))*0.5

      if (n1.eq.1.and.n2.eq.1) then 
      cr(neg) = r_si_si
      cs(neg) = s_si_si
      elseif (n1.eq.2.and.n2.eq.2) then 
      cr(neg) = r_c_c
      cs(neg) = s_c_c
      else
      cr(neg) = r_si_c
      cs(neg) = s_si_c
      endif

      IF(N1.EQ.N2) THEN
      CHIJ(NEG) = 1.E0
      ELSE
      CHIJ(NEG) = CHI_SI_C
      ENDIF

      IF(R1.LT.CR(NEG)) THEN
       FC(NEG)   = 1.E0
      ELSE
       FC(NEG)=0.5*(1.E0+COS(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG))))
      ENDIF
 600  CONTINUE

      DO 700 J = 1,NEG
      JPT = LIST2(J)
      N22=NTYPE(JPT)
      ZET(J) = 0.E0
      DO 800 K = 1,NEG
      IF(K.EQ.J) GOTO 800
      KP = LIST2(K)
      RMU = (RX(J)*RX(K)+RY(J)*RY(K)+RZ(J)*RZ(K))/R(J)/R(K)
      RMU1 = RMU-PH(N1)
      RMU2 = RMU1*RMU1
      RMU3 = PD(N1)+RMU2
      GMU  = C2D2(N1)-PC(N1)/RMU3
      ZET(J) = ZET(J)+FC(K)*GMU
 800  CONTINUE
 700  CONTINUE


      DO 704 J=1,NEG
      JPT = LIST2(J)
      N22 = NTYPE(JPT)
      ELM2  = EXP(-RLM2(J)*R(J))
      ELM2B = ELM2*CB(J)
      ELM1  = EXP(-RLM1(J)*R(J))
      ELM1A = ELM1*CA(J)
      IF(ZET(J).NE.0.E0) THEN
       BEZET  = 1.E0+ PBT(N1)*(ZET(J)**PN(N1))
       BIJ    = (BEZET**PWR(N1))*CHIJ(J)
      ELSE
       BIJ = 1.E0
      ENDIF

C     CALCULATE POTENTIAL ENERGY PER PARTICLE
      VIJR = FC(J)*ELM1A
      VIJA = -FC(J)*BIJ*ELM2B
      VIJ  = VIJR+VIJA
      EP(IPT) = EP(IPT) + VIJ/4.E0
      EP(JPT) = EP(JPT) + VIJ/4.E0

 704  CONTINUE
 500  CONTINUE

      POTE = 0.
      DO 400 I =1,N
 400  POTE  = POTE+EP(I)

      RETURN
      END


      SUBROUTINE RANDOM
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (K=5701,J=3612,M=566927,RM=566927.E0)
      COMMON SEED

      IX = INT(SEED*RM)
      KR = J*IX + K
      IRAND = MOD(KR,M)
      SEED = (FLOAT(IRAND)+0.5)/RM
      RETURN
      END





