Problem with NEMD Simulation of Rigid Particles in a LJ Solvent System

Dear LAMMPS Developers and Users:

Hope you all are doing well!

Currently, I am trying to run an NEMD simulation of a system comprising of LJ particles (background solvent) and rigid particles (collection of LJ particles resembling a linear rod of L = 6σ). I went over the LAMMPS manual for doing this particular task. I found that the nvt/sllod command must be only applied to the solvent that will thermostat the system at the desired temperature and the fix deform command (with the remap v option) to all the particles in the system (including the rigid bodies). I also used fix/rigid command for the rigid bodies to only allow for the center-of-mass translation and rotation.

When I ran the simulations (for both moderate and low shear rates upto 0.0008), I got the errors ‘rigid body atoms missing’ or ‘atoms lost’ after many time-steps (I tried using comm_modify too and it did not work). Particularly, I get this error for all the shear rates when the tilt factor reached half the box length (maximum allowed), 25σ in my case (i.e. my box length is 50σ, so the tilt factor will vary between -25σ to 25σ). The energy, maximum force (Fmax), temperature, and pressure is in control as desired. The manual also says that the velocities are not remapped for the rigid particles even though they are in the group defined for the fix/deform command. I am not sure whether the remapping of velocities of the rigid bodies, the real issue for missing atoms when the tilt factor becomes 25 and tries to vary in the other direction (say 25, 24…,0,…,-25,…0,…25).

Also, I tried to reduce the time step since the tilt factor is calculated as (T0 + L0eratedt). I tried this because it helped me in getting the simulations rightly done at high shear rates for the pure LJ system, say erate = 0.7. To be specific, at erate = 0.7 and dt = 0.003, I got the error ‘Atoms missing’. But, when I reduced the time-step to dt = 0.0015. I was able to run the simulations at this shear rate. I believe that allowing the tilt factor to change at a slower pace does matter a lot but it did not work for the systems with rigid particles.

Further, I find an answer by Matt Peterson (Oct 22, 2007, forwarded by Dr. Plimpton) about the velocity remapping of rigid bodies in the old mailing list (see below, in bold text) which is contradictory to the present LAMMPS manual. He says that the rigid body velocities is also remapped upon using the remap v option (not what the manual says). Everything else is in agreement with what I have done. He also says that he saw some instabilities at high shear rates (not sure of what that means). Please correct me if I am missing anything on my side.

<b>Here was Matt Peterson's answer:

Hi David.  So, to answer your questions; you should apply fix deform to
all particles and remap velocities, use the sllod thermostat for the
solvent only, and *don't* apply fix NVE to the rigid bodies.

The deformation of the box using slodd requires that all particles
(including the rigid bodies) have their velocities remapped when
crossing the sheared boundary.  Particles at this boundary "see" the
particles across the boundary moving relative to the shear rate, and
without the remapping would crash into the opposing streaming particles
when they cross.  The remapping adjusts their velocity to match the
streaming velocity across the boundary.  The rigid fix implies an NVE
integrator, so none is needed for the rigid bodies - slodd should only
be applied to the solvent.  The sheared boundaries and velocity
remapping will take care of shearing the rigid bodies, and the sllod
thermostatted solvent should be enough to control any heating due to the
shearing.  I've found some instabilities at high shear rates; maybe due
to a poorly chosen timestep, or maybe because the system is started so
far from the assumed velocity profile.  Also, take care to adjust the
degrees of freedom if you are using liner rigid bodies.  The rigid
documentation explains how an why.  And use neighbor modify to exclude
interactions between particles in a rigid body.  It wont change the
dynamics, but the pressure and temperature could be wrong.  Feel free to
email me if you have more questions.</b>

Could someone please help me get answers for the questions posted above. Is the velocity not being remapped for the rigid particles the issue or is it something else. That being said, is reverse NEMD (Muller-Plathe technique) a good option for these systems?

Any suggestions or answers to these queries are highly appreciated.

Have a great day!

Thanks and Regards,

Dinesh

1 Like