F77 Code Fragment for Lennard-Jones Forces

       inputs: ndim,natoms,r,ell2,ell,rcut2,sigma,epsilon,potcut
       output: potential,force

	real force(ndim,natoms),r(ndim,natoms),dx(ndim)
        real potential,ell2(ndim),ell(ndim),rcut2,sigma,epsilon,potcut
        integer ndim,natoms,k,i,j
        real r2,r2i,r6i,rforce

c Computes total potential and forces.
c Zero out forces and potential.
	potential = 0
	do k=1,ndim
	   do i=1,natoms
	      force(k,i) = 0
	   enddo
	enddo

c Loop over all pairs of atoms.
	do i=2,natoms
	   do j=1,i-1
c Compute minimum distance between i and j.
	      r2 = 0
	      do k=1,ndim
		 dx(k) = r(k,i) - r(k,j)
c Periodic boundary conditions.
		 if(dx(k).gt.ell2(k))  dx(k) = dx(k) - ell(k)
		 if(dx(k).lt.-ell2(k)) dx(k) = dx(k) + ell(k)
	         r2 = r2 + dx(k)*dx(k)
	      enddo
c Only compute for pairs inside radial cutoff.
	      if(r2.lt.rcut2) then
		 r2i=sigma2/r2
		 r6i=r2i*r2i*r2i
c Shifted Lennard-Jones potential.
		 potential = potential+epsilon*ri6*(ri6-1)- potcut
c Radial force.
		 rforce = epsilonr*r6i*r2i*(2*r6i-1)
		 do k = 1 , ndim
		    force(k,i)=force(k,i) + rforce*dx(k)
		    force(k,j)=force(k,j) - rforce*dx(k)
		 enddo
	      endif
	   enddo
	enddo