      SUBROUTINE INITIALIZE_TERSOFF_CONSTANTS
#     include "mulcp.fh"
C     SET OF PARAMTERS FOR Si. 
      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
      PD(1)   = (1.6217E1)**2
      PH(1)   = -5.9825E-1
C     PR(1)   = 2.70
C     PS(1)   = 3.00
      PR(1)   = 2.36
      PS(1)   = 2.56
      C2D2(1) = ONE+PC(1)/PD(1)
      PWR(1)  = -ONE/(2*PN(1))
      PWR1(1) = PWR(1)-ONE
      PWR2(1) = PN(1)-ONE

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
      PD(2)   = (7.0340)**2
      PH(2)   = -3.9953E-1
C     PR(2)   = 1.80
C     PS(2)   = 2.10
      PR(2)   = 2.36
      PS(2)   = 2.56
      C2D2(2) = ONE+PC(2)/PD(2)
      PWR(2)  = -ONE/(2*PN(2))
      PWR1(2) = PWR(2)-ONE
      PWR2(2) = PN(2)-ONE

      CHI_SI_C = 1.0086
      PM(1) = 28.086
      PM(2) = 12.011 
      RETURN
      END
      
      
      SUBROUTINE TERSOFF_DYNAMICAL
#     include "mulcp.fh"
      
      DIMENSION ZET2(NM,NX,NX),ZET3(NM,NX,NX),FC(NM),FCP(NM),
     A  ZET12(NX,NX),ZET13(NX,NX),ZET22(NX,NX),ZET23(NX,NX),
     A  ZET33(NX,NX),RX(NM,NX),RY(NM,NX),RZ(NM,NX),R(NM,NX),
     A  R2(NM,NX),RXT(NX),RYT(NX),RZT(NX),RJK2(NX,NX),
     A  RJK(NX,NX),RJKX(NX,NX),RJKY(NX,NX),FF11(NM,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  PF2(NX),PF3(NM,NX),FRST(NX),FRST2(NX),FF1(NM,NX),
     A  FF22(NM,NX,NX),FF23(NM,NX,NX),FF13(NM,NX,NX),FF12(NM,NX,NX),
     A  FF3(NM,NX,NX),FF2(NM,NX,NX),FF33(NM,NX,NX),FCPP(NM),
     A  LIST3(1+NX+NX*NX)
#ifdef matrix_generation_precision_is_real   
      REAL
#else
      DOUBLE PRECISION
#endif      
     a     lloopxx,lloopxy,nloopxx,nloopxy,
     a     lloop2xx,lloop2xy,nloop2xx,nloop2xy,
     a     lloopyx,nloopyx,lloop2yx,
     a     nloop2yx,lloopyy,nloopyy,lloop2yy,
     a     nloop2yy,lloopzz,lloopxz,nloopzz,
     a     nloopxz,lloop2zz,lloop2xz,nloop2zz,
     a     nloop2xz,lloopyz, lloopzx, lloopzy,
     a     nloopyz, nloopzx, nloopzy, lloop2yz,
     a     lloop2zx, lloop2zy, nloop2yz,nloop2zx,
     a     nloop2zy,
     a     nloop2,lloop3xx,lloop3xy,lloop3xz,
     a     lloop3yx,lloop3yy,lloop3yz,lloop3zx,
     a     lloop3zy,lloop3zz,lloop4xx,lloop4xy,
     a     lloop4xz,lloop4yx,lloop4yy,lloop4yz,
     a     lloop4zx,lloop4zy,lloop4zz,
     a     nloop3xx,nloop3xy,nloop3xz,nloop3yx,
     a     nloop3yy,nloop3yz,nloop3zx,nloop3zy,
     a     nloop3zz,nloop4xx,nloop4xy,nloop4xz,
     a     nloop4yx,nloop4yy,nloop4yz,nloop4zx,
     a     nloop4zy,nloop4zz,
     a     lloop5xx,lloop5yy,lloop5zz,lloop5xy,
     a     lloop5xz,lloop5yx,lloop5yz,lloop5zx,
     a     lloop5zy

      DO 10 IPT=1,N
      FX(IPT) = ZERO
      FY(IPT) = ZERO 
      FZ(IPT) = ZERO
      EP(IPT) = ZERO
 10   CONTINUE

      DO 500 IPT = 1,N
      N1 = NTYPE(IPT)
      SXI = SX(IPT)
      SYI = SY(IPT)
      SZI = SZ(IPT)
      NEG = 0

      DO 600 J = 1,INDEX(IPT)
      JPT = LIST(IPT,J)
      N2  = NTYPE(JPT)
      RSMAX = PS(N1)*PS(N2)
      SXJI  = SX(JPT) - SXI
      SYJI  = SY(JPT) - SYI
      SZJI  = SZ(JPT) - SZI
      IF(SXJI.GT.HALF)  SXJI=SXJI-ONE
      IF(SXJI.LT.-HALF) SXJI=SXJI+ONE
      IF(SYJI.GT.HALF)  SYJI=SYJI-ONE
      IF(SYJI.LT.-HALF) SYJI=SYJI+ONE
      IF(SZJI.GT.HALF)  SZJI=SZJI-ONE
      IF(SZJI.LT.-HALF) SZJI=SZJI+ONE
      RXX  = H(1,1)*SXJI + H(2,1)*SYJI + H(3,1)*SZJI
      RYY  = H(1,2)*SXJI + H(2,2)*SYJI + H(3,2)*SZJI
      RZZ  = H(1,3)*SXJI + H(2,3)*SYJI + H(3,3)*SZJI
      RR   = RXX*RXX+RYY*RYY+RZZ*RZZ      
      IF(RR.GT.RSMAX) GOTO 600

      NEG = NEG + 1
      IF (NEG.GT.NX) THEN 
      PRINT*,'Neighbourlist of atom',ipt,' more than mulcp can handle.'
      PRINT*,'Increase _NX in mulcp.fh and recompile.'
      STOP
      ENDIF
      
      LIST2(IPT,NEG) = JPT
      R1  = SQRT(RR)
      RX(IPT,NEG)  = RXX
      RY(IPT,NEG)  = RYY
      RZ(IPT,NEG)  = RZZ
      R2(IPT,NEG)  = RR
      R(IPT,NEG)   = R1

      RXT(NEG) = RXX/R1
      RYT(NEG) = RYY/R1
      RZT(NEG) = RZZ/R1

      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))

