Dear all,
My system has N particles, in which N1 particles(group fixed) are fixed with “fix setforce” command and the other N-N1 particles are free(group free).
Now I am trying to do run dynamics in NPT ensemble. In “fix NPT” the temperature is compute by: “compute ID all temp”, so fixed N1 atoms also are counted for computing temperature, but in reality, temperature should be compute by “compute ID free temp”.
Does anyone know how to solve this problem? Thanks!
Wei
Dear all,
My system has N particles, in which N1 particles(group fixed) are fixed with
"fix setforce" command and the other N-N1 particles are free(group free).Now I am trying to do run dynamics in NPT ensemble. In "fix NPT" the
temperature is compute by: "compute ID all temp", so fixed N1 atoms also are
counted for computing temperature, but in reality, temperature should be
compute by "compute ID free temp".Does anyone know how to solve this problem? Thanks!
this is well documented issue.
both the thermostatting howto
http://lammps.sandia.gov/doc/Section_howto.html#howto_16
and the fix npt documentation mention how you
http://lammps.sandia.gov/doc/fix_nh.html
can use fix_modify
http://lammps.sandia.gov/doc/fix_modify.html
to replace the default temperature compute that
is created by the fix with a custom one.
so please re-read the documentation.
axel.
Thanks! axel.
Following the above question.
In the fix NPT, I only want to control the pressure from part of the atoms. But the “compute pressure” can only compute the pressure for all atoms. I am not able to modify the pressure computation. One possible way I am thinking is to use “compute stress/atom” along with “compute reduce sum”, but it does not generate a pressure scalar. Do you have any idea about this? Thanks!
Wei
Thanks! axel.
Following the above question.
In the fix NPT, I only want to control the pressure from part of the atoms.
But the "compute pressure" can only compute the pressure for all atoms. I am
...and for a good reason.
not able to modify the pressure computation. One possible way I am thinking
is to use "compute stress/atom" along with "compute reduce sum", but it does
not generate a pressure scalar. Do you have any idea about this? Thanks!
yes, i have: you haven't read the documentation well enough.
but even if you'd do the computation, it would result in a bogus
result for the pressure. even if atoms don't move, they still
contribute to the pressure. the pressure of a subgroup of atoms
is an ill-defined property. this have been discussed on this
mailing list only a few hundred times. please have a look at
the mailing list archives and think about how you define
and measure pressure macroscopically.
axel.
Thanks for your reply. I do not totally agree with you.
According to virial theorem:
p=NkT/V + sum_on_i{<r_i * f_i>/3V}, where summation is on atoms in subgroup, and f_i is force between subgroup atoms. This is a well defined quantity. Of course, macroscopically, the pressure comes from the interaction between subgroup atoms and others.
Wei
Thanks for your reply. I do not totally agree with you.
According to virial theorem:
p=NkT/V + sum_on_i{<r_i * f_i>/3V}, where summation is on atoms in subgroup,
and f_i is force between subgroup atoms. This is a well defined quantity. Of
no. you have a well defined formula, yes, but you also
have one big problem there: what is the volume?
how do you define the volume of a single atom.
the only well define volume that you can have in
an MD is the volume of the entire box.
as i was saying, this topic has been discussed
many times. please have the courtesy and look
up the discussions in the mailing list archives as
i recommended to you. we don't have to go over
the exact same arguments again.
axel.
1) There is no need to use setforce. You can apply "fix npt" or "fix
nvt" to the "free" (mobile) group of atoms instead of group "all".
2) At constant pressure I recommend using "neigh_modify exclude group
fixed fixed" improves (but does not completely fix) the virial. (But
doing that and using "fix rigid" seems to work well.)
3) In a reply to a recent post last week:
http://sourceforge.net/mailarchive/message.php?msg_id=30018547
I posted an input script which immobilized a large solute:
http://www.moltemplate.org/examples/translocation/run.in.npt
In particular, the relevant portion of that script is below:
----------- excerpt: -----------
fix 1 fixed rigid single force * off off off torque * off off off
#(Neither Trung or Steve Plimpton use fix rigid for immobilizing objects, but
# I noticed that at NPT, it does a significantly better job of maintaining
# the correct volume)
neigh_modify exclude group fixed fixed
compute tempFree free temp
compute pressFree all pressure tempFree
thermo_style custom step c_tempFree c_pressFree temp press vol
# Set temp=300K, pressure=500bar, anisotropic equilibration
fix 2 free npt temp 300 300 100 aniso 500 500 1000.0 dilate free
fix_modify 2 temp tempFree
-------- end of excerpt -------
(Small changes: in the excerpt above I renamed the groups to "fixed"
and "free", to match the immobile and mobile groups in your example.
Credit goes to Trung who instructed me to use "fix_modify" to use the
temperature computed from the mobile group of atoms. I hope I am
reproducing it correctly here.)
Cheers
Andrew
P.S.
---- Regarding the virial: ----
If you want to keep the coordinates of the rigid atoms from being
expanded or compacted as the volume is increased or decreased, you
should use "dilate free", as I did above. In that case, only the
positions of the free atoms will be rescaled.
Under the hood, this requires a modification to the virial: all of the
forces acting on the immobilized "fixed" atoms must be excluded from
the virial. (Or equivalently, the constraint forces used to keep the
object motionless should be included in the virial. Perhaps this is
more intuitive way to think about the issue: All the forces acting on
all the particles must be included in the virial, including the
constraint forces. I attached a quick proof I wrote to motivate this
to myself a year ago but it's not very well written. In any case,
"neigh_modify exclude group fixed fixed" only eliminates forces
between atoms within the "fixed" group. This helps but it does not
eliminate contributions to the virial from the remaining mobile "free"
atoms acting on the fixed atoms.)
However the above combination of "neigh_modify exclude group fixed
fixed" AND "fix rigid" seems to work.
The post at:
http://sourceforge.net/mailarchive/message.php?msg_id=30018547
includes a discussion of the pros and cons of this approach.
Cheers again
Andrew
virial_immobile_atoms_sm.pdf (101 KB)