!     Evaluation of forces for the Stillinger Weber potential for Silicon
!     [F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262-5271 (1985)]
!     MSE 6270, Leonid Zhigilei and Avinash Dongare

      SUBROUTINE F_SW()
      INCLUDE 'common.h'
      LOGICAL :: keySTOP=.FALSE.
      REAL *8 LAM

      outer_loop: DO I=1,NAN        !loop over ALL particles
        ktypei=KTYPE(I)

        inner_loop: DO KJ=1,NNG(I)  !loop over known possible neighbors of i
          J=NNNG(KJ,i)              !extract the atom number
          KTYPEJ = KTYPE(J)

          DXIJ=XD(1,J)-XD(1,I)      ! distances Xij, Yij and Zij
          DYIJ=XD(2,J)-XD(2,I)
          DZIJ=XD(3,J)-XD(3,I)

!         Periodicity *************
          IF(LIDX.EQ.1) THEN
            IF(DXIJ.GT.XLHALF) DXIJ=DXIJ-XL
            IF(DXIJ.LT.-XLHALF) DXIJ=DXIJ+XL
          ENDIF
          IF(LIDZ.EQ.1) THEN
            IF(DZIJ.GT.ZLHALF) DZIJ=DZIJ-ZL
            IF(DZIJ.LT.-ZLHALF) DZIJ=DZIJ+ZL
          ENDIF
            IF(LIDY.EQ.1) THEN
            IF(DYIJ.GT.YLHALF) DYIJ=DYIJ-YL
            IF(DYIJ.LT.-YLHALF) DYIJ=DYIJ+YL
          ENDIF
!         End of periodicity ******

          IJ = ijindex(KTYPEI, KTYPE(J))     ! potential type for the pair
          RRIJ=DXIJ*DXIJ+DYIJ*DYIJ+DZIJ*DZIJ
          RDIJ=SQRT(RRIJ)                    ! Distance Rij

          IF (RDIJ.GE.RM(IJ)) CYCLE inner_loop

!         We are using full neighbor list.
!         Only do pair potential if j>i

          JGTI: IF (J.GT.i) THEN
            P=RDIJ*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/RDIJ
            DX=P*DXIJ                        ! x component of force
            FD(1,I)=FD(1,I)-DX
            FD(1,J)=FD(1,J)+DX
            DZ=P*DZIJ                        ! z component of force
            FD(3,I)=FD(3,I)-DZ
            FD(3,J)=FD(3,J)+DZ
            DY=P*DYIJ                        ! y component of force
            FD(2,I)=FD(2,I)-DY
            FD(2,J)=FD(2,J)+DY

          END IF JGTI

!         Now we add the three body part due to the neighbors of i

          inner_loop2: DO KK=KJ+1,NNG(I)  ! Loop over all neighbors of atom I

            K=NNNG(KK,I)
            KTYPEK = KTYPE(K)
            DXIK=XD(1,K)-XD(1,I)
            DYIK=XD(2,K)-XD(2,I)
            DZIK=XD(3,K)-XD(3,I)
!           Periodicity *************
            IF(LIDX.EQ.1) THEN
              IF(DXIK.GT.XLHALF) DXIK=DXIK-XL
              IF(DXIK.LT.-XLHALF) DXIK=DXIK+XL
            ENDIF
            IF(LIDZ.EQ.1) THEN
              IF(DZIK.GT.ZLHALF) DZIK=DZIK-ZL
              IF(DZIK.LT.-ZLHALF) DZIK=DZIK+ZL
            ENDIF
            IF(LIDY.EQ.1) THEN
              IF(DYIK.GT.YLHALF) DYIK=DYIK-YL
              IF(DYIK.LT.-YLHALF) DYIK=DYIK+YL
            ENDIF
!           End of periodicity ******

            IK = ijindex(KTYPEI, KTYPE(K))
            RRIK=DXIK*DXIK+DYIK*DYIK+DZIK*DZIK
            RDIK=SQRT(RRIK)
            IF (RDIK.GE.RM(IK)) CYCLE inner_loop2

!           Calculations of the distance between J and K atoms
            DXJK=XD(1,K)-XD(1,J) ! with atom I. Hence K>J required
            DYJK=XD(2,K)-XD(2,J)
            DZJK=XD(3,K)-XD(3,J) ! distances Xik, Yik and Zik

