temperature dof vs pressure dof

Dear all,

I am currently messing with the compute_temp_* and compute_press codes
and, inspired by the thread started by Giacomo Saielli
http://lammps.sandia.gov/threads/msg35483.html , I think there's
something not yet clear about the degrees of freedom used for
temperature and pressure computations.

When computing its global scalar, the ComputePress class solves the
pressure as
temp->dof / 3.0 * force->boltz * temp->scalar / volume (+virial terms)
I think this was intended to correctly take into account the
constrained
translational degrees of freedom due to fix_rigid or alike.
But when the temperature takes into account extra degrees of freedom
such as rotations (temp/(a)sphere), the formula above will not yield the
_real_ pressure. For an indeal gas PV=NKT independently of the
rotational
and vibrational degrees of freedom of its molecules!

I think the ComputePess class must be provided a
ComputePress::compute_dof() method to compute its own translational dof
and it must be possible to modify its extra constrained dof by means of
the fix_modify extra command applied directly to the pressure compute
instead of the termperature compute. Then the pressure will be the same
whatever dof are used by whatever temperature compute if the computed
temperature is the same.

It could simply be translational_dof = 3 * natoms - extra_dof - fix_dof,
just like the simple ComputeTemp does. I would do it myself, but I
wouldn't know how to handle rigid bodies.

What do you think?

I’ll let Aidan comment on this.
I thought that the kinetic contribution to pressure
was consistent since both the temperature
itself and the DOF were taken from the temperature
class. Note that you assign a different temp compute
to compute pressure if you like, e.g. one that only
calculates translational temp.

Steve

Hi D.

Even though both temperature and pressure are intensive properties, their statistical mechanics definitions involve different normalizations. Temperature is normalized by the number of degrees of freedom, whereas pressure is normalized by the volume. Hence the quantity Ndof is not required when computing pressure. It appears in the formula you cite below only for the purpose of unnormalizing the temperature, which recovers the desired extensive quantity: kinetic energy. In the presence of constraints such as fix rigid, the constrained degrees of freedom do not contribute to the kinetic energy, and the correct ideal gas law is observed. That’s the theory. In practice, it seems to be working right in LAMMPS for the standard molecular constraint methods such as SHAKE. If you have a counter-example that you can demonstrate where:

PV != NkT, as N/V —> 0

let us know and I’m sure we can “fix” it.

Aidan

Whenever the dof accounted for temperature are not the only
translational dof of the compute pressure group the pressure will not be
correct. In particular
- the default behaviour of npt/asphere will lead to PV=2NKT
- if a subgroup of M atoms is used to estimate the temperature of the
all N group then PV=MKT
- when atoms are bonded into molecules by any means so some dof are
constrained but every molecule keeps 5 or 6 dof PV=5NKT/3 or 2NKT

That is: only
compute myT all temp
compute myP all pressure myP
will lead to real pressure but only if the atoms/molecules do not have
not translational degrees of freedom the temperature accounts for.

El dv 12 de 04 de 2013 a les 15:09 +0000, en Thompson, Aidan va escriure:

For an ideal gas of flexible molecules, the internal
kinetic energy contribution is balanced (in the ensemble average) by the
internal force virial contribution. Alternatively, the internal kinetic
energy contribution and force virial contribution can also be omitted,
yielding a numerically distinct definition of the pressure that yields the
same ensemble average.

The point is, there is more than one way to compute
a correct pressure, as long as it is done consistently. The same is true
for constrained systems.

In the examples you cite below, it would be helpful if you could
demonstrate your assertion that the pressure is wrong by providing output
from a small test. Also, better not use npt, but rather nvt. NPT with
rigid constraints introduces additional complications.

Aidan

For the asphere problem I point you to the thread started by Saielli http://lammps.sandia.gov/threads/msg35483.html , where Michael Brown pointed the problem when he said:
"npt/asphere creates a temp/asphere compute and uses this for pressure computation. Your thermodynamic output is probably using a pressure compute that does not include angular terms in the temperature."

I'll skip the subgroup issue since it is a quite clear example and can be solved by chosing the all group to compute temperature, although I think the feature to sample temperature from subgroups should be allowed (and it will when the ComputePress class has its own compute_dof() method).

For the bonded atoms issue, I have recompiled lammps with the RIGID package and extended the example in examples/rigid/in.rigid in the lammps distribution as follows:

### INPUT ###

# Simple rigid body system

units lj
atom_style atomic

pair_style lj/cut 2.5

read_data data.rigid

velocity all create 100.0 4928459

# unconnected bodies

group clump1 id <> 1 9
group clump2 id <> 10 18
group clump3 id <> 19 27
group clump4 id <> 28 36
group clump5 id <> 37 45
group clump6 id <> 46 54
group clump7 id <> 55 63
group clump8 id <> 64 72
group clump9 id <> 73 81

fix 1 all rigid group 9 clump1 clump2 clump3 clump4 clump5 clump6 clump7 clump8 clump9

