Immobilized atoms with PPPM in NPT

Hello,

I am trying to simulate a box of SPC/E with an immobilized hydrophobic solute in the NPT, using PPPM for long-range electrostatics. What is the best way to immobilize the hydrophobic atoms? Is it sufficient to only integrate the mobile atoms using “fix 1 mobile_group npt …”?

I have read in some locations that for the NPT ensemble, a “neigh_modify" command is required to turn off force calculations between immobilized atoms. However, the neigh_modify page of the LAMMPS manual says that neigh_modify and PPPM together will lead to incorrect dynamics.

Thanks,

Michelle

Hello,

I am trying to simulate a box of SPC/E with an immobilized hydrophobic
solute in the NPT, using PPPM for long-range electrostatics. What is the
best way to immobilize the hydrophobic atoms? Is it sufficient to only
integrate the mobile atoms using “fix 1 mobile_group npt …”?

​the easiest way to immobilize atoms in your case would be to set their
velocity​ to 0.0 0.0 0.0 and use fix setforce 0.0 0.0 0.0 on them. however,
there is still an issue with with fact that fix npt with expand and shrink
the box and that will include the immobilized particles. if this is not a
problem (e.g. if you only have a single layer of immobilized atoms and you
couple fix npt only to the direction orthogonal to it), then you don't need
to worry. if this is a problem, then you can use fix npt only on the
mobile group and use the dilate partial option, but in this case you have
to choose the origin around which the box is dilated very carefully.
overall, it may be much easier to use a fixed volume with a non-periodic
direction and use a soft or reflective walls to contain the mobile atoms.

I have read in some locations that for the NPT ensemble, a “neigh_modify"
command is *required* to turn off force calculations between immobilized
atoms.

​where?​

However, the neigh_modify page of the LAMMPS manual says that neigh_modify
and PPPM together will lead to incorrect dynamics.

​neigh_modify exclude will remove pairs of atoms from the pairwise force
computation, but that also removes their contribution to the stress tensor
and thus the pressure. for fix npt, the pressure is always computed for the
entire system, so this can be beneficial, if the forces between the
immobile particles are unphysical, but if those particles have a partial
charge, they will still be included in the PPPM calculation, and then they
will still cause a residual force due to the pppm calculation, and those
forces will be dependent on the coulomb cutoff (and they won't without
neigh_modify exclude).

again, probably the best approach for your setup is to try and avoid fix
npt altogether.

​axel.​

You said the molecule you are immobilizing is hydrophobic. If you can
approximate it as completely charge-free (all atoms have charge 0),
then I suspect that the WARNING message about long-range
electrostatics and fix rigid can be ignored. However I doubt this
applies to your molecule.

I have read in some locations that for the NPT ensemble, a “neigh_modify"
command is *required* to turn off force calculations between immobilized
atoms.

That might have been me who said that.

First of all, if you are running NVT simulations, then it does not
matter. (In that case, we all agree that there's no need to use
"neigh_modify exclude" or "fix rigid".) As Axel said, running NVT
simulations with

If you are running NPT simulations, then it is not clear.

Several people have tried using "fix setforce immobile_group 0 0 0"
without success. (Me too.) Using neigh_modify exclude seemed to be
necessary in my case. In my case I had enormous repulsive Lennard
Jones interactions between immobile particles which were dominating
the virial. Using "neigh_modify exclude" was one way to turn those
forces off which helped a lot. (I also found "fix rigid" useful.)

I recommend you try it both ways, with and without "neigh_modify
exclude", and also with and without "fix rigid". (Let us know what
works best! If you actually get "fix setforce" to work, let us know
as well.)

-------------- history -------------

There was some controversy regarding how to run simulations of
immobilized objects in LAMMPS in NPT conditions. In spite of Steve
and Trung patiently listening and trying to convince me otherwise, I
was advocating that people use both "neigh_modify exclude" and "fix
rigid" on the group of atoms you wish to immobilize.

Under NPT conditions, Trung and I ran some tests in 2011. We both
agreed that "neigh_modify exclude" was necessary. However, when I
tried it, my results suggested that "fix rigid" on the immobile atoms
seemed to do a slightly better job of preserving the volume of the
simulation, before and after immobilization.

These test results are shown here in these 4 plots:

http://sourceforge.net/p/lammps/mailman/attachment/CAEbHWTTLawUHd0h0n95Wc-cncguNWCvvH2%2BAGRjdQmLhOW6vMg%40mail.gmail.com/1/

http://sourceforge.net/p/lammps/mailman/attachment/CAEbHWTTLawUHd0h0n95Wc-cncguNWCvvH2%2BAGRjdQmLhOW6vMg%40mail.gmail.com/2/

(Trung's approach which does not use "fix rigid" is the blue curve.
My approach which is identical except that it uses "fix rigid" is the
green curve. Both approaches use "neigh_modify exclude".)

However the "fix rigid" command has changed since 2011. So everything
I am saying could be irrelevant. That's why I suggest trying it
yourself and see which way works best.

Here is a long list of posts on this topic:

http://sourceforge.net/p/lammps/mailman/message/30018547/
http://sourceforge.net/p/lammps/mailman/message/28151383/
http://lammps.sandia.gov/threads/pdfzpIPCFTvWW.pdf
http://lammps.sandia.gov/threads/msg42104.html
http://lammps.sandia.gov/threads/msg23401.html
http://lammps.sandia.gov/threads/msg42125.html
http://lammps.sandia.gov/threads/msg23341.html
http://lammps.sandia.gov/threads/msg29868.html

Cheers
Andrew

If you simply use fix npt with a group that includes the mobile atoms,
then only those will be time integrated. All the atoms will be remapped
as the box dilates/contracts, unless you use the dilate keyword.

None of this has anything to do with PPPM.

Steve

Hello,

Thanks for the tips. Andrew - yes, I got the template for NPT settings from moltemplate’s examples, so the suggestions must have been yours. I am currently approximating the solute as charge-free, so I’d guess that the PPPM/neigh_modify warning is okay to ignore. I am running everything without the “fix rigid,” and I had run some simulations with “fix rigid” before and didn’t find that it made a difference.

Since my particles have no charge, I also tried running without the neigh_modify and found that the volume conservation still looks fine.

I have not tried anything with “fix setforce" yet, and I think that with these charge-free particles there is no need. However, I would like to run with charged particles later on, so that is something I’ll explore.

Thanks,
Michelle