Dear Lammps users,
I am simulating a system of water and hydrogen molecules. The force field of water is TIP4P and the force field of hydrogen is based on papers.
Equilibration: NPT
Production: NVT
The results of the center of mass of water and hydrogen are uploaded. The center of mass is nearly constant in equilibration (slight increase), However, it decreases in the NVT ensemble.
The articles have considered only the first two nanoseconds of the production run for MSD analysis. I plotted the first 2 ns COM and also a plot of 8ns COM to see whether COM is constant or change versus time. when I continue the simulation till 8ns, COM reaches negative values! I don’t really know what the meaning of negative COM is. also, I don’t see a specific problem with visualizing the system.
A short video of the system is attached:
when I simulate water molecules alone, the COM does not change versus time in the NVT ensemble. However, I don’t know what is the problem with hydrogen molecules that changes the center of mass.
I would appreciate any help.
I feel like Axel’s answer to this post applies perfectly here as well.
This make no sense, the COM of a fluid will always move with time, unless you are applying extra constraint to the system. Without your input file, its impossible to tell.
From the PDF this seems like a classic case of the “flying ice cube” effect, since you have a linearly accelerating water centre of mass. But without seeing your input script and your LAMMPS version it’s hard to even start giving you any concrete advice about what to do.
Please note that any force field can only accurately describe a mixture of molecules when it is fitted for that mixture of molecules. You cannot expect any generic force field for molecular hydrogen to accurately describe molecular hydrogen in water. It’s also highly likely that your timestep is too long, depending on your hydrogen model.
Many thanks for your reply
The “flying ice cube” could be caused by a floating point match with finite precision. Could you please guide me on how to increase the floating point number to check whether this behavior is related to that or if something else is wrong?
You cannot change the floating point precision (unless you are running with GPU or INTEL package acceleration).
What you can do instead is to reset any residual center of mass momentum after equilibration with velocity all zero linear. A net center of mass momentum can form most likely at the beginning of a simulation when there are large forces.
If there is still a small buildup of a center of mass momentum during a run (most common with inhomogeneous systems) you can cancel it very few thousand steps or so with fix momentum. Calling fix momentum too often will suppress fluctuations that are required for correct sampling of the desired statistical mechanical ensemble.
When I search for “flying ice cube effect” on DuckDuckGo, every result that comes up on the first page warns about using the wrong algorithms and I don’t see anything about floating point accuracy. So I would be very surprised if it is a matter of switching precision.
As in my previous answer, I strongly suspect that you are using a timestep that is too large for accurate modelling of hydrogen (which is the lightest element and therefore would have larger accelerations). A simple diagnostic would be to set the hydrogen masses to a suitably large number, like 12 amu instead of 1 amu, and see what happens. (The problem of fast hydrogen motions is why water molecules are so often modelled as rigid bodies.)
Also, again as in my previous answer, it is almost impossible for anyone to know what’s going on in your simulation without seeing at least your simulation script. For example, are you using a rigid bond length in your hydrogen molecules or a movable bond length?