Dear Prof. Kohlmeyer
Thanks very much for your reply.
There are still some questions that confused me. You said that “force of
bond and angle do not contribute to the forces on the rigid bodies as a
whole and if they are not exactly zero they will contribute to the virial”
Is this means that the force of bond and angle in rigid molecules have
influence on the fix rigid/nvt?
no. only your computed pressure will be off. when you run with
constant volume, it is only a diagnostic property and will not impact
the trajectory.
1. In the manual, it says that “For computational efficiency, you may wish
to turn off pairwise and bond interactions within each rigid body, as they
no longer contribute to the motion. The neigh_modify exclude and
delete_bonds commands are used to do this” I don’t know whether the
computation of bond and angle force in rigid molecule only affect the
computational efficiency, or it will influent the temperature control of the
system?
if you look at the breakdown of the computational cost, you will see,
that the cost of the bonded forces in your specific system is
negligible.
2. Whether can I eliminate the bonds and angles in the data file directly
for rigid molecule?
you can do this, but in your case, it would be much simpler to only
set the force constants to zero.
3. After turn off bond interaction, I use the “neigh_modify exclude molecule
rigid” to turn off the pairwise interaction of atoms in the same molecule,
but a warning occurs because “Excluding pairwise interactions will not work
correctly when also using a long-range solver via the kspace_style command”
I want to know how to deal with this warning and whether the computation of
pairwise interaction in rigid molecules affect the simulation.
this is a warning and not an error, because the code cannot know,
whether a neighbor list exclusion is proper or not. when you define
bonds or angles or dihedrals in your topology, the corresponding
interaction will be excluded from neighbor list for pairwise
interactions. since you have no dihedrals defined, there is going to
be a small difference, as your exclusion of all intra-molecular pairs
will also exclude the interactions between the two hydrogens. however,
if you remove/delete all bonds and angles, you *have* to use that
neigh_modify command, otherwise your intra-molecular contribution of
the virial will be quite from the pair-wise interactions from bonds
and angles.
in short, you either have to do what you did with deleting/removing
bonded interactions and "neigh_modify exclude molecule", or you keep
the bonds and angles and just set their force constants to zero.
computational efficiency is not affected when using kspace. if you
care about performance, you should use one of the rigid/small fixes.
4. My simulation follows a published paper, in which the npt ensemble is
used to simulate a single particle. I hope I have a comparison of the
an NPT ensemble is only possible for a bulk system. so this has to be
either a typo or the authors of the paper (and its reviewers) are
apparently lacking proper knowledge in statistical mechanics.
unfortunately, these days the fact that something got published does
not automatically mean that it has seen a thorough review and can be
trusted in every detail.
simulation, so that my simulation is convincing. If I use nvt ensemble, I
don’t know whether it still is convincing?
you won't have an NVT ensemble with your simulation either. that one
also requires a bulk system. the purpose of the nose-hoover thermostat
is not only to maintain temperature at the desired value to mimic the
embedding of the simulated system into a bulk of the same material,
but also to induce fluctuations that are due to the bulk system. the
latter affects the partition function and thus makes it that you can
get the same statistical sampling as in a bulk system, which in turn
makes it an NVT ensemble.
what you are doing is to run an isolated cluster/droplet with a
nose-hoover thermostat. even using the thermostat beyond equilibration
is questionable, since for an isolated object, where is the heat
reservoir (i.e. the embedding bulk) that it couples to?
as i outlined in the previous e-mail, in my personal opinion, the most
correct way to compute the interactions for an isolated
droplet/cluster is:
- use shrinkwrap non-periodic boundaries
- do not use kspace, but lj/coul/cut with a large enough cutoff. that
will make it so both lj and coulomb computed *exactly*.
- use either "neigh_modify exclude molecule" with bonds/angles removed
or set all force constants for bonds/angles to zero
- use "fix momentum 100 linear 1 1 1 angular rescale" or something
similar to avoid the flying/rotating ice cube syndrome.
using periodic boundaries and kspace/pppm is usually acceptable, too.
provided the distance between the periodic images is large enough. but
since LAMMPS doesn't have a poisson solver to decouple the coulomb
interactions between the periodic images, there is a risk of (small)
artifacts, depending on the magnitude of the multipole interactions
that forms between the droplet and its periodic images. if the
droplet/cluster is much larger, then using the full cutoff will become
computationally too demanding and kspace is the only viable choice.
axel.
