!     Evaluation of forces for pair potentials with neighbor list
!     Multiple potentials can be used
!     Potential energy is calculated
!     Frenkel & Smit, p.59
!     MSE 6270, Leonid Zhigilei
      SUBROUTINE F_pair()
      INCLUDE 'common.h'
      LOGICAL :: keySTOP=.FALSE.
      outer_loop: DO I=1,NAN
        ktypei=KTYPE(I)
        inner_loop: DO K=1,NNG(I)
          J=NNNG(k,i)
          DX=XD(1,J)-XD(1,I)
          DY=XD(2,J)-XD(2,I)
          DZ=XD(3,J)-XD(3,I)
!         A more elegant implementation of periodic boundary condition
!         DX=DX-XL*NINT(DX/XL) is time consuming because of the nearest
!         integer function NINT() in FORTRAN.  See the footnote on p. 59
!         of Frenkel & Smit, about the inefficiency of NINT
!         Periodicity *************
          IF(LIDX.EQ.1) THEN
            IF(DX.GT.XLHALF) DX=DX-XL
            IF(DX.LT.-XLHALF) DX=DX+XL
          ENDIF
          IF(LIDZ.EQ.1) THEN
            IF(DZ.GT.ZLHALF) DZ=DZ-ZL
            IF(DZ.LT.-ZLHALF) DZ=DZ+ZL
          ENDIF
          IF(LIDY.EQ.1) THEN
            IF(DY.GT.YLHALF) DY=DY-YL
            IF(DY.LT.-YLHALF) DY=DY+YL
          ENDIF
!         End of periodicity ******
          IJ = ijindex(KTYPEI, KTYPE(J))
          RR=DX*DX+DY*DY+DZ*DZ
          RD=SQRT(RR)
          IF (RD.GE.RM(IJ)) CYCLE inner_loop

!         If your run crashes, you can try to uncomment this part
          IF (RD.LE.RM(IJ)/NT) THEN
            Print *, 'Molecules ',I,' and ',J,' are too close.', &
                     'Program will stop.'
            keySTOP=.TRUE.
          ENDIF

          P=RD*DXR(IJ)
          KF=P
          P=P-KF

!         Potential energy (give half the potential to each particle)
          PENER=.5*(UT(IJ,KF)+(UT(IJ,KF+1)-UT(IJ,KF))*P)
          POT(I)=POT(I)+PENER
          POT(J)=POT(J)+PENER
!         End of potential energy

          P=FT(IJ,KF)+(FT(IJ,KF+1)-FT(IJ,KF))*P
          P=P/RD
          DX=P*DX                          ! x component of force
          FD(1,I)=FD(1,I)-DX
          FD(1,J)=FD(1,J)+DX
          DZ=P*DZ                          ! z component of force
          FD(3,I)=FD(3,I)-DZ
          FD(3,J)=FD(3,J)+DZ
          DY=P*DY                          ! y component of force
          FD(2,I)=FD(2,I)-DY
          FD(2,J)=FD(2,J)+DY
        End do inner_loop
      End do outer_loop
      if(keySTOP) STOP
      RETURN
      END