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