Fluctuating bond lengths during rigid/npt fix


I am simulating rigid CO2 molecules equilibrated in the NVE ensemble followed by NVT and NPT ensembles. I am employing the rigid/nve, rigid/nvt and rigid/npt fixes for this purpose.

I output the average bond lengths between carbon and oxygen in my log file.
I find that the average bond length equals the imposed value in the NVE and NVT ensemble,
however there are small deviations in the NPT ensemble. These deviations appear in the fourth or fifth digit after the decimal. I would have expected that the bond length equals the imposed value irrespective of the ensemble.

I have uploaded the input files and a sample log file.

Any leads on what might be on here ?


CO2.txt (449 Bytes)
simProg.log (49.8 KB)
simulation.in (4.2 KB)

Hi @thejas_hc,

As far as the documentation tells, all atoms are rescaled during NPT simulations even with the rigid version of the integrator. This means that distances between atoms are changed and so are angles (note that you see no change in angles because you have linear molecules). There are barostat algorithms in the literature that rescale molecule positions only by computing the center of mass displacements and moving the whole molecules, but as far as I know none is implemented in LAMMPS at the moment. I also do not know how those could be correctly combined with rigid body dynamics.

Note that making your simulation run longer and looking at the time average of the bond length (your compute 3), it is very close to 1.60 A as you initially set.

Hi @Germain,

many thanks for your reply.

Is there a reason why the rescaling of the bonds and angles are applied only in the NPT ensemble, and not the NVT or NVE ensemble ?

Indeed, you are correct in pointing out that the time average of the bond length is very close to the imposed bond length.


Is there a reason why the rescaling of the bonds and angles are applied only in the NPT ensemble, and not the NVT or NVE ensemble ?

Yes, in the NPT ensemble pressure is corrected every time step by adjusting the dimensions of the box. It is not a rescaling “of the bonds” but a general rescaling of all the atoms, changing all their relative distances. This is why I mentioned the rescaling “during NPT simulations”, it is a general feature. You’ll find more information about this in the fix npt documentation or the relevant literature. As this is pretty basic in modern MD simulation, I suggest you either look for textbooks on the topic or anyone directly able to advise you on how to set up your MD simulations.

The difference between standard fix npt and fix rigid/npt relies essentially on the equations of motions for the atoms/rigid bodies and the way they are coupled to the dynamics of the box, but the way the box changes are handled and atoms are rescaled is the same.

@Germain, indeed it is a rescaling of all atoms, not the bonds themselves !
Thank you for the explanation and the references.