Re: [ESPResSo-devel] harmonic potential

[Markus Deserno]

Dear all,

I'm not sure whether this email gets through to the ESPResSo list;
if it doesn't, would please someone relay?

> I vote for Axel's second solution.
> Default cutoff 2*r_0 and negative cutoff means no cutoff.

I vote against (sorry, Axel!), because I think a default cutoff of
2*r_0 does not make physical sense at all.  Let us think a bit more
carefully what is actually going on here.  I think we're talking about
potentials of the form

   U(r) = (1/2) * k * (r - r_0)^2

Evidently, there is nothing special about the point r = 2*r_0.

What we really want is to avoid bonds stretching so much that their
ENERGY becomes unphysical.  I think therefore that an energy criterion
is much more natural.  Let us say that we permit a maximum energy
E_max before the bond fails.  This then easily translates to a cutoff
at which this happens:

   (1/2) * k * (r_cut - r_0)^2 = E_max

implies

   r_cut = r_0 + sqrt[ 2 * E_max / k ]

You can either ask the user to supply a value of E_max, or you can
guess a value which will unlikely happen.  As you know, from StatPhys
the average energy of such a harmonic bond is 1/2 kT.  Choosing E_max
= 10kT might thus be good, but I suspect it would still break way too
often just by accident (because exp(-10)~4.5E-5).  Choosing E_max=20kT
is probably safe, since exp(-20)~2E-9 and will thus not happen during
a simulation.  The main problem with this is that for soft springs
this will lead to a quite substantial increase of the cutoff radius.
Let me say a few more words about that.

The latter problem does not seem to permit a simple numerical
solution.  Harmonic bonds are notorious for being infinitely
stretchable (ask all the Gaussian Polymer Field Theorists!).  Rather
than becoming algorithmically creative, my suggestion would be to
reconsider the model which required the harmonic bond.  Evidently
stiff harmonic bonds are not much of a problem because they do not
force one to increase the cutoff very much.  Soft harmonic bonds are
the issue.  But then, is it really physically necessary to have such a
soft bond, or is it just pure laziness?  Specifically, would it be
physically equally acceptable to have a bond of the form

  (1/2) * k * (r - r_0)^2 + (1/4) * k' * (r - r_0)^4

such that for short and medium elongations the quartic term does not
matter, but for elongations getting into the "many kT range" the
quartic term kicks in and begins to suppress the occurrence of such a
large elongation?  Evidently, many other ways are possible for doing
something like that.

What I'm saying is that my suggested energy-criterion cutoff makes
physical sense, but might yield uncomfortably large cutoffs for soft
springs.  I then push the burden back to the modeler and ask, whether
it really has to be that soft all the way out.

ESPResSo could use the standard criterion E_max = 20kT (or anything
else which a bit more thinking seems to agree with), from there work
out what the cutoff is, and then tell the user the answer.  If this
gets too large by some other fluffy standards, ESPResSo might issue a
warning and recommend that the user artifically stiffen the bond at
large elongations in order to rescue numerical efficiency.

Best,

Markus

-- 
Dr. Markus Deserno
Associate Professor of Physics    ++1-412-268-4401 (office)
Carnegie Mellon University        ++1-412-681-0648 (fax)
5000 Forbes Avenue                ++1-412-268-8367 (Donna Thomas)
Pittsburgh, PA 15213              address@hidden