Hi all.
I was wondering about which quantity is supposed to be conserved in NPH dynamics.
One piece of the total energy of the system is computed and available as a “temp” in older versions of LAMMPS and available as a “compute” in current versions.
Hi again.
Let me rephrase my question. I have now upgraded to the latest LAMMPS release which has an “enthalpy” option in the custom thermo output. I was expecting this “enthalpy” to be conserved under the NPH dynamics.
As a test, I modified the Lennard-Jones melt example, “in.melt” with some changes. Basically, I just let the crystal contract under the Lennard-Jones attraction and then rebound.
The system goes unstable at the end of the first cycle of contraction/expansion. This doesn’t seem unreasonable. The thing I’m worried about is that the enthalpy is not even approximately conserved. Does anyone know why?
Thanks,
Craig
--------------------input file --------------------------
3d Lennard-Jones melt
units lj
atom_style atomic
#lattice fcc 0.8442
lattice fcc 0.83
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1
mass 1 1.0
#velocity all create 3.0 87287
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix nphFix all nph aniso 0.0 0.0 NULL NULL NULL NULL 10.0 drag 0.0
thermo_style custom step ke pe etotal enthalpy pxx lx
thermo_modify flush yes
dump id all atom 10 dump.melt
thermo 10
run 50000
There is also a fix_modify thermo option for fix nph which will
print a quantity (or add it to the pot energy) that might be closer
to what is conserved with NPH (or NPT).
You can see the formula for it at the bottom of fix_nph.cpp
Steve