Has anyone successfully used LAMMPS to calculate sensible enthalpy and entropy of a substance? I am attempting to calculate sensible enthalpy for a propellant I am researching, but I am getting calculated results that are an order of magnitude off. My simulations predict the condensed-phase density of the substance fairly accurately, but are significantly off on enthalpy.
I am calculating the enthalpy of a substance (using an NPT run) by getting the average E_total (from the LAMMPS output) and adding the average pressure divided by the average density (h = u + P/rho). I do this for two different P-T combinations and calculate the difference in enthalpy between the two runs. This difference is typically almost an order of magnitude (or more) too small as compared with experimental data.
I assume that the units on E_total are in kcal/mol, the units on temperature are in K, the units on pressure are in atmospheres, and the units on volume are in cubic angstroms. In my input file I set the units to “real”.
Any help would be appreciated.
Thanks,
Tim Kokan
Here are some comments from Aidan Thompson at Sandia, who
knows stat mech better than I. His comment about energy units
is correct - for "real" units the reported energy is for all the
atoms in the box (and a mole of the box). I just noticed that
the doc page on thermo_modify is incorrect in stating the defaults
for this. Norm is default "yes" for LJ units, but default "no" for
real units, so energy is not per-atom.
I have never attempted to compute enthalpy differences using MD, but in
should be possible to do it (entropy is much harder to do). The equation
you use seems correct, but remember, LAMMPS reports energy as kcal/mol, but
the mole here refers to moles of the simulation box i.e. kcal/(6.0e23
simulation boxes). There is no scaling by the number of particles in the
box. Hence, I would write:
h = U/N + P/rho
Just a related question to this, is the calculated total energy is always U or it depends upon the ensemble . For example is it Gibbs free energy in NPT.
Thanks
Arnab
Steve Plimpton wrote:
The total energy displayed is always U + KE, regardless
of the ensemble. However, note that fixes nvt and npt have an option
to include a contribution that is added to U. In the case of NPT,
I think it is PV, so that what is then displayed is enthalpy = U + PV,
the quantity that should be conserved as NPT integration proceeds.
The doc for fix_npt explains how to enable this option.
Steve
Steve raises a good point. However, if you use fix npt, then you should
make sure that the energy of the barostat is *not* added to the total
energy. I think this is the default mode, but you should perhaps check by
turning it on and then off. This energy contribution, as well as a
P(V-Vinitial) term, also contains kinetic energy of the barostat, which is
useful for checking energy conservation, but is not relevent to
thermodynamics of the system.
With "energy no" you can then compute enthalpy as:
<H> = <Total Energy + Pressure*Volume>
or maybe
<H> = <Total Energy> + <Pressure><Volume>
Where <> indicates average over the run.
Aidan