The bond contribution to the virial pressure at low density limit

Dear Sir,
I am using the July/19/2011 version of Lammps and I am encountering a problem in a simple simulation of propane molecules. I searched the mailing list but I couldn’t find anything in this respect.
I am running NVT simulation of 500 propane molecules at T=300 K in the framework of TraPPE potential. In the actual TraPPE potential the bond length is kept constant for n-alkanes but in Lammps the SHAKE algorithm (or any other similar methods) cannot be well parallelized, as it is described in the manual, so I used the harmonic potential with an arbitrary force constant for bond vibration. Also, I chose the cutoff distance for each i-i interaction equal to the sigma of the site i so that I can completely ignore the attractive part of the potential and I don’t have any electrostatic interactions.
For each run, I can calculate compressibility factor (Z) with having N, V, T and the pressure calculated by Lammps given by thermo_style command. The problem is even at very low density limit (with the length of the simulation box > 140 A) the calculated Z doesn’t go to 1 and instead, it gets smaller than one. After seeing this problem, I calculated different contribution of virial pressure an I observed that all contribution of virial pressure go to zero at low density limit (as physically expected) but the bond contribution. The bond contribution of virial pressure reaches asymptotically to a negative value (~-10) which cause the Z factor to be less than 1. The same quantity is positive at high density limit which does make sense physically but at low densities it must go to zero as the angle contribution of virial pressure does. I tried several values (high, moderate and low) for bond force constant but I got, again, a negative value for P of bond. I also calculated different components of pressure tensor for bond contribution, the xy, xz and yz components are close to zero but the diagonal components are not.
CH3-CH2 bond length and CH3-CH2-CH3 bond angle are fixed at 1.54 A and 114.0 degress, respectively, based on TraPPE potential. I have monitored the contributions of the total potential seperately as the simulation evolves. The simulation reaches to equilibrium without any problem and the velocity rescaling method (in conjunction with a fix NVE) can control the temperature just fine.
What is wrong with this system? Why the bond contribution of virial pressure gets to a negative number?
Thank you in advance for your help.

I don't see why the virial contribution from bonds should go to 0
at the low-density limit. The bonded atoms are not getting farther
apart, so however much virial contribution is in one bond will
just stay constant and get divided by the growing volume
as the box size grows.

Maybe Aidan wants to comment.

Steve

Thank you for your fast reply.
Your statement about bond contribution should be true for angle contribution as well but it is not the case. Angle contribution goes to zero and by “contribution” I mean the virial term divided by the volume. It is the same quantity that Lammps calculates when you use the compute pressure command with additional keys like “bond” or “angle”.
On the other hand, as we increase the volume of the box (which I run separate simulations for doing this, I don’t actually expand the simulation box in this case), we should retrieve the ideal gas behavior and Z factor must go to one. I plotted Log (Z) for several densities and it is clearly not tending toward 1 but instead, it is smaller than one. I don’t understand why the simulation of a simple gas with a simplified interaction potential (LJ with just the repulsive part) at very large volumes should result in a Z factor significantly smaller than one (~ 0.6).
I believe either I am doing something wrong or there is a bugs in the code.

What Steve says is correct. In fact, the total bonded virial of an
isolated must be non-zero in order to balance out the kinetic contribution
from interal degrees of freedom.

For N molecules in volume V, LAMMPS computes the pressure as:

PV = N*m*kT - 3V<dU/dV>

where m is the number of atoms in the molecule and dU/dV is the virial.
But the ideal gas law says:

PV = NkT

In the V->infinity limit, these two equation should agree, and so:

<dU/dV> = N*(m-1)*kT

So we expect the bonded virial to be non-zero. However, there may also be
problems with ergodically sampling dilute gases using MD, because energy
may not efficiently transfer between translational and vibrational modes.
You should pay close attention to the effect of the timestep, spring
constants, and thermostatting. I suggest Langevin.

Aidan

Dear Aidan,
When I rearrange your equation, I get:
<dU/dV> = N*(m-1)*kT.../3V
Why does the "3V" -> 1 when V-> infinity?

JRE

Yes, you are correct. I forgot a factor of V. As V->Infinity, the
extensive virial will remain finite, but the intensive quantity <dU/dV> ->
(m-1)kT*(N/V) -> 0. So the compressibility factor PV/NkT should approach a
value of 1 as V->Infinity. It may be the sampling ergodicity that is
causing the problem. You could try measuring the temperature of differnent
modes to see how well equipartition is satisfied.