!           Periodicity *************
            IF(LIDX.EQ.1) THEN
              IF(DXJK.GT.XLHALF) DXJK=DXJK-XL
              IF(DXJK.LT.-XLHALF) DXJK=DXJK+XL
            ENDIF
            IF(LIDZ.EQ.1) THEN
              IF(DZJK.GT.ZLHALF) DZJK=DZJK-ZL
              IF(DZJK.LT.-ZLHALF) DZJK=DZJK+ZL
            ENDIF
            IF(LIDY.EQ.1) THEN
              IF(DYJK.GT.YLHALF) DYJK=DYJK-YL
              IF(DYJK.LT.-YLHALF) DYJK=DYJK+YL
            ENDIF
!           End of periodicity ******

            RRJK = DXJK*DXJK + DYJK*DYJK + DZJK*DZJK
            RDJK = SQRT(RRJK)
            JK = ijindex(KTYPE(J), KTYPE(K))

!           This now forms the triplet of I, J and K atoms
            eij = SW_SIG(IJ)/(RDIJ - RM(IJ))
            eik = SW_SIG(IK)/(RDIK - RM(IK))

            deij = -SW_SIG(IJ)/(RDIJ-RM(IJ))**2
            deik = -SW_SIG(IK)/(RDIK-RM(IK))**2

 !          Use the cosine rule to calculate cos(theta)
            costh = (RRIK + RRIJ - RRJK)/(2.0D+00*RDIJ*RDIK)

            dcosdrij = 1.0D+00/RDIK - costh/RDIJ
            dcosdrik = 1.0D+00/RDIJ - costh/RDIK
            dcosdrjk = -RDJK/(RDIJ * RDIK)

            cosTHIRD = costh + 1.0D+00/3.0D+00
            IF(cosTHIRD ==0.0D+00) THEN
              TWcTH = 0.0D+00
            ELSE
              TWcTH = 2.0D+00/cosTHIRD
            ENDIF

            EPS = SQRT(SW_EPS(IJ) * SW_EPS(IK))
            LAM = (SW_LAM(KTYPEJ)*((SW_LAM(KTYPEI))**2.0d0)*SW_LAM(KTYPEK))**0.25d0
            potegy= EPS*LAM*exp(SW_GAM(IJ)*eij+SW_GAM(IJ)*eik)*cosTHIRD**2


            PEGY   = potegy/ENUNIT
!           Arbitarily giving all the three body energy to the central atom
            POT(I)= POT(I)+PEGY

!           Calculate the Forces

!           Rij derivatives:
            FFIJ = -potegy * (SW_GAM(IJ)*deij + TWcTH * dcosdrij)/RDIJ
            FFIJ =FFIJ/ENUNIT               ! Program units
            FD(1,I) = FD(1,I) - FFIJ*DXIJ
            FD(1,J) = FD(1,J) + FFIJ*DXIJ

            FD(3,I) = FD(3,I) - FFIJ*DZIJ
            FD(3,J) = FD(3,J) + FFIJ*DZIJ

            FD(2,I) = FD(2,I) - FFIJ*DYIJ
            FD(2,J) = FD(2,J) + FFIJ*DYIJ

!           Rik derivatives:
            FFIK = -potegy * (SW_GAM(IK)*deik + TWcTH * dcosdrik)/RDIK
            FFIK =FFIK/ENUNIT               ! Program units
            FD(1,I) = FD(1,I) - FFIK*DXIK
            FD(1,K) = FD(1,K) + FFIK*DXIK

            FD(3,I) = FD(3,I) - FFIK*DZIK
            FD(3,K) = FD(3,K) + FFIK*DZIK

            FD(2,I) = FD(2,I) - FFIK*DYIK
            FD(2,K) = FD(2,K) + FFIK*DYIK

!           RJk derivatives:
            FFJK =  -potegy * (TWcTH * dcosdrjk)/RDJK
            FFJK =FFJK/ENUNIT               ! Program units
            FD(1,J) = FD(1,J) - FFJK*DXJK
            FD(1,K) = FD(1,K) + FFJK*DXJK

            FD(3,J) = FD(3,J) - FFJK*DZJK
            FD(3,K) = FD(3,K) + FFJK*DZJK

            FD(2,J) = FD(2,J) - FFJK*DYJK
            FD(2,K) = FD(2,K) + FFJK*DYJK

          END DO inner_loop2
        END DO inner_loop
      END DO outer_loop
      RETURN
      END SUBROUTINE F_SW