Thanks for your quick reply.
I have to do it numerically for two reasons:2. My final goal is to calculate surface tension and I don't think there is
a built-in command in Lammps to do that. Reproducing the Pz was just and
example and also good practice before I tackle surface tension. To get
surface tension, I should calculate dU/dA where A is the cross-section area
of the simulation box while I keep the volume constant.Have you got a suggestion for my problem yet?
Calculating solvation energy or surface tension is tricky because it
involves calculating a free energy difference between two different
systems with different hamiltonians. There are probably multiple ways
to do this.
One way to do this is to run several simulations with interactions
between solvent and solute turned on and off. (Note: when
interactions are off, the solvent and solute atoms may pass through
each other. When on, they won't.)
Then you would use LAMMPS's new "rerun" command
http://lammps.sandia.gov/doc/rerun.html
...to take these two existing trajectories (on and off) and
recalculate the energy of every frame in the trajectory using the two
different hamiltonians (the two different ettings with the
interactions between solvent switched on and off). You will generate
4 lists of energies from the two trajectories:
The energy of systemON evaluated using hamiltonianON at time t
The energy of systemON evaluated using hamiltonianOFF at time t
The energy of systemOFF evaluated using hamiltonianON at time t
The energy of systemOFF evaluated using hamiltonianOFF at time t
Then you would use a free-energy calculation program (such as the
Weighted Histogram Analysis Method "WHAM" algorithm) to calculate the
free-energy difference between the system with the interactions
switched on and the interactions switched off. Here's a paper
describing WHAM:
http://jcp.aip.org/resource/1/jcpsa6/v129/i12/p124105_s1
Most implementations of WHAM can not handle this kind of energy data,
however Shirt's and Chodera's MBAR probably can:
https://simtk.org/home/pymbar
--- Details ---
In the third case above:
"The energy of systemOFF evaluated using hamiltonianON"
...you may want to place a "cap" on the maximum energy of a
conformation. Atoms which overlap in non-interacting simulations may
have huge energies or NANs when the interaction is switched back on.
This is not supposed to be a problem If you use a robust WHAM
analysis program to analyze your energies., these conformations should
receive a zero statistical weight which should not influence the
result, however some implementations of WHAM crash when they see this.
--- NVT or NPT ---
I can't remember if you can get away with running these simulations at
NVT. At some point I thought you could, but now I have doubt. If you
decide to run these simulations at NPT, then I think you want to
weight the enthalpies not the energies. So I think you would want to
add "P*V(t)" to each energy in the 4 lists above before you run WHAM.
(When using WHAM, remember never to use the Berendsen thermostat for
NPT or NVT, and do not to use the "drag" keyword with fix npt.)
Hope this helps.
Andrew