Aidan

I forgot to mention that I have already tried a wide range of time steps from 0.01 fs to 5 fs and the problem was still there. I am using the “real” system units of Lammps in the simulations and based on TraPPE potential the angle force constant is around 60. So, as I said before, I have already tried a bond force constant range from 60 to 1000 in combination with different time steps and I got no significant change in the results. The Z factor is still around 0.5
As you suggested, I tried temp/langevin as well as temp/berendsen for thermostating instead of temp/rescale and this time the problem is getting worse. Z factor now is around 0.3 for different force constant and time steps. Also, I gave up on combining these commands with a “fix nve” command and used simply a “fix nvt” command to thermostat but again no success. I should mention that I used these commands with different tdamp values, namely 10, 100 and 1000.
On a separate try, in the calculation of Z factor I ignored the bond contribution of virial and I got Z factors around 1.04 for the same system which perfectly makes sense. Although I use very large simulation box, the system is not quite an ideal gas so we expect a number close to 1 and a little bit higher than that. With this way of Z calculation, changing the method of thermostating, time step and bond force constant doesn’t change the Z factor that much. It stays around 1.05 but it is obvious that the problem is not solved, I am just cheating here…
This experiment along with the fact that the angle contribution of virial pressure is fine and the angle term goes to zero as the volume increases shows that, I think, there is some thing wrong with the calculated bond contribution of virial pressure. If you think it is necessary, I can send you my input file and an abbreviated version of my data file. Also, I have an Excel file with a plot of log(Z) vs. density for propane at 300 K that I can send here if it is necessary.

I think you are jumping to conclusions in pinning the problem on the way
the bond virial is calculated. One piece of counter-evidence is that you
get a differ value for Z using different MD methods. If the virial was
wrong and the MD methods were all correct, you would get the same answer
for all MD methods. I suspect that not all vibrational modes are at the
same temperature. Equipartition is something that we can take for granted
in condensed phases, but in a gas, you need to verify it.

The fact that you get Z=1.05 by removing the bond virial is a good thing.
That eliminates a lot of possible error sources.

If you really want to test the bond virial, independent of ergodicity
issues, the way to do it is to build a dimer at zero temperature in a 3D
periodic cell with no nonbond interaction. Make a small deformations of
the cell and see if the LAMMPS pressure is consistent with deltaU/deltaV.
To be really safe, you should do this for the full stress tensor. This is
similar to how I compute elastic constants.

http://lammps.sandia.gov/doc/Section_howto.html#howto_18

Aidan

      Aidan P. Thompson
      01425 Advanced Device Technologies
      Sandia National Laboratories
      PO Box 5800, MS 1322 Phone: 505-844-9702
      Albuquerque, NM 87185 Fax : 505-845-7442
      E-mail:[email protected] Cell : 505-550-2614

Dear Aidan,
We decided to follow your suggestion but the simulation details seems a little bit ambiguous to me. I read the link that you sent and also I studied the input files in example/elastic directory of Lammps. In that example, the initial configuration is first relaxed by box/relax command, then it is minimized. After that, the box length is change followed by another minimization to calculate elastic force constants. Do you mean I should use the same input files in my case? ( of course with appropriate changes… )
I guess I need to run the simulation for a while at 0 K to find out if the P bond is calculated correctly or not. If I do that and force the system to have 0.0 total energy at 0 K, I fall into a glassy condition where nothing changes in an NVE ensemble simulation and the whole system freezes. On the other hand, if i relax the box, then run minimization followed by an NVE simulation, the very small bond energies in the system fluctuates between translational and internal degrees of freedom which is physically expected but because of the fluctuations (which has a constant amplitude and wavelength since I am using harmonic oscillator as the bond potential) one cannot calculate exactly the effect of V change on Ubond. So, I am confused.
I have tried to follow the strategy of examples/elastic input files but after second minimization command the energy of my system is exactly 0 and I don’t trust it. Although, I used that to calculate -3*dU(bond)/dV but I am getting a different number comparing to calculated P bond by Lammps. On a separate try, I came back to my propane system and this time with no thermostating method, I did NVE simulation and the problem was still there. So I didn’t even touch the T and let it go naturally but still the bond contribution was significantly far from 0 at very large volumes. I also tried morse potential instead of harmonic potential for bond contribution and again I got non-zero energies for bond term.

