##########################################################
#   Summary:
# The following example uses fix wall/region to apply an inward harmonic force
#   -2*k*r
# to every particle (group "all") in the simulation, 
# where "r" is the distance from the origin (in this example).
#
# The "k" parameter is an argument of the fix wall/region command.
# (currently set to 1.0e-3 in the example below).
#
# You must set the "rc" and "rc0" parameters (defined below) equal to each other
# and large enough that the particles will never exceed this distance from the
# center of the sphere.  (This must take into account their kinetic/thermal
# energy.  In this example, Uext(rc) = k*rc^2 = 90.0 >> kB*T in most cases.)
#
##########################################################
# Now define "rSphere", a spherical region (used by LAMMPS fix wall/region)
#
#                     sphere_center radius
#                      x0  y0   z0   rc0
#                      |    |    |    |
#                     \|/  \|/  \|/  \|/
#                      V    V    V    V

region rSphere sphere  0    0    0  300.0   side in

# This sets the parameter "rc0" equal to 300.0.
# The "300.0" means we assume no particle will ever be further than 300.0
# from the center of the sphere.  If you want Uext(r) = k*r^2,
# (proportional to r^2), then make sure that "rc0" equals the "rc" parameter.
# ("rc" is the 3rd argument in the "fix wall/region" command below.)



# Now apply an inward force applied to every particle.
# The energy of each particle (due to that force) is:
#
# Uext(r) = k*(rsurf-rc)^2
#       where "rsurf" = the distance from the particle to the surface = (rc0-r)
#              = k*((rc0-rc) - r)^2    (if (rc0-rc) < r < rc0,  0 otherwise)
#
#
#                                               k    ignore   rc
#                                               |      |      |
#                                              \|/    \|/    \|/
#                                               V      V      V

fix fxWall all wall/region rSphere harmonic   1.0e-3  0.0  300.0

# Try playing with the "k" value.


# In the specific case when rc0 = rc, this simplifies to:
#
# Uext(r) = k * r^2     (0 < r < rc0)
#
# This is like a simple spring (rest-length 0) that pulls every particle inward.
# The force is non-zero everywhere except at r=0
#    (and also r>rc0.  However rc0 is very big, so I ignore this possibility.)


# For an explanation of these commands go here:
# http://lammps.sandia.gov/doc/fix_wall_region.html
# http://lammps.sandia.gov/doc/region.html


# ALTERNATE INTERPRETATION:
# This applies an inward harmonic force to every particle.
# If you like, this can be interpreted as a pairwise attractive force between
# all PAIRS of particles.  The energy of interaction between each pair would be:
#        Upair(r) = (k/2N) * (r_i - r_j)^2
# where "k" is the parameter specified above in the "fix wall/region" command.
# To prove this, sum this over all pairs of particles, and assume the average
# particle position is located at the origin.  While not necessarily true,
# we can use fix recenter to enforce this. To do that, uncomment the line below:
#
#fix fxCen all recenter 0 0 0
#
# For more details, see http://en.wikipedia.org/wiki/Radius_of_gyration