I have a hybrid system that includes both rigid (water molecules) and non-rigid body styles. I was trying to perform NPT dynamics. As the LAMMPS document suggested, there are three options to do this:
Use one of the 4 NPT or NPH styles for the rigid bodies. Use the dilate all option so that it will dilate the positions of the non-rigid particles as well. Use fix nvt (or any other thermostat) for the non-rigid particles.
Use fix npt for the group of non-rigid particles. Use the dilate all option so that it will dilate the center-of-mass positions of the rigid bodies as well. Use one of the 4 NVE or 2 NVT rigid styles for the rigid bodies.
Use fix press/berendsen to compute the pressure and change the box dimensions. Use one of the 4 NVE or 2 NVT rigid styles for the rigid bodies. Use fix nvt (or any other thermostat) for the non-rigid particles.
I followed these suggestions. Here is my command:
fix integrationWater water rigid/nve/small molecule
fix integrationNPT mobile npt temp {temp} {temp} {tdamp} x {pressure} {pressure} {Pdamp} y {pressure} {pressure} {Pdamp} z {pressure} {pressure} {Pdamp} drag 0.1 dilate all
After the simulation is done, I checked the coordinates of rigid water molecules. Found that the internal structure of rigid water molecule has been changed, the atom distance between O and H is not fixed at 0.9572 ( TIP4P/ICE),. I have tried all the suggested opptions. None of them work.
Thanks for your reply. In my system, some molecules are non-rigid, some are rigid. I didn’t manage to use fix shake for the four-site water clusters (the massless charge points M are explicitly added, use pair_style lj/cut/coul/long correspondingly). The simulation is just suspended. Then I try to use fix rigid for water molecules.
I think I have fixed this error. When the command line fix npt for non-rigid bodies used before fix rigid/nve/small for rigid bodies, the problem is solved.
You need to use pair style lj/cut/tip4p/long instead to realize a TIP4P style water model. This is much more efficient than using fix rigid as you can have a much larger timestep and have fewer particles in the neighbor list and thus faster computation.
Since my model includes other molecules apart from TIP4P water molecules, I have tried to use lj/cut/tip4p/long for all pair interactions, but this is more expensive actually.
It cannot be as expensive as having to use a much shorter timestep. For fix rigid and water molecules you may be required to use a 5-10x shorter timestep than for use with fix shake.
Thanks Alex. I’ll do some more tests on this. But using timestep of 1 fs and 2fs works well. In my understanding, when I use lj/cut/tip4p/long for all pair interations , this will make other pair atoms search massless charge point as well, similar to tip4p water. This is probably more expensive.
I doubt that you will get good energy conservation. Integration of the rotational degrees of freedom is much more sensitive because the moments of inertia differ in different directions.
That is not how the TIP4P pair styles work. The position of the M site is computed on-the-fly and that information is cached and it is only accessed for coulomb interactions and only if one of the participating atoms is a water hydrogen or oxygen. For everything else, the pair style operates exactly as lj/cut/coul/long. When using an explicit site, the point M will be included in all neighbor lists in addition to the oxygen site and thus distances are always computed even though the oxygen does not carry a charge and the M site has a LJ interaction with epsilon = 0.