neigh_modify exclude group clump1 clump1
neigh_modify exclude group clump2 clump2
neigh_modify exclude group clump3 clump3
neigh_modify exclude group clump4 clump4
neigh_modify exclude group clump5 clump5
neigh_modify exclude group clump6 clump6
neigh_modify exclude group clump7 clump7
neigh_modify exclude group clump8 clump8
neigh_modify exclude group clump9 clump9

compute virial all pressure thermo_temp virial
compute lammps_press all pressure thermo_temp ke
# ( 9 rigid bodies * 3 translational dof - 3 center-of-mass dof ) / 3 dimensions = 8
variable real_press equal 8*c_thermo_temp/vol
variable press_dof equal c_lammps_press*vol/c_thermo_temp

timestep 0.0001
thermo 50
thermo_style custom step temp epair emol etotal press c_virial c_lammps_press v_real_press v_press_dof
run 10000

### OUTPUT ###

LAMMPS (18 Feb 2013)
Reading data file ...
   orthogonal box = (-12 -12 -12) to (12 12 12)
   1 by 1 by 1 MPI processor grid
   81 atoms
9 atoms in group clump1
9 atoms in group clump2
9 atoms in group clump3
9 atoms in group clump4
9 atoms in group clump5
9 atoms in group clump6
9 atoms in group clump7
9 atoms in group clump8
9 atoms in group clump9
9 rigid bodies with 81 atoms
Setting up run ...
Memory usage per processor = 1.86701 Mbytes
Step Temp E_pair E_mol TotEng Press virial lammps_p real_pre press_do
        0 115.29439 5235.9179 0 5272.2142 -2.7403788 -2.8821616 0.14178275 0.066721293 17
...
      100 16298.442 136.66184 0 5267.653 16.444229 -3.5987048 20.042934 9.431969 17
...
     1000 16738.586 -0.051135144 0 5269.5036 14.496229 -6.0879693 20.584198 9.6866814 17
...
    10000 16738.49 0 0 5269.5246 11.395075 -9.1890059 20.584081 9.6866262 17
Loop time of 0.241194 on 1 procs for 10000 steps with 81 atoms

Whenever the dof accounted for temperature are not the only
translational dof of the compute pressure group the pressure will not be
correct. In particular

sorry - I still don’t understand your point. The compute pressure
group is always “all”, there is no other choice.

you have the freedom to select what compute temperature is used
with compute pressure (4th argument). you have the freedom
to choose how many dof are used with the compute temperature
you select (compute modify extra). So if you want to
define a pressure that uses only translational DOF with the
pressure it computes, you can specify that.

So are you saying you don’t like what LAMMPS is computing
as the default for aspherical particles? Or that you cannot adjust
things to compute what you want?

Steve

For the asphere problem I point you to the thread started by Saielli http://lammps.sandia.gov/threads/msg35483.html , where Michael Brown pointed the problem when he said:
“npt/asphere creates a temp/asphere compute and uses this for pressure computation. Your thermodynamic output is probably using a pressure compute that does not include angular terms in the temperature.”

I’ll skip the subgroup issue since it is a quite clear example and can be solved by chosing the all group to compute temperature, although I think the feature to sample temperature from subgroups should be allowed (and it will when the ComputePress class has its own compute_dof() method).

I went back to the thread you mention and your previous posts. The thread is advocating computing
the pressure as:

compute  rot all temp/asphere
group      spheroid type 1
variable     dof equal count(spheroid)*3+3 # For BIAXIAL ellipsoids, for uniaxial, 2 degrees of freedom per ellipsoid
compute_modify rot extra ${dof}
compute my_pressure all pressure rot

Are you saying that is the same as

compute  rot all temp
compute my_pressure all pressure rot

i.e. you want to compute the pressure using only
the translational dof of the aspherical particles?

If so, then why not just use those 2 commands?
Since you get to choose what temperature (and its dof()) is
used with compute pressure, I don’t see why it is better
to have a single compute pressure dof() (unchoosable), instead of letting
the user select which temperature dof() to use,
if as Aidan says, there is some ambiguity, as to what is the
right choice.

Steve

Concerning the asphere issue:

compute rot all temp/asphere
compute_modify rot extra ${dof}
compute my_pressure all pressure rot

This would yield the correct pressure but would also turn the temperature compute wrong.

compute rot all temp
compute my_pressure all pressure rot

This will yield the correct pressure and temperature. But what is the point for having a wrong default behaviour for the np*/asphere ? And I don't feel as if this is warned enough in the doc so any end-user is aware of this issue.

I think it would be clearer if the ComputePress class had its own dof count and could be modified, say

compute myT all temp/asphere
compute myP all pressure myT
compute_modify myP extra ${dof} # this is not yet possible

This would allow to forget about the dof used for temperature and even to indirectly sample temperature on a different group of atoms.

But don't miss the rigid bodies simulation I attached before! The only solution there would be

compute fakeT all temp
compute_modify fakeT extra ${dof}
compute realP all pressure fakeT
compute realT all temp