c     if (n1.eq.1.and.n2.eq.1) then 
c     cr(neg) = 1.6
c     cs(neg) = 2.9
c     elseif (n1.eq.2.and.n2.eq.2) then 
c     cr(neg) = 1.8
c     cs(neg) = 2.1
c     else
c     cr(neg) = 2.2
c     cs(neg) = 2.5
c     endif

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

      IF(R1.LT.CR(NEG)) THEN
       FC(NEG)   = ONE
       FCP(NEG)  = ZERO
       FCPP(NEG) = ZERO
      ELSE
       Fc(NEG)=HALF*(ONE+COS(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG))))
       FCP(NEG)=-HALF*PI*
     A     SIN(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG)))/(CS(NEG)-CR(NEG))
       FCPP(NEG)=(HALF-FC(NEG))*PI*
     A     PI/(CS(NEG)-CR(NEG))
     A     /(CS(NEG)-CR(NEG))
      ENDIF
 600  CONTINUE
      INDEX2(IPT) = NEG

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

      ZET(J) = ZERO
      SZET1(J)  = ZERO
      SZET11(J) = ZERO
      SFTM2X(J) = ZERO
      SFTM2Y(J) = ZERO
      SFTM2Z(J) = ZERO
      SFTM3X(J) = ZERO
      SFTM3Y(J) = ZERO
      SFTM3Z(J) = ZERO

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

      RMU = (RX(IPT,J)*RX(IPT,K)+RY(IPT,J)*RY(IPT,K)+
     a       RZ(IPT,J)*RZ(IPT,K))/R(IPT,J)/R(IPT,K)
      RJK2(J,K) = R2(IPT,J)+R2(IPT,K)
     a            -2*RMU*R(IPT,J)*R(IPT,K)
      RJK(J,K)  = SQRT(RJK2(J,K))
      RJKX(J,K) = RX(IPT,J)-RX(IPT,K)
      RJKY(J,K) = RY(IPT,J)-RY(IPT,K)
      RJKZ(J,K) = RZ(IPT,J)-RZ(IPT,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*RMU1*PC(N1)/RMU3**2

      COM1 = FCP(K)*GMU
      COM2 = FC(K)*GMUP
      COM3 = FCPP(K)*GMU
      ZET(J) = ZET(J)+FC(K)*GMU
      
      TIP1 = RMU/R(IPT,J)-ONE/R(IPT,K)
      TIP2 = RMU/R(IPT,K)-ONE/R(IPT,J)
      TIP3 = RJK(J,K)/R(IPT,J)/R(IPT,K)
      
      ZET1 = COM2*(-TIP1)
      SZET1(J)  = SZET1(J)+ZET1
      ZET2(IPT,J,K) = COM1+COM2*(-TIP2)
      ZET3(IPT,J,K) = COM2*(-TIP3)

      GMUPP = 2*PC(N1)/RMU3**2-8.E0*PC(N1)*RMU2/RMU3**3
      COM4  = FCP(K)*GMUP
      COM5  = FC(K)*GMUPP

      TIP11 = 2*RMU/R2(IPT,J)-ONE/R(IPT,J)/R(IPT,K)
      TIP12 = RMU/R(IPT,J)/R(IPT,K)-ONE/R2(IPT,K)-ONE/R2(IPT,J)
      TIP13 = RJK(J,K)/R2(IPT,J)/R(IPT,K)
      TIP22 = 2*RMU/R2(IPT,K)-ONE/R(IPT,J)/R(IPT,K)
      TIP23 = RJK(J,K)/R2(IPT,K)/R(IPT,J)
      TIP33 = -ONE/R(IPT,J)/R(IPT,K)

      ZET11 = COM5*TIP1*TIP1+COM2*TIP11
      SZET11(J) = SZET11(J) + ZET11
      ZET12(J,K) = COM4*(-TIP1)+COM5*TIP2*TIP1+COM2*TIP12
      ZET13(J,K) = COM5*TIP3*TIP1+COM2*TIP13
      ZET22(J,K) = COM3+2*COM4*(-TIP2)+COM5*TIP2*TIP2+COM2*TIP22
      ZET23(J,K) = COM4*(-TIP3)+COM5*TIP2*TIP3+COM2*TIP23
      ZET33(J,K) = COM5*TIP3*TIP3+COM2*TIP33

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

 800  CONTINUE
 700  CONTINUE

      DO 704 J=1,NEG

      JPT = LIST2(IPT,J)
      N22 = NTYPE(JPT)

      ELM2  = EXP(-RLM2(J)*R(IPT,J))
      ELM2B = ELM2*CB(J)
      ELM1  = EXP(-RLM1(J)*R(IPT,J))
      ELM1A = ELM1*CA(J)
      FCD1  = FCP(J)-RLM1(J)*FC(J)
      FCD2  = FCP(J)-RLM2(J)*FC(J)
      FCD1P = FCPP(J)-2*RLM1(J)*FCP(J)+RLM1(J)**2.*FC(J)
      FCD2P = FCPP(J)-2*RLM2(J)*FCP(J)+RLM2(J)**2.*FC(J)
 
      IF(ZET(J).NE.ZERO) THEN
       BEZET  = ONE+ PBT(N1)*(ZET(J)**PN(N1))
       BIJ    = (BEZET**PWR(N1))*CHIJ(J)
       BIJ1   = HALF*PBT(N1)*(ZET(J)**PWR2(N1))
       BIJP   = BIJ1*(BEZET**PWR1(N1))*CHIJ(J)
       PF(J)  = FC(J)*ELM2B*BIJP
       PF2(J) = FCD2*ELM2B*BIJP
       PF3(IPT,J) = PF(J)*(PWR2(N1)/ZET(J)-
     A      (1.+2.*PN(N1))*BIJ1/BEZET)
      ELSE
       BIJ = ONE
       PF(J)  = ZERO
       PF2(J) = ZERO
       PF3(IPT,J) = ZERO
      ENDIF

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

      FF1(IPT,J)  = FRST(J)+PF(J)*SZET1(J)

      FF11(IPT,J) = FRST2(J)+2*PF2(J)*SZET1(J)
     A     +PF3(IPT,J)*SZET1(J)**2+PF(J)*SZET11(J)

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 AND THREE BODY STRESSES
      DO 840 K=1,NEG
      
      KPT = LIST2(IPT,K)
      IF (K.EQ.J) GOTO 840

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

      FF12(IPT,J,K) = PF2(J)*ZET2(IPT,J,K)+PF3(IPT,J)*SZET1(J)
     A                *ZET2(IPT,J,K)+PF(J)*ZET12(J,K)
      FF13(IPT,J,K) = PF2(J)*ZET3(IPT,J,K)+PF(J)*ZET13(J,K)
     A                +PF3(IPT,J)*SZET1(J)*ZET3(IPT,J,K)
      FF22(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)**2+PF(J)*ZET22(J,K)
      FF23(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)*ZET3(IPT,J,K)
     A                +PF(J)*ZET23(J,K)
      FF33(IPT,J,K) = PF3(IPT,J)*ZET3(IPT,J,K)**2+PF(J)*ZET33(J,K)

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

C     CALC FORCE ON PARTICLE I
      FIX = ((FRST(J)+PF(J)*SZET1(J))*RXT(J)
     A      +PF(J)*SFTM2X(J))*HALF
      FIY = ((FRST(J)+PF(J)*SZET1(J))*RYT(J)
     A      +PF(J)*SFTM2Y(J))*HALF
      FIZ = ((FRST(J)+PF(J)*SZET1(J))*RZT(J)
     A      +PF(J)*SFTM2Z(J))*HALF
      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)
     A     +SFTM3X(J)*PF(J))*HALF
      FJY=-((FRST(J)+PF(J)*SZET1(J))*RYT(J)
     A     +SFTM3Y(J)*PF(J))*HALF
      FJZ=-((FRST(J)+PF(J)*SZET1(J))*RZT(J)
     A     +SFTM3Z(J)*PF(J))*HALF
      FX(JPT) = FX(JPT)+FJX
      FY(JPT) = FY(JPT)+FJY
      FZ(JPT) = FZ(JPT)+FJZ

 704  CONTINUE

 500  CONTINUE

      POTE = ZERO
      DO 400 I =1,N
      POTE  = POTE  + EP(I)
 400  CONTINUE
