[lammps-users] Issue with mixing fix rigid with fix nvt (gases with polymers)

Hi there,

I’m working on simulating gas diffusion within polymer materials and I’m seeing bizarre results that I can’t quite figure out with just the manual or past forum posts. I’m attaching the inputs I’ve been using at the bottom of this post.

I’m using the OPLS force field to represent my polymer systems and the TraPPE potential to represent nitrogen gas. I understand that the massless particle at the center of mass of the TraPPE model is problematic and I took the suggestion of a previous forum post and set the mass to be very small (1e-20) instead of just zero.

My procedure thus far has been to apply a fix rigid nvt thermostat to the nitrogen atoms of the three-site nitrogen molecule that I added into the system with the ‘create_atoms’ command and the typical npt fix to the polymer in order to equilibrate the volume of the system. I then switched the npt fix on the polymer to an nvt so now all atoms are thermostated under nvt dynamics. But when I run this I see very strange behavior over pretty short time scales as well as get useless mean-squared displacement data for the gas. The entire system (polymer included) almost moves in unison by a large amount over even a single nano-second instead of the expected behavior of the gas molecules moving through (or relative) to the polymer.

I’m not sure if there’s something obvious I’m missing or if there’s some extra considerations to be made when mixing rigid and non-rigid thermostats in this particular case but any help would be greatly appreciated!

Thanks,
Micah Welsch
University of Kansas

What you describe sounds like the so-called “flying icecube” effect. Somehow your initialization and time integration creates a center of mass momentum and that will not go away.
This can be removed when restarting by using the velocity command to remove to rescale velocities with the “mom yes” keyword/value pair.
I would also suggest to restart after equilibration not from a binary restart but from the data file (so you will re-initialized the nose-hoover chains).

the combination of a rigid with a per-atom integrator is making things difficult. But probably also your timestep. 1fs may be too large for the rigid body time integration. That one is far more sensitive to instabilities than per-atom time integration. plus using fix npt is not making things easier.

the first “run” with only 100 steps is miniscule and thus will not achieve anything.