      SUBROUTINE FORCE

      IMPLICIT CHARACTER ( A-Z )

      COMMON/COORS/RX,RY,RZ
      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ
      COMMON/ENERGY/V,K,W
      COMMON/ATOMS/NATM
      COMMON/LJPAR/SIGMA
      COMMON/ESHIFT/VRCUT,DVRCUT,DVRC12
      COMMON/CUTOFF/RCUT
      COMMON/FORCES/FX,FY,FZ

C    *******************************************************************
C    ** ROUTINE TO COMPUTE FORCS AND POTENTIAL USING A LINK LIST     **
C    *******************************************************************

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ

      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL)

      REAL        RX(NMAX), RY(NMAX), RZ(NMAX)
      REAL        FX(NMAX), FY(NMAX), FZ(NMAX)
      REAL        RCUT, SIGMA
      REAL        VRCUT,DVRCUT,DVRC12
      REAL        V,K,W
      REAL        RXI, RYI, RZI, FXIJ, FYIJ, FZIJ, FIJ, RCUTSQ
      REAL        SIGSQ, FXI, FYI, FZI, SR2, SR6, VIJ, WIJ
      REAL        RIJ, RIJSQ, RXIJ, RYIJ, RZIJ
      REAL        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      LOGICAL     EDGEI, EDGEJ

      INTEGER     HEAD(NCELL), LIST(NMAX), MAP(MAPSIZ), MX, MY,MZ
      INTEGER     ICELL, JCELL0, JCELL, I, J, NABOR, NATM

      SIGSQ  = SIGMA**2
      RCUTSQ = RCUT**2

C    ** ZERO FORCS AND POTENTIAL **
      DO 10 I = 1, NATM
        FX(I) = 0.0
        FY(I) = 0.0
        FZ(I) = 0.0
10      CONTINUE

      V = 0.0
      W = 0.0

C    ** LOOP OVER ALL CELLS **

      DO 5000 ICELL = 1 , MX*MY*MZ

        I = HEAD(ICELL)

C**     TEST TO SEE IF CURRENT CELL IS AN EDGE CELL

        EDGEI=(I.GT.NATM)

C       ** LOOP OVER ALL MOLECULES IN THE CELL **

1000    IF ( I .GT. 0 ) THEN
          RXI = RX(I)
          RYI = RY(I)
          RZI = RZ(I)

C**       IF AN EDGE CELL OMIT FORCE CALCULATIONS BETWEEN ATOMS WITHIN
C**       CURRENT CELL

          IF (EDGEI) GOTO 2200
          FXI = FX(I)
          FYI = FY(I)
          FZI = FZ(I)
C          ** LOOP OVER ALL MOLECULES BELOW I IN THE CURRENT CELL **
          J = LIST(I)
2000      IF ( J .GT. 0 ) THEN
            RXIJ  = RXI-RX(J)
            RYIJ  = RYI-RY(J)
            RZIJ  = RZI-RZ(J)
            RIJSQ = RXIJ*RXIJ+RYIJ*RYIJ+RZIJ*RZIJ
            IF (RIJSQ.LT.RCUTSQ) THEN
              RIJ   = SQRT(RIJSQ)
              SR2   = SIGSQ/RIJSQ
              SR6   = SR2*SR2*SR2
              VIJ   = SR6*(SR6-1.0)-VRCUT-DVRC12*(RIJ-RCUT)
              WIJ   = SR6*(SR6-0.5)+DVRCUT*RIJ
              V     = V+VIJ
              W     = W+WIJ
              FIJ   = WIJ/RIJSQ
              FXIJ  = FIJ*RXIJ
              FYIJ  = FIJ*RYIJ
              FZIJ  = FIJ*RZIJ
              FXI   = FXI+FXIJ
              FYI   = FYI+FYIJ
              FZI   = FZI+FZIJ
              FX(J) = FX(J)-FXIJ
              FY(J) = FY(J)-FYIJ
              FZ(J) = FZ(J)-FZIJ
              ENDIF
            J = LIST(J)
            GO TO 2000
            ENDIF
 2200     CONTINUE

C          ** LOOP OVER NEIGHBOURING CELLS **

        JCELL0 = 13*(ICELL-1)
        DO 4000 NABOR = 1, 13

C             ** LOOP OVER ALL MOLECULES IN NEIGHBOURING CELLS **

          JCELL = MAP(JCELL0+NABOR)
          J = HEAD(JCELL)

C**     TEST TO SEE IF NEIGHBOURING CELL IS AN EDGE CELL

          EDGEJ=(J.GT.NATM)

C**       IF BOTH EDGE CELLS OMIT FORCE CALCULATIONS BETWEEN PAIRS
C**       OF EDGE CELLS

          IF ((EDGEI).AND.(EDGEJ)) GOTO 4000
3000      IF (J.NE.0) THEN
            RXIJ  = RXI-RX(J)
            RYIJ  = RYI-RY(J)
            RZIJ  = RZI-RZ(J)
            RIJSQ = RXIJ*RXIJ+RYIJ*RYIJ+RZIJ*RZIJ
            IF (RIJSQ.LT.RCUTSQ) THEN
              RIJ   = SQRT(RIJSQ)
              SR2   = SIGSQ/RIJSQ
              SR6   = SR2*SR2*SR2
              VIJ   = (SR6*(SR6-1.0)-VRCUT-DVRC12*(RIJ-RCUT))*0.5
              WIJ   = SR6*(SR6-0.5)+DVRCUT*RIJ
              FIJ   = WIJ/RIJSQ
              WIJ   = 0.5*WIJ
              FXIJ  = FIJ*RXIJ
              FYIJ  = FIJ*RYIJ
              FZIJ  = FIJ*RZIJ

C**   IF CURRENT CELL IS NOT AN EDGE CELL SUM FORCES AND POTENTIAL

              IF (.NOT.(EDGEI)) THEN
                V     = V+VIJ
                W     = W+WIJ
                FXI   = FXI+FXIJ
                FYI   = FYI+FYIJ
                FZI   = FZI+FZIJ
                ENDIF

C**   IF NEIGBOURING CELL IS NOT AN EDGE CELL SUM FORCES AND POTENTIAL

              IF (.NOT.(EDGEJ)) THEN
                V     = V+VIJ
                W     = W+WIJ
                FX(J) = FX(J)-FXIJ
                FY(J) = FY(J)-FYIJ
                FZ(J) = FZ(J)-FZIJ
                ENDIF
              ENDIF
            J = LIST(J)
            GO TO 3000
            ENDIF
4000      CONTINUE

C**   IF CURRENT CELL IS NOT AN EDGE CELL RETURN SUMMED FORCES TO ARRAYS

          IF(.NOT.EDGEI)THEN
            FX(I) = FXI
            FY(I) = FYI
            FZ(I) = FZI
            ENDIF
          I = LIST(I)
          GO TO 1000
          ENDIF
5000    CONTINUE

C    ** INCORPORATE ENERGY FACTORS **

      V = 4.0 *V
      W = 48.0*W/3.0

      DO 50 I = 1, NATM
        FX(I) = 48.0*FX(I)
        FY(I) = 48.0*FY(I)
        FZ(I) = 48.0*FZ(I)
50      CONTINUE

      RETURN
      END