c     print *, 'pote =', pote/n,n
c      print *,' '

C     CALCULATE THE DYNAMICAL MATRIX
C     THE NOTATION FOLLOWS THE ANALYTIC DERIVATION BY L.J.PORTER
C     I.E, LLOOP AND LLOOP2 REFER TO THE FIRST SUM OVER L
C     AND THE SECOND SUM OVER L FOUND IN THOSE NOTES.
C
C     THE DERIVATIVE NOTATION IS A "HYBRID EINSTEIN NOTATION".
C     FOR EXAMPLE, DIJDIJVIL = THE SECOND DERIVATIVE OF VIL WITH
C     RESPECT TO RIJ. DIJDIKVIJ = THE SECOND DERIVATIVE OF VIJ WITH
C     RESPECT TO RIJ AND RIK. ETC

C     zero the dynamical matrix
      do 17 i = 1,3*n
      idx(i) = 0
 17   continue

C     MEIJIE DEFINED XIJ,YIJ,ZIJ AS XJ - XI, ETC
C     NOTE ALSO THAT BELOW, I HAVE QUANTITIES LIKE
C     RJNX = RX(I,N) - RX(I,J) BECAUSE AFTER THE
C     LOOPS BELOW, RX(I,J) WILL  = XI - XJ.

      do 21 ipt= 1,n
      do 22 j = 1,INDEX2(ipt)
      RX(ipt,j) = -RX(ipt,j)
      RY(ipt,j) = -RY(ipt,j)
      Rz(ipt,j) = -Rz(ipt,j)
 22   continue
 21   continue

      do 1002 ipt = 1,n

c     find the possible neighbours that could have non-zero 
c     dynamical matrix element.
      ipw = 0
      do 257 k=1,index2(ipt)
      ipw = ipw+1
      list3(ipw)=list2(ipt,k)
      jpt = list2(ipt,k)
      do 257 kk=1,index2(jpt)
      ipw = ipw+1
      list3(ipw)=list2(jpt,kk)
 257  continue
      do 45 i=1,ipw
      if (list3(i).eq.0) goto 45
      do 48 j=i+1,ipw
      if (list3(j).eq.list3(i)) list3(j) = 0
 48   continue
 45   continue
      
      do 1100 jkk = 1,ipw
      jpt = list3(jkk)
      if (jpt.eq.0.or.jpt.eq.ipt) goto 1100

      firstxx = ZERO
      firstyy = ZERO
      firstzz = ZERO
      firstxy = ZERO
      firstxz = ZERO
      firstyx = ZERO
      firstyz = ZERO
      firstzx = ZERO
      firstzy = ZERO

      lloopxx = ZERO
      lloopyy = ZERO
      lloopzz = ZERO
      lloopxy = ZERO
      lloopxz = ZERO
      lloopyx = ZERO
      lloopyz = ZERO
      lloopzx = ZERO
      lloopzy = ZERO
     
      lloop2xx = ZERO
      lloop2yy = ZERO
      lloop2zz = ZERO
      lloop2xy = ZERO
      lloop2xz = ZERO
      lloop2yx = ZERO
      lloop2yz = ZERO
      lloop2zx = ZERO
      lloop2zy = ZERO

      lloop5xx = ZERO
      lloop5yy = ZERO
      lloop5zz = ZERO
      lloop5xy = ZERO
      lloop5xz = ZERO
      lloop5yx = ZERO
      lloop5yz = ZERO
      lloop5zx = ZERO
      lloop5zy = ZERO

      negij = 0
      do 654 iw = 1,index2(ipt)
 654  if (list2(ipt,iw).eq.jpt) negij = iw

      negji = 0
      do 659 iw = 1,index2(jpt)
 659  if (list2(jpt,iw).eq.ipt) negji = iw

      if (negij.eq.0) goto 1800

      rij  = r(ipt,negij)
      rijx = RX(ipt,negij)
      rijy = RY(ipt,negij)
      rijz = RZ(ipt,negij)
      rij2 = r2(ipt,negij)
      rij3 = rij2*rij

      dijvij = ff1(ipt,negij)
      dijdijvij = ff11(ipt,negij)
      dijvji = ff1(jpt,negji)
      dijdijvji = ff11(jpt,negji)

      firstxx = (rijx*rijx/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijx*(-rijx)/rij2)*dijdijvij + (rijx*(-rijx)/rij2)*dijdijvji
c
      firstyy = (rijy*rijy/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijy*(-rijy)/rij2)*dijdijvij + (rijy*(-rijy)/rij2)*dijdijvji
c
      firstxy = (rijx*rijy/rij3)*(dijvij + dijvji) +
     a  (rijx*(-rijy)/rij2)*dijdijvij + (rijx*(-rijy)/rij2)*dijdijvji
c
      firstyx = (rijy*rijx/rij3)*(dijvij + dijvji) +
     a  (rijy*(-rijx)/rij2)*dijdijvij + (rijy*(-rijx)/rij2)*dijdijvji
c
      firstzz = (rijz*rijz/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijz*(-rijz)/rij2)*dijdijvij + (rijz*(-rijz)/rij2)*dijdijvji
c
      firstxz = (rijx*rijz/rij3)*(dijvij + dijvji) +
     a  (rijx*(-rijz)/rij2)*dijdijvij + (rijx*(-rijz)/rij2)*dijdijvji
c
      firstyz = (rijy*rijz/rij3)*(dijvij + dijvji) +
     a  (rijy*(-rijz)/rij2)*dijdijvij + (rijy*(-rijz)/rij2)*dijdijvji
c
      firstzx = (rijz*rijx/rij3)*(dijvij + dijvji) +
     a  (rijz*(-rijx)/rij2)*dijdijvij + (rijz*(-rijx)/rij2)*dijdijvji
