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.