Dear LAMMPS users,
Conceptually, is the use of box/relax equivalent to enthalpy minimization? From how I understand fix box/relax, it allows you to minimize energy while allowing V as a degree of freedom, and constraining P. Enthalpy minimization, on the other hand, minimizes H = U + PV, so that the condition dH/dV = dU/dV + P = 0 sets the pressure, dH/dr_i = - forces, H ~= U for incompressible liquid. What I don’t understand is how box/relax constrains the pressure. With enthalpy minimization, dH/dL can be derived analytically or approximated with finite differences.
I have run into some trouble using box/relax because the line search alpha -> 0 with default minimization settings far before a minimum is reached. Attempting to use quadratic line search led to increases in U, which isn’t acceptable. In order to get a decent minimum with box/relax, I wrote a script (following the note in the manual) which alternately minimizes with box/relax to decent energy tolerance (1e-8), then minimizes constant volume for 1 step. After alternating this several hundred times, I can finally reach a state where minimizing for both box/relax and constant volume get zero line search alpha. Why does box/relax lose steam before reaching a minimum, and has anyone found a better way to overcome this problem?
Thank you,
Harold Hatch
Aidan can answer this better, but the issue is that for various kinds of
box relaxations there is no formulation (mathematically) of the min problem
that satisfies all the constraints exactly at the energy min. Except when
your starting point is the right answer. So iterating, as you are doing, is
not a bad approach, if you want a hi-precision answer.
Steve
Dear Steve,
Thank you for your informative and timely reply. I see that the analytical
form of dH/dL can get very tricky, and I was only considering a cubic box.
I was also only considering isotropic pressure. So box/relax must be easier
to implement for more general box relaxations. Although dH/dL could be
approximated by finite difference from calculating enthalpy change upon a
small contraction/expansion, which may be slower. But Middleton and Wales
(JCP vol 118, pg 4583, 2003) say the analytic derivatives didn't speed-up
the calculation.
Still, I do not understand how the pressure is constrained with box/relax,
and I eagerly await any help Dr. Thompson may offer.
Regards,
Harold Hatch
P.S. I haven't forgotten about our discussion of the stress tensor, which I
am still working out. Thank you for the help you both have offered in this
regard.
Harold,
For the xyz style, which you are using, the LAMMPS formulation corresponds
to the enthalpy-minimization equations you wrote in your original e-mail.
So, the minimize command will try to converge to a local energy minimum with
the following properties:
dU/dr = -F = 0
dU/dV = -p = p_target
Because the linesearch prohibits uphill moves, this state will also be a
minimum of U+p_target*V.
I am not surprised that the default linesearch style failed for you.
However, I think you were too quick to reject the quadratic linesearch
style. You may need to tryo harder and play with the parameters,
particularly vmax.
Attempting to use quadratic line search led to increases in U, which isn¹t
acceptable.
The potential energy may increase, but this is to be expected, if you are
converging to a high-pressure state. I have found that for many crystalline
systems, the quadratic linesearch works very well.
Aidan
Dear Aidan,
Thank you for clarifying box/relax. I have realized that the 'offending'
energy increases I mentioned still led to enthalpy decreases, because the
energy changes were on the order of delU ~ p_target*dV/N. You obviously
understood this better than I did, so thank you for pointing it out. I
erroneously considered the instantaneous pressure in the calculation of
enthalpy instead of the target pressure when I decided quadratic line search
doesn't work. But the instantaneous pressure is only the thermodynamic
pressure when appropriately ensemble averaged.
Quadratic line search works quite well!
Best,
Harold
Dear All,
I've tested fix box/relax as carefully as I can and here is a summary of
what I have found.
For constant pressure minimization (enthalpy minimization) of binary lennard
jones system, minima found by lammps precisely agree with those cited by
Middleton and Wales (JCP vol 118, pg 4583, 2003). I wrote the style for the
potential for the peculiar cutoff that I can share should anyone need it
(U=U^LJ + A*r^2 + B; A,B chosen such that U,U'->0 as r->rc). Quadratic line
search was very efficient, with good convergence, and there were no uphill
moves in enthalpy. Awesome!
For enthalpy minimization of a protein in water using real units (no shake),
however, I observed uphill moves in enthalpy (see log.real attached). If I
let the simulation run longer, where it is closer to converging, I calculate
a series of enthalpies (perhaps erroneously?) as follows (see log.real2).
-52251.45904
-52251.46304
-52251.45505
-52251.44505
-52251.23311
-52251.41006
-52251.39606
-52251.40806
-52251.39406
-52251.40306
-52251.39106
-52251.39806
-52251.38707
-52251.39406
-52251.37807
-52251.39006
-52251.37607
-52251.38507
-52251.36107
-52251.38107
-52251.35707
where pv2e = 1.458397427E-05, H = U + pv2e*P_target*V (extensive)
I realize that most users simply use minimization to remove overlaps and
such, but for those few who wish to generate inherent structures and get
local minimum precisely, one should be very concerned by uphill moves in
enthalpy (if minimizing at constant pressure). There is also the acute
possibility that I did something wrong, and there are no uphill moves in
enthalpy. I realize these tests are very tedious, so I would be very
appreciative if someone checks these calculations. Please let me know if
you need any of the data/input files.
On a good note, the enthalpy minimum which the simulation appears to be
converging toward agrees qualitatively with a simulation program written in
our lab for a rigid protein (but no exact quantitative comparison between
these two programs seems to be possible).
Regards,
Harold
log.real (2.83 KB)
log.real2 (3.42 KB)
Glad to here that you were able to reproduce Wales' LJ results. That is a
very stringent test. I am guessing that for the water calculation, these
small uphill moves are due to the approximate method used to calculate the
electrostatic energy and forces. If you bump up the kspace precision, you
should see a decrease in the magnitude of the enthalpy increases.
Aidan
I was concerned with the accuracy of the kspace as well, so I had tried
this, but it actually led to greater enthalpy increases. For kspace_style
ewald 1e-7;
Step TotEng enthalpy Press Volume
0 -39946.12 -39944.533 -4262.1338 108825.95
1 -39955.695 -39954.117 123.95047 108194.39
2 -43137.037 -43135.462 1501.5945 108000.72
3 -43144.638 -43143.059 -44.479202 108221.29
4 -44270.145 -44268.571 -173.5008 107939.02
5 -44267.816 -44266.243 626.37408 107824.75
In theory once the potential is specified, however, regardless of how well
it models electrostatics it shouldn't lead to inaccuracy in calculating the
minimum. Furthermore, even with the protein-water system, I don't see
enthalpy increases with backtrack line search nor energy increases without
fix/box. This hints that its a fix/box problem and not an ewald problem.
I think I saw a line in the code to reject uphill moves, so why is lammps
allowing these?
I thought another source of inaccuracy could be from using real units. So I
changed my protein system to lj units (with the same data file that is
supposed to be in real units) and no longer see enthalpy increases! To test
the robustness, I kept increasing P_target to 10 (see log.lj10) and still no
enthalpy increases. But at P_target = 100, the enthalpy spikes at the first
step (see log.lj100). I'm not too concerned by this, in fact, the system
seems to prefer lj units. So perhaps there is a unit conversion inaccuracy
(most likely by me). It's possible that the enthalpy I am checking for
uphill moves isn't the same enthalpy that the program is checking for uphill
moves. Indeed, I was toying with the source code and found that this check
(for the energy, with a pdV term "eng") was passing even though the enthalpy
appears to be going up.
Harold
log.lj10 (2.79 KB)
log.lj100 (2.8 KB)
Okay, so if things work right for LJ units, then the apparent enthalpy
increase you observed is probably just a slight disagreement in how you and
LAMMPS convert from PV units to energy units.
The actual energy compute by fix box/relax is:
eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*vol0;
and
force->nktv2p = 68568.415;
Note that small differences in the value of V will also shift the enthalpy
around. If you are using V from the log file, it is not exactly the same as
V stored internally. I have not tried this, but you may be able to get the
true enthalpy used by the minimizer by reporting the LAMMPS thermo quantity
"enthalpy". Don't forget to turn off the kinetic part, or set the velocities
to zero.
We plan on making the fix/box relax energy available to users as a scalar
quantity referenced as say f_myfix, but this is not in the current version
of LAMMPS.
Aidan
Dear Aidan,
regarding enthalpy minimization, we know that
H=E_pot + E_kin + (P_virial + P_kin)V = H_pot + H_kin
E_kin is a constant and does not play any role in minimization of H. But we cannot forget the kinetic part of pressure and just minimize H_pot (box/relax does this exactly).
in other words,
for H_min we should have,
dE_pot/dr = -F = 0
dE_pot/dV = - (P_virial + P_kin).
if you plot H_pot during minimization, you’ll see that it always decreases. However H is not monotonically decreasing (see the attached plot for a soft sphere system.) and as a result the converged configuration is not necessarily in local enthalpy minimum.
adding an option to the box/relax fix to optimize total pressure instead of virial pressure would be an easy way to have enthalpy minimizer (maybe not the best way).
Best,
Majid.
H.pdf (29.9 KB)
Dear Majid,
From what I understand of minimization, velocities play no explicit role and there are no kinetic terms. When minimizing an ensemble of configurations, however, velocities (or temperature) may be important in generating these configurations which will be minimized. But during minimization itself, there are no velocities.
Regards,
Harold
Majid,
It would be pretty easy to do what you are requesting, but I am having a
hard time seeing why it is useful. Also, I am confident that it including
that as an option would make the documentation more confusing to novice
users. Without getting in to a big discussion about the "true" definition of
enthalpy for a minimized structure (I doubt that argument will ever be
settled), can you explain why it is useful to include the kinetic part of
the pressure in fix box/relax.
Aidan
Dear all,
adding the possibility to be able to have total pressure in box/relax fix is very useful, since it enables the user to minimize the true enthalpy (H=E_total + V*P_total) and to study enthalpy landscape of the system.
Currently it is only possible to minimize the configurational part of enthalpy via box/relax (configurational is not a good phrase here. because in enthalpy, kinetic part of pressure and volume are coupled). The resulting configuration does not always corresponds to the minimum enthalpy. So one has to consider the kinetic part in optimization procedure which involves volume changes in addition to particle displacements.
Best,
Majid.
Since enthalpy is given by the total pressure (H=E+PV). So in order to be able to minimize
This is not a compelling argument, because you are requesting that we
combine two things that are not mutually consistent. The two things are:
-Hconf: Configurational enthalpy (depends on x,V)
-Hkin : Kinetic enthalpy (depends on T or velocities, with trivial V
scaling)
These are not consistent because:
Hkin is based on some set of velocities specified by the user, presumably
from some finite-temperature simulation. The minimizer and fix box/relax
have no control over these velocities.
Hconf is for a geometry that is a local potential energy minimum w.r.t x. It
could also be a higher-order stationary point, but this seems unlikely. In
any case, this is not a configuration that is representative of any finite
temperature state. It really is a zero-temperature state. Therefore, the
value of Hconf obtained is not a good estimate of the average Hconf that
would be measured in a finite temperature simulation at the same volume. In
some cases e.g. harmonic soid, average Hconf may be insensitive to
temperature, but there will be other cases where the temperature effect will
be large.
Dear Aidan,
Thanks for the clarification.
Now I see your point. you expect to have no contribution from kinetic
part in the final enthalpy just because that configuration corresponds
to T=0. So the pressure in that state is totally structural and should
be equal to total pressure of the the finite-temperature initial
configuration (both kinetic and virial parts). The starting
configuration is a representative of the NPT ensemble with a fixed
average enthalpy H(T) and of course we should have H_min < H(T).
Majid.