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.
Hi Axel,
Apologies for leaving this thread hanging for so long.
First, I agree that that 100 step run was pointless and I’m not quite sure why I had that set so low but that’s fixed.
But also, I’ve tried out some of the suggestions you’ve given without much success. I started from a data file instead of binary restart file and added the “mom yes” keyword/value pair to my velocity create command and the system still drifted with the 1fs tilmestep. I then attempted the same thing with a 0.1 fs timestep and still got the drift on approximately the same timescale as using 1 fs. I could go lower but I’m not sure it’ll help.
Something I have found from vaguely related threads is potentially using something along the lines of "fix lnmom all momentum 100 linear 1 1 1 " to reset the linear momentum of the system every so many time steps. Does this seem valid?
Thanks,
Micah
What about using fix shake for your rigid molecules instead of fix rigid/nvt?
The other possibility is using another thermostat – I would try temp/csvr to start with. That would at least rule out it being the force field rather than the thermostat. If it’s the thermostat’s problem and you still want a deterministic thermostat you could try adjusting Nose-Hoover parameters – time constant, chain thermostat etc.
Cheers
Shern