c
      firstzy = (rijz*rijy/rij3)*(dijvij + dijvji) +
     a  (rijz*(-rijy)/rij2)*dijdijvij + (rijz*(-rijy)/rij2)*dijdijvji

      do 1150 l = 1, INDEX2(ipt)
      if (l .eq. negij) goto 1150
      ril = r(ipt,l)
      ril2 = ril*ril
      rilx = RX(ipt,l)
      rily = RY(ipt,l)
      rilz = RZ(ipt,l)
      rmujl = (RX(ipt,negij)*RX(ipt,l) + RY(ipt,negij)*RY(ipt,l) + 
     a     RZ(ipt,negij)*RZ(ipt,l))/(rij*ril)
      rjl2 = rij2 + ril2 -2.0*rmujl*rij*ril
      rjlx = RX(ipt,l)-RX(ipt,negij)
      rjly = RY(ipt,l)-RY(ipt,negij)
      rjlz = RZ(ipt,l)-RZ(ipt,negij)
      rjl = sqrt(rjl2)

      nloopxx = ZERO
      nloopxy = ZERO
      nloopyx = ZERO
      nloopyy = ZERO
      nloopzz = ZERO
      nloopxz = ZERO
      nloopyz = ZERO
      nloopzx = ZERO
      nloopzy = ZERO
     
      do 1200 m = 1,INDEX2(ipt)
      if ((m .eq. negij).or.(m.eq.l)) goto 1200
      rinx = Rx(ipt,m)
      riny = RY(ipt,m)
      rinz = RZ(ipt,m)
      rin  = r(ipt,m)

      dindjlvij = pf3(ipt,negij)*zet2(ipt,negij,m)*zet3(ipt,negij,l)
      dildijvin = pf3(ipt,m)*zet2(ipt,m,negij)*zet2(ipt,m,l)
      dindjlvil = pf3(ipt,l)*zet2(ipt,l,m)*zet3(ipt,l,negij)

      nloopxx = nloopxx + (rilx*(-rijx)/(rij*ril))*dildijvin +
     a     (rinx*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyy = nloopyy + (rily*(-rijy)/(rij*ril))*dildijvin +
     a     (riny*rjly/(rin*rjl))*(dindjlvij + dindjlvil)
c     
      nloopxy = nloopxy + (rilx*(-rijy)/(rij*ril))*dildijvin +
     a     (rinx*rjly/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyx = nloopyx + (rily*(-rijx)/(rij*ril))*dildijvin +
     a     (riny*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopzz = nloopzz + (rilz*(-rijz)/(rij*ril))*dildijvin +
     a     (rinz*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopxz = nloopxz + (rilx*(-rijz)/(rij*ril))*dildijvin +
     a     (rinx*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyz = nloopyz + (rily*(-rijz)/(rij*ril))*dildijvin +
     a     (riny*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopzx = nloopzx + (rilz*(-rijx)/(rij*ril))*dildijvin +
     a     (rinz*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c     
      nloopzy = nloopzy + (rilz*(-rijy)/(rij*ril))*dildijvin +
     a     (rinz*rjly/(rin*rjl))*(dindjlvij + dindjlvil)

1200  continue

          dijvil = ff2(ipt,l,negij)
          dijdijvil = ff22(ipt,l,negij)
          dijdjlvij = ff13(ipt,negij,l)
          dijdjlvil = ff23(ipt,l,negij)
          dijdilvil = ff12(ipt,l,negij)
          dijdilvij = ff12(ipt,negij,l)
          djldilvij = ff23(ipt,negij,l)
          djldilvil = ff13(ipt,l,negij)

          lloopxx = lloopxx + (rijx*rijx/rij3 - 1.0/rij)*dijvil +
     a           (rijx*(-rijx)/rij2)*dijdijvil +
     a           (rijx*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxx
c
          lloopyy = lloopyy + (rijy*rijy/rij3 - 1.0/rij)*dijvil +
     a           (rijy*(-rijy)/rij2)*dijdijvil +
     a           (rijy*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyy
c
          lloopxy = lloopxy + (rijx*rijy/rij3)*dijvil +
     a           (rijx*(-rijy)/rij2)*dijdijvil +
     a           (rijx*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxy
c
          lloopyx = lloopyx + (rijy*rijx/rij3)*dijvil +
     a           (rijy*(-rijx)/rij2)*dijdijvil +
     a           (rijy*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyx
c
          lloopzz = lloopzz + (rijz*rijz/rij3 - 1.0/rij)*dijvil +
     a           (rijz*(-rijz)/rij2)*dijdijvil +
     a           (rijz*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzz
c
          lloopxz = lloopxz + (rijx*rijz/rij3)*dijvil +
     a           (rijx*(-rijz)/rij2)*dijdijvil +
     a           (rijx*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxz
c
          lloopyz = lloopyz + (rijy*rijz/rij3)*dijvil +
     a           (rijy*(-rijz)/rij2)*dijdijvil +
     a           (rijy*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyz
c
          lloopzx = lloopzx + (rijz*rijx/rij3)*dijvil +
     a           (rijz*(-rijx)/rij2)*dijdijvil +
     a           (rijz*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzx
c
          lloopzy = lloopzy + (rijz*rijy/rij3)*dijvil +
     a           (rijz*(-rijy)/rij2)*dijdijvil +
     a           (rijz*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzy

 1150   continue

        do 1250 l = 1,INDEX2(jpt)
        lpt = LIST2(jpt,l)
        if(lpt .eq. ipt) goto 1250
        rjl = r(jpt,l)
        rjl2 = rjl*rjl
        rjlx = RX(jpt,l)
        rjly = RY(jpt,l)
        rjlz = RZ(jpt,l)
        rmuil = (RX(jpt,negji)*RX(jpt,l) + RY(jpt,negji)*RY(jpt,l) +
     a       RZ(jpt,negji)*RZ(jpt,l))/(rij*rjl)      
        ril2 = rij2 + rjl2 - 2.0*rmuil*rij*rjl
        ril = sqrt(ril2)
        rilx = RX(jpt,l)-RX(jpt,negji)
        rily = RY(jpt,l)-RY(jpt,negji)
        rilz = RZ(jpt,l)-RZ(jpt,negji)

        nloop2xx = ZERO
        nloop2xy = ZERO
        nloop2yx = ZERO
        nloop2yy = ZERO
        nloop2zz = ZERO
        nloop2xz = ZERO
        nloop2yz = ZERO
        nloop2zx = ZERO
        nloop2zy = ZERO
     
        do 1300 m = 1, INDEX2(jpt)
        npt = LIST2(jpt,m)
        if ((npt.eq.ipt) .or.(m.eq.l)) goto 1300
        rjn = r(jpt,m)
        rjn2 = rjn*rjn
        rmuin = (RX(jpt,negji)*RX(jpt,m) + RY(jpt,negji)*RY(jpt,m) +
     a       RZ(jpt,negji)*RZ(jpt,m))/(rij*rjn) 
        rin2 = rij2 + rjn2 - 2.0*rmuin*rij*rjn
        rin = sqrt(rin2)
        rinx = RX(jpt,m)-RX(jpt,negji)
        riny = RY(jpt,m)-RY(jpt,negji)
        rinz = RZ(jpt,m)-RZ(jpt,negji)

        djldijvjn = pf3(jpt,m)*zet2(jpt,m,negji)*zet2(jpt,m,l)
        djldinvji = pf3(jpt,negji)*zet2(jpt,negji,l)*zet3(jpt,negji,m)
        djldinvjn = pf3(jpt,m)*zet2(jpt,m,l)*zet3(jpt,m,negji)

        nloop2xx = nloop2xx + (rijx*rjlx/(rij*rjl))*djldijvjn 
     a       +(rinx*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2yy = nloop2yy + (rijy*rjly/(rij*rjl))*djldijvjn 
     a       +(riny*rjly/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2xy = nloop2xy + (rijx*rjly/(rij*rjl))*djldijvjn 
     a       +(rinx*rjly/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2yx = nloop2yx + (rijy*rjlx/(rij*rjl))*djldijvjn 
     a       +(riny*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2zz = nloop2zz + (rijz*rjlz/(rij*rjl))*djldijvjn 
     a       +(rinz*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2xz = nloop2xz + (rijx*rjlz/(rij*rjl))*djldijvjn 
     a       +(rinx*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2yz = nloop2yz + (rijy*rjlz/(rij*rjl))*djldijvjn 
     a       +(riny*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2zx = nloop2zx + (rijz*rjlx/(rij*rjl))*djldijvjn 
     a       +(rinz*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2zy = nloop2zy + (rijz*rjly/(rij*rjl))*djldijvjn 
     a       +(rinz*rjly/(rin*rjl))*(djldinvji + djldinvjn)

 1300   continue

        dijvjl = ff2(jpt,l,negji)
        djldijvji = ff12(jpt,negji,l)
        dijdijvjl = ff22(jpt,l,negji)
        dildjlvjl = ff13(jpt,l,negji)
        dijdjlvjl = ff12(jpt,l,negji)
        dildjlvji = ff23(jpt,negji,l)
        dijdilvji = ff13(jpt,negji,l)
        dijdilvjl = ff23(jpt,l,negji)
c
        lloop2xx = lloop2xx +(rijx*rijx/rij3 - 1.0/rij)*dijvjl +
     a       (rijx*rjlx/(rij*rjl))*djldijvji +
     a       (rijx*(-rijx)/rij2)*dijdijvjl +
     a       (rilx*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xx
c
        lloop2yy = lloop2yy +(rijy*rijy/rij3 - 1.0/rij)*dijvjl +
     a       (rijy*rjly/(rij*rjl))*djldijvji +
     a       (rijy*(-rijy)/rij2)*dijdijvjl +
     a       (rily*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjly/(rij*rjl))*dijdjlvjl +
     a       (rily*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yy
c
        lloop2xy = lloop2xy +(rijx*rijy/rij3)*dijvjl +
     a       (rijx*rjly/(rij*rjl))*djldijvji +
     a       (rijx*(-rijy)/rij2)*dijdijvjl +
     a       (rilx*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjly/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xy
c     
        lloop2yx = lloop2yx +(rijy*rijx/rij3)*dijvjl +
     a       (rijy*rjlx/(rij*rjl))*djldijvji +
     a       (rijy*(-rijx)/rij2)*dijdijvjl +
     a       (rily*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rily*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yx
c
        lloop2zz = lloop2zz +(rijz*rijz/rij3 - 1.0/rij)*dijvjl +
     a       (rijz*rjlz/(rij*rjl))*djldijvji +
     a       (rijz*(-rijz)/rij2)*dijdijvjl +
     a       (rilz*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zz
c
        lloop2xz = lloop2xz +(rijx*rijz/rij3)*dijvjl +
     a       (rijx*rjlz/(rij*rjl))*djldijvji +
     a       (rijx*(-rijz)/rij2)*dijdijvjl +
     a       (rilx*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xz
c
        lloop2yz = lloop2yz +(rijy*rijz/rij3)*dijvjl +
     a       (rijy*rjlz/(rij*rjl))*djldijvji +
     a       (rijy*(-rijz)/rij2)*dijdijvjl +
     a       (rily*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rily*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yz
c
        lloop2zx = lloop2zx +(rijz*rijx/rij3)*dijvjl +
     a       (rijz*rjlx/(rij*rjl))*djldijvji +
     a       (rijz*(-rijx)/rij2)*dijdijvjl +
     a       (rilz*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zx
c
        lloop2zy = lloop2zy +(rijz*rijy/rij3)*dijvjl +
     a       (rijz*rjly/(rij*rjl))*djldijvji +
     a       (rijz*(-rijy)/rij2)*dijdijvjl +
     a       (rilz*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjly/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zy

 1250   continue
 1800  continue

       do 1700 l = 1,INDEX2(ipt)
       lpt = LIST2(ipt,l)
       if (lpt .eq. jpt) goto 1700

       negjl = 0
       do 644 iw = 1,index2(jpt)
 644   if (list2(jpt,iw).eq.lpt) negjl = iw
       if (negjl .eq. 0) goto 1700

       negli = 0
       do 634 iw = 1,index2(lpt)
 634   if (list2(lpt,iw).eq.ipt) negli = iw
       neglj = 0
       do 624 iw = 1,index2(lpt)
 624   if (list2(lpt,iw).eq.jpt) neglj = iw

       ril = r(ipt,l)
       rilx = RX(ipt,l)
       rily = RY(ipt,l)
       rilz = RZ(ipt,l)
       ril2 = ril*ril
       rlj = r(lpt,neglj)
       rlj2 = rlj*rlj
       rljx = RX(lpt,neglj)
       rljy = RY(lpt,neglj)
       rljz = RZ(lpt,neglj)

c     The following is a calculation of the distance rij if j is not
c     a nearest neighbor of i. From above, l must be a neighbor of both
c     i and j, so we can get rij through l.

       if (negij .eq. 0) then
       rmuij=(RX(lpt,negli)*RX(lpt,neglj)+RY(lpt,negli)*RY(lpt,neglj)
     a        +RZ(lpt,negli)*RZ(lpt,neglj))/(ril*rlj)
        rij2 = ril2+rlj2-2.0*rmuij*ril*rlj
        rij = sqrt(rij2)
        rij3 = rij2*rij
        rijx = RX(lpt,neglj)-RX(lpt,negli) 
        rijy = RY(lpt,neglj)-RY(lpt,negli)
        rijz = RZ(lpt,neglj)-RZ(lpt,negli)
       endif 

       nloop3xx = ZERO
       nloop3xy = ZERO
       nloop3xz = ZERO
       nloop3yx = ZERO
       nloop3yy = ZERO
       nloop3yz = ZERO
       nloop3zx = ZERO
       nloop3zy = ZERO
       nloop3zz = ZERO
     
       do 1750 m = 1,INDEX2(lpt)
       npt = LIST2(lpt,m)
       if((npt.eq.ipt).or.(npt.eq.jpt)) goto 1750
       rln = r(lpt,m)
       rln2 = rln*rln
       rmuin = (RX(lpt,negli)*RX(lpt,m) + RY(lpt,negli)*RY(lpt,m) +
     a      RZ(lpt,negli)*RZ(lpt,m))/(ril*rln)
       rin2 = ril2 + rln2 - 2.0*rmuin*ril*rln
       rin = sqrt(rin2)
       rinx =  RX(lpt,m)-RX(lpt,negli)
       riny =  RY(lpt,m)-RY(lpt,negli)
       rinz =  RZ(lpt,m)-RZ(lpt,negli)
       rmujn = (RX(lpt,neglj)*RX(lpt,m) + RY(lpt,neglj)*RY(lpt,m) +
     a      RZ(lpt,neglj)*RZ(lpt,m))/(rlj*rln)           
       rjn2 = rln2 + rlj2 - 2.0*rmujn*rln*rlj
       rjn = sqrt(rjn2)
       rjnx = RX(lpt,m)-RX(lpt,neglj)
       rjny = RY(lpt,m)-RY(lpt,neglj)
       rjnz = RZ(lpt,m)-RZ(lpt,neglj)

       dindjlvli = pf3(lpt,negli)*zet2(lpt,negli,neglj)
     a      *zet3(lpt,negli,m)
       dindjlvln = pf3(lpt,m)*zet2(lpt,m,neglj)*zet3(lpt,m,negli)
       dildjnvlj= pf3(lpt,neglj)*zet2(lpt,neglj,negli)
     a      *zet3(lpt,neglj,m)
       dildjnvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet3(lpt,m,neglj)
       dindijvli = pf3(lpt,negli)*zet3(lpt,negli,m)
     a      *zet3(lpt,negli,neglj)
       dildjlvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet2(lpt,m,neglj)
       dijdjnvlj = pf3(lpt,neglj)*zet3(lpt,neglj,negli)
     a      *zet3(lpt,neglj,m)
       dindjnvln = pf3(lpt,m)*zet3(lpt,m,negli)*zet3(lpt,m,neglj)
     
       nloop3xx = (rinx*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a      (rilx*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a      (rinx*(-rijx)/(rin*rij))*dindijvli +
     a      (rilx*(-rljx)/(ril*rlj))*dildjlvln +
     a      (rijx*rjnx/(rij*rjn))*dijdjnvlj +
     a      (rinx*rjnx/(rin*rjn))*dindjnvln + nloop3xx
c
        nloop3yy = (riny*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a      (rily*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a      (riny*(-rijy)/(rin*rij))*dindijvli +
     a      (rily*(-rljy)/(ril*rlj))*dildjlvln +
     a      (rijy*rjny/(rij*rjn))*dijdjnvlj +
     a      (riny*rjny/(rin*rjn))*dindjnvln + nloop3yy
c
        nloop3zz = (rinz*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijz)/(rin*rij))*dindijvli +
     a       (rilz*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijz*rjnz/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjnz/(rin*rjn))*dindjnvln + nloop3zz
c     
        nloop3xy = (rinx*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilx*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinx*(-rijy)/(rin*rij))*dindijvli +
     a       (rilx*(-rljy)/(ril*rlj))*dildjlvln +
     a       (rijx*rjny/(rij*rjn))*dijdjnvlj +
     a       (rinx*rjny/(rin*rjn))*dindjnvln + nloop3xy
c
        nloop3xz = (rinx*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilx*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinx*(-rijz)/(rin*rij))*dindijvli +
     a       (rilx*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijx*rjnz/(rij*rjn))*dijdjnvlj +
     a       (rinx*rjnz/(rin*rjn))*dindjnvln + nloop3xz
c
        nloop3yx = (riny*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rily*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (riny*(-rijx)/(rin*rij))*dindijvli +
     a       (rily*(-rljx)/(ril*rlj))*dildjlvln +
     a       (rijy*rjnx/(rij*rjn))*dijdjnvlj +
     a       (riny*rjnx/(rin*rjn))*dindjnvln + nloop3yx
c     
        nloop3yz = (riny*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rily*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (riny*(-rijz)/(rin*rij))*dindijvli +
     a       (rily*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijy*rjnz/(rij*rjn))*dijdjnvlj +
     a       (riny*rjnz/(rin*rjn))*dindjnvln + nloop3yz
c     
        nloop3zx = (rinz*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijx)/(rin*rij))*dindijvli +
     a       (rilz*(-rljx)/(ril*rlj))*dildjlvln +
     a       (rijz*rjnx/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjnx/(rin*rjn))*dindjnvln + nloop3zx
c
        nloop3zy = (rinz*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijy)/(rin*rij))*dindijvli +
     a       (rilz*(-rljy)/(ril*rlj))*dildjlvln +
     a       (rijz*rjny/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjny/(rin*rjn))*dindjnvln + nloop3zy
     
 1750    continue

         dijvli = ff3(lpt,negli,neglj)
         dijvlj = ff3(lpt,neglj,negli)
         dildjlvli = ff12(lpt,negli,neglj)
         dildjlvlj = ff12(lpt,neglj,negli)
         dildijvli = ff13(lpt,negli,neglj)
         dildijvlj = ff23(lpt,neglj,negli)
         dijdjlvli = ff23(lpt,negli,neglj)
         dijdjlvlj = ff13(lpt,neglj,negli)
         dijdijvli = ff33(lpt,negli,neglj)
         dijdijvlj = ff33(lpt,neglj,negli)

         lloop5xx = (rijx*rijx/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rilx*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xx + lloop5xx
c
         lloop5yy = (rijy*rijy/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rily*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yy + lloop5yy
c
         lloop5zz = (rijz*rijz/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rilz*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zz + lloop5zz
c
         lloop5xy = (rijx*rijy/rij3)*(dijvli + dijvlj) +
     a             (rilx*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xy + lloop5xy
c
         lloop5xz = (rijx*rijz/rij3)*(dijvli + dijvlj) +
     a             (rilx*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xz + lloop5xz
c
         lloop5yx = (rijy*rijx/rij3)*(dijvli + dijvlj) +
     a             (rily*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yx + lloop5yx
c
         lloop5yz = (rijy*rijz/rij3)*(dijvli + dijvlj) +
     a             (rily*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yz + lloop5yz
c
         lloop5zx = (rijz*rijx/rij3)*(dijvli + dijvlj) +
     a             (rilz*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zx + lloop5zx
c
         lloop5zy = (rijz*rijy/rij3)*(dijvli + dijvlj) +
     a             (rilz*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zy + lloop5zy
c
 1700  continue
        
        smij = sqrt(pmass(ipt)*pmass(jpt))
      
        ibb = 3*(ipt-1)+1
        jbb = 3*(jpt-1)+1

        DXXij = HALF*(firstxx 
     a       + lloopxx + lloop2xx + lloop5xx)/smij
      
        if (abs(DXXij).gt.TOL) then
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb
        hess(idp,ibb) = DXXij
        endif

        DXYij = HALF*(firstxy 
     a       + lloopxy + lloop2xy + lloop5xy)/smij
        if (abs(DXYij).gt.TOL) then
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+1
        hess(idp,ibb) = DXYij
        endif

        DXZij = HALF*(firstxz 
     a       + lloopxz + lloop2xz + lloop5xz)/smij
        if (abs(DXZij).gt.TOL) then
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+2
        hess(idp,ibb) = DXZij
        endif

        DYXij = HALF*(firstyx
     a       + lloopyx + lloop2yx + lloop5yx)/smij
        if (abs(DYXij).gt.TOL) then
        ibb = ibb+1
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb
        hess(idp,ibb) = DYXij
        ibb = ibb-1
        endif

        DYYij = HALF*(firstyy 
     a      + lloopyy + lloop2yy + lloop5yy)/smij
        if (abs(DYYij).gt.TOL) then
        ibb = ibb+1
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+1
        hess(idp,ibb) = DYYij
        ibb = ibb-1
        endif

        DYZij = HALF*(firstyz 
     a      + lloopyz + lloop2yz + lloop5yz)/smij
        if (abs(DYZij).gt.TOL) then
        ibb = ibb+1
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+2
        hess(idp,ibb) = DYZij
        ibb = ibb-1
        endif

        DZXij = HALF*(firstzx 
     a      + lloopzx + lloop2zx + lloop5zx)/smij
        if (abs(DZXij).gt.TOL) then
        ibb = ibb+2
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb
        hess(idp,ibb) = DZXij
        ibb = ibb-2
        endif

        DZYij = HALF*(firstzy 
     a      + lloopzy + lloop2zy + lloop5zy)/smij
        if (abs(DZYij).gt.TOL) then
        ibb = ibb+2
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+1
        hess(idp,ibb) = DZYij
        ibb = ibb-2
        endif

        DZZij = HALF*(firstzz 
     a    + lloopzz + lloop2zz + lloop5zz)/smij
        if (abs(DZZij).gt.TOL) then
        ibb = ibb+2
        idx(ibb) = idx(ibb)+1
        idp = idx(ibb)
        id(idp,ibb) = jbb+2
        hess(idp,ibb) = DZZij
        ibb = ibb-2
        endif
 1100   continue

C  THIS PART CALCULATES DXX(Ipt,I) DXY(Ipt,I) ETC, IE, J = I

       sumkxx = ZERO
       sumkyy = ZERO
       sumkzz = ZERO
       sumkxy = ZERO
       sumkxz = ZERO
       sumkyx = ZERO
       sumkyz = ZERO
       sumkzx = ZERO
       sumkzy = ZERO

       do 1500 k = 1,INDEX2(ipt)
       kpt = LIST2(ipt,k)
       negki = 0
       do 604 iw = 1,index2(kpt)
 604   if (list2(kpt,iw).eq.ipt) negki = iw
       rik = r(ipt,k)
       rik2 = rik*rik
       rik3 = rik2*rik
       rikx = RX(ipt,k)
       riky = RY(ipt,k)
       rikz = Rz(ipt,k)
       dikvik = ff1(ipt,k)
       dikvki = ff1(kpt,negki)
       dikdikvik = ff11(ipt,k)
       dikdikvki = ff11(kpt,negki)

       firstxx = (1./rik - rikx*rikx/rik3)*(dikvik + dikvki) +
     a      (rikx*rikx/rik2)*(dikdikvik + dikdikvki)
c
       firstyy = (1./rik - riky*riky/rik3)*(dikvik + dikvki) +
     a      (riky*riky/rik2)*(dikdikvik + dikdikvki)
c     
       firstzz = (1./rik - rikz*rikz/rik3)*(dikvik + dikvki) +
     a      (rikz*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstxy = -(rikx*riky/rik3)*(dikvik + dikvki) +
     a      (rikx*riky/rik2)*(dikdikvik + dikdikvki)
c     
       firstxz = -(rikx*rikz/rik3)*(dikvik + dikvki) +
     a      (rikx*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstyx = -(riky*rikx/rik3)*(dikvik + dikvki) +
     a      (riky*rikx/rik2)*(dikdikvik + dikdikvki)
c     
       firstyz = -(riky*rikz/rik3)*(dikvik + dikvki) +
     a      (riky*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstzx = -(rikz*rikx/rik3)*(dikvik + dikvki) +
     a      (rikz*rikx/rik2)*(dikdikvik + dikdikvki)
c
       firstzy = -(rikz*riky/rik3)*(dikvik + dikvki) +
     a      (rikz*riky/rik2)*(dikdikvik + dikdikvki)

      lloop3xx = ZERO
      lloop3yy = ZERO
      lloop3zz = ZERO
      lloop3xy = ZERO
      lloop3xz = ZERO
      lloop3yx = ZERO
      lloop3yz = ZERO
      lloop3zx = ZERO
      lloop3zy = ZERO

      do 1550 l = 1, INDEX2(ipt)
      if (l.eq.k) goto 1550
      ril = r(ipt,l)
      rilx = RX(ipt,l)
      rily = RY(ipt,l)
      rilz = RZ(ipt,l)
      nloop2 = ZERO
      do 1580 m = 1, INDEX2(ipt)
      if ((m.eq.k).or.(m.eq.l)) goto 1580
      dildikvin = pf3(ipt,m)*zet2(ipt,m,l)*zet2(ipt,m,k)
      nloop2 = nloop2 + dildikvin
 1580 continue

      dikvil = ff2(ipt,l,k)
      dikdikvil = ff22(ipt,l,k)
      dikdilvik = ff12(ipt,k,l)
      dikdilvil = ff12(ipt,l,k)
    

      lloop3xx = lloop3xx + (1./rik - rikx*rikx/rik3)*dikvil +
     a     (rikx*rikx/rik2)*dikdikvil +
     a     (rikx*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yy = lloop3yy + (1./rik - riky*riky/rik3)*dikvil +
     a     (riky*riky/rik2)*dikdikvil +
     a     (riky*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zz = lloop3zz + (1./rik - rikz*rikz/rik3)*dikvil +
     a     (rikz*rikz/rik2)*dikdikvil +
     a     (rikz*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3xy = lloop3xy  - (rikx*riky/rik3)*dikvil +
     a     (rikx*riky/rik2)*dikdikvil +
     a     (rikx*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3xz = lloop3xz  - (rikx*rikz/rik3)*dikvil +
     a     (rikx*rikz/rik2)*dikdikvil +
     a     (rikx*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yx = lloop3yx  - (riky*rikx/rik3)*dikvil +
     a     (riky*rikx/rik2)*dikdikvil +
     a     (riky*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yz = lloop3yz  - (riky*rikz/rik3)*dikvil +
     a     (riky*rikz/rik2)*dikdikvil +
     a     (riky*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zx = lloop3zx  - (rikz*rikx/rik3)*dikvil +
     a     (rikz*rikx/rik2)*dikdikvil +
     a     (rikz*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zy = lloop3zy  - (rikz*riky/rik3)*dikvil +
     a     (rikz*riky/rik2)*dikdikvil +
     a     (rikz*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)

 1550 continue

      lloop4xx = ZERO
      lloop4yy = ZERO
      lloop4zz = ZERO
      lloop4xy = ZERO
      lloop4xz = ZERO      
      lloop4yx = ZERO
      lloop4yz = ZERO
      lloop4zx = ZERO
      lloop4zy = ZERO
      
      do 1600 l = 1, INDEX2(kpt)
      lpt = LIST2(kpt,l)
      if(lpt.eq.ipt) goto 1600
      rkl = r(kpt,l)
      rkl2 = rkl*rkl
      rklx = RX(kpt,l)
      rkly = RY(kpt,l)
      rklz = RZ(kpt,l)
      rmuil = (RX(kpt,negki)*RX(kpt,l) + RY(kpt,negki)*RY(kpt,l) +
     a     RZ(kpt,negki)*RZ(kpt,l))/(rik*rkl)
      ril2 = rik2 + rkl2 -2.0*rmuil*rik*rkl
      ril = sqrt(ril2)
      ril3 = ril2*ril
      rilx = -1.0*(RX(kpt,negki) - RX(kpt,l))
      rily = -1.0*(RY(kpt,negki) - RY(kpt,l))
      rilz = -1.0*(RZ(kpt,negki) - RZ(kpt,l))
      
      nloop4xx = ZERO
      nloop4yy = ZERO
      nloop4zz = ZERO
      nloop4xy = ZERO
      nloop4xz = ZERO
      nloop4yx = ZERO
      nloop4yz = ZERO
      nloop4zx = ZERO
      nloop4zy = ZERO 
      
      do 1650 m = 1, INDEX2(kpt)
      npt = LIST2(kpt,m)
      if ((npt.eq.ipt) .or.(m.eq.l)) goto 1650
      rkn = r(kpt,m)
      rkn2 = rkn*rkn
      rmuin = (RX(kpt,negki)*RX(kpt,m) + RY(kpt,negki)*RY(kpt,m) +
     a     RZ(kpt,negki)*RZ(kpt,m))/(rik*rkn)
      rin2 = rik2 + rkn2 - 2.0*rmuin*rik*rkn
      rin = sqrt(rin2)
      rinx = RX(kpt,m)-RX(kpt,negki)
      riny = RY(kpt,m)-RY(kpt,negki)
      rinz = RZ(kpt,m)-RZ(kpt,negki)
      
      dildinvki = pf3(kpt,negki)*zet3(kpt,negki,l)*zet3(kpt,negki,m)
       
      nloop4xx = nloop4xx + (rilx*rinx/(ril*rin))*dildinvki
      nloop4yy = nloop4yy + (rily*riny/(ril*rin))*dildinvki
      nloop4zz = nloop4zz + (rilz*rinz/(ril*rin))*dildinvki
      nloop4xy = nloop4xy + (rilx*riny/(ril*rin))*dildinvki
      nloop4xz = nloop4xz + (rilx*rinz/(ril*rin))*dildinvki
      nloop4yx = nloop4yx + (rily*rinx/(ril*rin))*dildinvki
      nloop4yz = nloop4yz + (rily*rinz/(ril*rin))*dildinvki
      nloop4zx = nloop4zx + (rilz*rinx/(ril*rin))*dildinvki
      nloop4zy = nloop4zy + (rilz*riny/(ril*rin))*dildinvki    
 1650 continue

      dikvkl = ff2(kpt,l,negki)
      dilvki = ff3(kpt,negki,l)
      dilvkl = ff3(kpt,l,negki)
      dildilvkl = ff33(kpt,l,negki)
      dildikvki = ff13(kpt,negki,l)
      dildikvkl = ff23(kpt,l,negki)
      dildilvki = ff33(kpt,negki,l)
      dikdikvkl = ff22(kpt,l,negki)

      lloop4xx = (1./rik - rikx*rikx/rik3)*dikvkl +
     a               (1./ril -rilx*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rilx*rikx + rikx*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*rikx/rik2)*dikdikvkl +
     a   (rilx*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xx + lloop4xx
c
      lloop4yy = (1./rik - riky*riky/rik3)*dikvkl +
     a               (1./ril -rily*rily/ril3)*(dilvki + dilvkl) +
     a   ((rily*riky + riky*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*riky/rik2)*dikdikvkl +
     a   (rily*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yy + lloop4yy
c
      lloop4zz = (1./rik - rikz*rikz/rik3)*dikvkl +
     a               (1./ril -rilz*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rilz*rikz + rikz*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*rikz/rik2)*dikdikvkl +
     a   (rilz*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zz + lloop4zz
c
      lloop4xy = -(rikx*riky/rik3)*dikvkl 
     a              -(rilx*rily/ril3)*(dilvki + dilvkl) +
     a   ((rilx*riky + rikx*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*riky/rik2)*dikdikvkl +
     a   (rilx*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xy + lloop4xy
c          
      lloop4xz = -(rikx*rikz/rik3)*dikvkl 
     a              -(rilx*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rilx*rikz + rikx*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*rikz/rik2)*dikdikvkl +
     a   (rilx*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xz + lloop4xz
c          
      lloop4yx = -(riky*rikx/rik3)*dikvkl 
     a              -(rily*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rily*rikx + riky*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*rikx/rik2)*dikdikvkl +
     a   (rily*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yx + lloop4yx
c          
      lloop4yz = -(riky*rikz/rik3)*dikvkl 
     a              -(rily*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rily*rikz + riky*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*rikz/rik2)*dikdikvkl +
     a   (rily*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yz + lloop4yz
c 
      lloop4zx = -(rikz*rikx/rik3)*dikvkl 
     a              -(rilz*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rilz*rikx + rikz*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*rikx/rik2)*dikdikvkl +
     a   (rilz*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zx + lloop4zx
c          
      lloop4zy = -(rikz*riky/rik3)*dikvkl 
     a              -(rilz*rily/ril3)*(dilvki + dilvkl) +
     a   ((rilz*riky + rikz*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*riky/rik2)*dikdikvkl +
     a   (rilz*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zy + lloop4zy
          
 1600 continue

      sumkxx = sumkxx + firstxx + lloop3xx + lloop4xx
      sumkyy = sumkyy + firstyy + lloop3yy + lloop4yy
      sumkzz = sumkzz + firstzz + lloop3zz + lloop4zz
      sumkxy = sumkxy + firstxy + lloop3xy + lloop4xy
      sumkxz = sumkxz + firstxz + lloop3xz + lloop4xz
      sumkyx = sumkyx + firstyx + lloop3yx + lloop4yx
      sumkyz = sumkyz + firstyz + lloop3yz + lloop4yz
      sumkzx = sumkzx + firstzx + lloop3zx + lloop4zx
      sumkzy = sumkzy + firstzy + lloop3zy + lloop4zy
 1500 continue

      smij = pmass(ipt)

      ibb = 3*(ipt-1)+1
      jbb = ibb

      DXXij = HALF*sumkxx/smij
      DYYij = HALF*sumkyy/smij
      DZZij = HALF*sumkzz/smij
      DXYij = HALF*sumkxy/smij
      DXZij = HALF*sumkxz/smij
      DYXij = HALF*sumkyx/smij
      DYZij = HALF*sumkyz/smij
      DZXij = HALF*sumkzx/smij
      DZYij = HALF*sumkzy/smij

      if (abs(DXXij).gt.TOL) then
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb
      hess(idp,ibb) = DXXij
      endif

      if (abs(DXYij).gt.TOL) then
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+1
      hess(idp,ibb) = DXYij
      endif

      if (abs(DXZij).gt.TOL) then
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+2
      hess(idp,ibb) = DXZij
      endif
      
      if (abs(DYXij).gt.TOL) then
      ibb = ibb+1
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb
      hess(idp,ibb) = DYXij
      ibb = ibb-1
      endif

      if (abs(DYYij).gt.TOL) then
      ibb = ibb+1
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+1
      hess(idp,ibb) = DYYij
      ibb = ibb-1
      endif

      if (abs(DYZij).gt.TOL) then
      ibb = ibb+1
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+2
      hess(idp,ibb) = DYZij
      ibb = ibb-1
      endif

      if (abs(DZXij).gt.TOL) then
      ibb = ibb+2
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb
      hess(idp,ibb) = DZXij
      ibb = ibb-2
      endif
      
      if (abs(DZYij).gt.TOL) then
      ibb = ibb+2
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+1
      hess(idp,ibb) = DZYij
      ibb = ibb-2
      endif

      if (abs(DZZij).gt.TOL) then
      ibb = ibb+2
      idx(ibb) = idx(ibb)+1
      idp = idx(ibb)
      id(idp,ibb) = jbb+2
      hess(idp,ibb) = DZZij
      ibb = ibb-2
      endif
 1002 continue

      DO I = 1,3*N
      IF (IDX(I).GT.ND) THEN 
      PRINT*,'Error (tersoff_dynamical): idx(',I,')=',idx(I),
     A       ' out of bound =',ND
      PRINT*,'Increase _ND in mulcp.fh and recompile.'
      STOP
      ENDIF
      ENDDO
      
C     print *,'INPUT ATOM NUMBER I and J'
c     read *, i,j
c     i = 1
c     j = 2
c     print *,idx(i),idx(j)
c     do 456 k=1,idx(i)
c     print *,hess(k,i)
c     if (id(k,i).eq.j) print *,'**** hess(k,i)'
c 456   continue

      RETURN
      END