So, here is a list of things that I have done so far that might be use full to summarize here:
Propane at 300 K, NVT simulation with N=500 and V=10^6 cubic angs. Units=real system
I use Packmol software to create initial configuration and test that using VMD package.
Interaction potential: LJ/cut/opt with a cutoff distance equals to sigma. This turns off the attractive forces with a Barker-Henderson splitting.
time steps : 0.01, 0.1, 1.0 and 5.0
force constants: 20, 60 100, 200, 500, 1000 kcal/mol
thermostating method: temp/rescale, temp/Berendsen, temp/langevin and “fix nvt” option. I also tried “fix nve” with no thermostating.
Tdamp: 1, 10, 100, 1000 and 10000
minimization: minimize the initial configuration using minimize command: ON and OFF
bond potential: Harmonic and Morse
Turn on the attractive forces with a cutoff=10 A
Tail correction: ON and OFF
Now Lets assume that I am successful in doing the 0 K simulation that you suggested and it is perfectly consistent with itself which demonstrates that Lammps is fine and there is sth wrong with my input files. Considering all things that I have tried, what do you think I should do next? I just want to simulate a near ideal gas condition for a simple gas like ethane or propane.

So sorry for several questions and a long post and thank you in advance.

Hi Ahmadreza,

I guess I should have been a little more specific. The only thing you
should copy from the elastic constant script is the general procedure for
estimating strain derivatives using finite deformations. In your case, the
derivative you are after is dU/dV.

Because this is relatively simple, you don't actually need to write a
complicated script. Just compute the energy before and after a small
expansion, uisng say:

displace_box all x scale 1.01 y scale 1.01 z scale 1.01

Aidan

Ahmadreza,

I performed the test that I was talking about. The numerical pressure
converges to the one given by LAMMPS' harmonic bond virial. Once again, I
think you need to examine the issues of ergodicity, equipartition, etc.
Unlike Monte Carlo, MD relies on purely physical mechanisms to transfer
energy between different modes. This is not normally a problem, but
paradoxically, it grows in importance as you approach the ideal gas limit.
If you only have a few molecules undergoing infrequent collisions, then
you are not sampling an equilibrium distribution of molecular
conformations, unless you take specific measures. In principle, the
Langevin thermostat should take care of this, but you will have to verify
that it is doing what you need.

Aidan

[[email protected]... bondvirial] mpirun \-np 1 \.\./\.\./src/lmp\_mac\_mpi &lt; in\.virial | grep \#\#\# \#\#\# Delta Pressure Error \#\#\# 1\.0e\-8 1\.010999999 \-2\.864e\-06 \#\#\# 1\.0e\-7 1\.011 \-2\.863e\-06 \#\#\# 1\.0e\-6 1\.010993994 \-8\.869e\-06 \#\#\# 1\.0e\-5 1\.010911191 \-9\.1672e\-05 \#\#\# 1\.0e\-4 1\.009965954 \-0\.001036909 \#\#\# 1\.0e\-3 1\.000630199 \-0\.010372664 \#\#\# 1\.0e\-2 0\.9081782379 \-0\.1028246251 \#\#\# 1\.0e\-1 0\.0678184609 \-0\.9431844021 \#\#\# 1\.0 \-3\.374479863 \-4\.385482726 \[athomps@\.\.\.2929\.\.\. bondvirial\] more in.virial
# A LAMMPS bond

print "### Delta Pressure Error"

atom_style bond
atom_modify sort 0 0.0
# Soft potential push-off

read_data data.bond
mass * 1.0
bond_style harmonic
bond_coeff 1 50.0 0.75

thermo 1
thermo_modify norm no

run 0

variable etmp equal etotal
variable etotal0 equal \{etmp\} variable vtmp equal vol variable vol0 equal {vtmp}
variable ptmp equal press
variable p0 equal ${ptmp}

label loop
variable delta index 1.0e-8 1.0e-7 1.0e-6 1.0e-5 1.0e-4 1.0e-3 1.0e-2
1.0e-1 1.0
variable scale equal 1+${delta}

displace_box all x scale \{scale\} y scale {scale} z scale ${scale}
run 0

variable pnum equal -\{vol0\}\*\({etmp}-\{etotal0\}\)/\({vtmp}-\{vol0\}\) variable perr equal {pnum}-\{p0\} print &quot;\#\#\# {delta} \{pnum\} {perr}"
next delta
jump in.virial loop