And it would only solve the scalar pressure problem *but* the pressure tensor would still be wrong! And, again, I am afraid nobody has this in mind when conducting simulations with lammps.

Obviously all this can be solved by successive temperature and pressure computes accounting for different dof to compute either temperature, either pressure, either scalar, either tensor. But I don't think this is a good solution.

If the ComputePress objects had their own dof count:
1) they would yield the correct pressure independently of the temperature compute dof and group
2) the asphere package would be right by default
3) only extra pressure dof would have to be manually constrained for atoms bonded into molecules, unless the ComputePress::compute_dof() method is really smart

Al 15/04/13 15:36, En Steve Plimpton ha escrit:

Hi D.

The fix rigid example gets to the heart of this issue very nicely. Thanks
for sharing it. After examining this output for a while and looking at the
code, I conclude that the LAMMPS pressure calculation is consistent with
the ideal gas law behavior i.e.

<Press> = N*k*T/V

where N is the number of rigid molecules. The reason that this works is
because LAMMPS computes the pressure based on the total kinetic energy of
all atoms, with out regard to constraints, but also includes a force
virial due to the constraint forces. The sum of kinetic and force virials
gives the correct pressure. This approach eliminates the need to
explicitly worry about the correct number of degrees of freedom. I have
not looked at the code for fix asphere, but if it is not handling pressure
the same way as fix rigid, it probably should.

Aidan

More detailed analysis

Hi Daniel,

Concerning the asphere issue:

compute rot all temp/asphere
compute_modify rot extra ${dof}
compute my_pressure all pressure rot

This would yield the correct pressure

This is wrong. The pressure computed by LAMMPS is unaffected by Ndof. Try
it. The instantaneous pressure depends on 3 things:

1. The volume
2. The kinetic energy = Sum[mi*vi*vi/2]
3. The force virial = Sum[Fi*Ri]

I don't see how Ndof can be included in this list.

Aidan

Actually, this was the solution proposed by Steve and Michael Brown in that thread. Note that the variable ${dof} only accounts for the non translational dof thus when removed from the compute, only the 3N translational dof stay.

Al 15/04/13 21:48, En Thompson, Aidan ha escrit:

Not at all. I am sorry I chosed poorly clear names for my variables. I'll try again:

Note that the actual pressure will be the sum of the pressure due to a molecular density at a given pressure (ideal gas pressure) plus the pressure due to intermolecular interactions (virial pressure)

In particular, in the simulation I attached, the neighbour settings are modified so that intramolecular interactions are not even computed at any timestep. Therefore the virial does not know anything about the bonding constrains.

What I am complaining is that the ideal pressure contribution is twice as much what it should be.

compute virial all pressure thermo_temp virial
compute wrongIdealP all pressure thermo_temp ke
# thermo_press = wrongIdealP + virial
variable wrongN equal c_wrongIdealP*vol/c_thermo_temp+1
variable rightIdealP equal 8*c_thermo_temp/vol
variable rightTotalP equal v_rightIdealP+c_virial

The fact that in this simulation rightIdealP is similar to wrongTotalP is merely casual and proves nothing, since the actual pressure will be rightIdealP+virial.

I have not completely understood your dof count, but note that the wrongN variable now is exactly the number of molecules the pressure compute thinks there are. Note it is exactly 18 while there are only 9 rigid bodies in that simulation.

Al 15/04/13 20:47, En Thompson, Aidan ha escrit:

Daniel,

Let me repeat two important things that you may not have picked up on the
first time:

1. With your fix rigid example, there *is* a contribution from the force
virial, even though intramolecular interactions are not computed. Read my
explanation below, or just look at the output from your simulation.

2. With fix asphere, the number of degrees of freedom specified for the
temperature compute using the compute_modify extra command has *no effect*
on the computed pressure. This is because the pressure calculation uses
kinetic energy, not kinetic energy divided by Ndof.

Aidan

Okay, I ran a longer simulation and saw what you said. I'll explain it
here so neither people reading this thread nor I get confused.

fix rigid computes a correction to pressure that can be accessed through
compute rigidP all pressure thermo_temp fix

Therefore, the ideal pressure due to this rigid bodies at a given
density and temperature is
compute idealP all pressure thermo_temp ke fix
instead of
compute wrongIdealP all pressure thermo_temp ke

So if the pressure is computed as a variable then the virial
contribution must not account the fix rigid correction:
compute virial all pressure thermo_temp pair
compute wrongVirial all pressure thermo_temp virial
variable idealP equal ($N-1)*temp/vol # say N is the number of bodies
variable totalP equal v_idealP+c_virial

And, as you said, despite the two ways are instantaneously different,
they will average the same (attached plot) and the question moves to
whether each fix or whatever mechanism to bond particles is computing
its correction to pressure, which can be checked through
compute fixP all pressure thermo_temp fix
compute bondP all pressure thermo_temp bond
And so forth.

So, Aidan and Steve, thanks a lot for your patience.

pressure.png

Hi Daniel,

That's good news. The plots looks better than I expected. I agree with
your explanation below. Thanks for persevering.

Aidan