[[email protected]... bondvirial]$ more data.bond
# A LAMMPS bond
           2 atoms
           1 bonds
           1 atom types
           1 bond types

  Atoms

      1 0 1 0.0 0.000 0.000 0 0 0
      2 0 1 0.5 0.4 0.3 0 0 0

Bonds

      1 1 1 2

Aidan,
OK. I could reproduce the pressure using the experiment that you suggested but with a slightly different procedure. First of all, I did it not only for one molecule but a simulation box filled with 500 molecules at 0 K and I got good results for both systems.
But, some thing still confuses me: you calculate dU/dV then multiply that by -V and recognize it as your P. How come your pressure has a unit of internal energy “U”? Besides, a couple of days ago, you mentioned that PV=…-3VdU/dV. So P=…-3dU/dV not P=…-3VdU/dV. What happened to the factor “3” in your calculation?
I could only reproduce P when I ignored the factor “3”, like you, and by changing the dimension of dU/dV to atm and I didn’t multiply by V at the end.
Here is my input file as you suggested:
units real
atom_style bond
atom_modify sort 0 0.0
bond_style harmonic
read_data 1ethane.data
compute Pke all pressure thermo_temp ke
compute Pbond all pressure thermo_temp bond
thermo_style custom temp press ke ebond vol c_Pke c_Pbond
run 0
displace_box all x scale 1.000001 y scale 1.000001 z scale 1.000001
run 0

and the data file for the case of one single molecule (I deleted the necessary “empty lines”):
LAMMPS Description
2 atoms
1 bonds
1 atom types
1 bond types
-10.0 10.0 xlo xhi
-10.0 10.0 ylo yhi
-10.0 10.0 zlo zhi
Bond Coeffs

Trappe Potential, MartinJPCB1998,

1 60.0 1.540
Masses
1 15.0350
Atoms
1 1 1 -0.424 -6.544 -2.897
2 1 1 0.883 -5.823 -2.521
Bonds
1 1 1 2

and the result:
Temp Press KinEng E_bond Volume Pke Pbond
0 0.36573017 0 2.8816027e-05 8000 0 0.36573017
Temp Press KinEng E_bond Volume Pke Pbond
0 0.36491708 0 2.8688158e-05 8000.024 0 0.36491708
If i calculate -dU/dV and change its unit [=] kcal/mol/A^3 to atm, I will get P=0.36532 which is good enough. Why don’t I need to multiply by “3” to reproduce the pressure?
I guess this experiment shows that my simulation has a problem, right? How can I find out the source of problem? I tried tmp/langevin but no good came out of that. Also, I tried respa instead off Verlet and I couldn’t improve the results.
This is getting a little bit frustrating. As you said, if the problem is related to equipartition and ergodicity, how can I fix that?

Your calculation looks correct. Pressure is defined as P = -dU/dV. My
script had an incorrect factor of V, which had no effect because V=1.0 in
my example. Also, you need to apply the correct unit conversion, which for
LJ units is 1. The factor of 3 that I mentioned previously was in a
different context.

So, we can conclude that the bond virial is correct. As I said from the
start, the underlying problem is probably a violation of equipartition. I
think that the Langevin thermostat should eliminate this issue, if it is
done right. But you should verify equipartition by computing the average
kinetic energy of each atom. Results will be sensitive to all the time
constants in the simulation: timestep, bond stiffness, Langevin drag
coefficient, mean collision time.

Dear Aidan,
In the Lammps’s manual on P. 240 I see that the pressure is computed using virial theorem through a force calculation step. This is consistent with what I find in Hensen&McDonald’s book (3rd edition) on P.19 for the pressure. Also I looked at the source code of Lammps (compute_pressure.cpp) and it seems to me that it follows what is mentioned in the manual. On the other hand, you said Lammps calculate P by:
P= … -dU/dV
After all, how did you bring this equation to the table? If this is the case, how does Lammps keep the entropy constant? Because this relation is only valid when the entropy is constant. I can understand how the entropy is constant for our test of 0 K simulation but what about other situations?

Ahmadreza,

You are asking me to teach you a few weeks of a graduate stat mech course
in an e-mail. I can't do it. But maybe I can sketch out an answer. P =
-dU/dV is only valid at T=0; U is the potential energy function. No
entropy in play. dU/dV is just another way of writing what you call the
force virial, modulo a minus sign and a factor of V. At non-zero T, a
kinetic energy term must be added in. Look this stuff up in any good book
on atomistic simulation. I recommend Allen and Tildesley or Frenkel and
Smit.

Aidan