How to make fix rigid/nvt work compatibly with compute temp/com command, or any other alternative options in non-equilibrium simulations.

        I am running non-equilibrium molecular dynamics (NEMD) simulations for CO₂ flowing in the carbon nanotubes, with CO₂ being treated as rigid three site linear molecule. The carbon nanotubes are treated as rigid as well. An eternal force is applied on each CO₂ molecule to drive the fluid flow, and there hence exists a non-zero streaming velocity for the fluid. The lammps version I am using is 15 August 2015.

      Although I am using the commands below to adjust the temperature of CO₂, I know this is an unjustified method to follow, as the streaming velocity must be deducted before we thermostat the fluid.

fix 2 CO2 rigid/nvt  molecule temp 300.0 300.0 100.0

fix 200 CO2 addforce 0.0 0.0 9.31E-04 every 1 # (for the external driving force)

However, it is given explicitly in the manual that the 2 NVT rigid fixes do not use any external compute to compute instantaneous temperature. This actually emphasis the fact that we cannot use compute temp/com and fix_modify to deduct the streaming velocity when we thermostat the fluid. So, are there any other possible ways to thermostat the fluid properly?

Nevertheless, I tried to change the command fix rigid/nvt list above to the following commands:

fix 2 CO2 rigid/nve molecule

fix 3 CO2 temp/berendsen 300.0 300. 100.0

fix 200 CO2 addforce 0.0 0.0 9.31E-04 every 1 # (for the external driving force)

compute CO2_tempcom CO2 temp/com

fix_modify 3 temp CO2_tempcom

I found, when I adjust the temperature of CO₂ in this, the temperature increases with simulation time. More specifically the temperature goes up to thousands Kelvins when the simulation runs for 100 ns. However, the streaming velocity in this case stays steadily around 40 m/s (i.e. 4.0E-4 fs/A). So, I speculate that the combination of fix rigid/nve and berendsen thermostat do not work well in NEMD simulations, although they may work for equilibrium simulations.

One may suspect it is the fix_modify command that drives the temperature away. The truth is, even if I remove the fix_modify command, the temperature of the system faces the similar problem.

Considering the way I conducted the NEMD simulations, which I believe should have been widely used in many NEMD simulations, it therefore will benefit a wide range of MD simulations if someone can come up with a better solution regarding this. Certainly, if I haven’t made any conceptual mistakes in my simulations, developing the fix rigid source code to cove this type of simulations is the best solution.

I appreciate all the people who spend any effort on this issue, and I look forward for any suggestions.

Best Regards

Lang

> So, are there any other possible ways to thermostat the fluid properly?

The current answer is no, I don't think there is. Most of the rigid body work we have done for flowing systems is for rigid bodies in a background fluid, e.g. LJ spheres. In those cases thermostatting the solvent is typically sufficient. The fix rigid doc page talks about those kinds of scenarios.

Also, if you were using fix shake (not possible for CO2), e.g. for rigid water, I think it should work fine with fix nvt or nvt/sllod and using a temp compute that subtracts out the flow velocity bias.

The fix rigid/small variants allow for thermostatting (both COM translation and rotation) by NVT or Langevin.  Which is fine when there is no bias.

However for a pure rigid body system with collective motion, what is needed (as you say) is a way to subtract out the collective motion (bias) from the COM velocity of each particle before applying the thermostat.

We've thought about this - what is needed I think is a compute temp/rigid that could be assigned to fix rigid/small that it would use to subtract the chosen bias(es). This is somewhat complicated b/c of how the COM info is stored and computed (in parallel) for lots of small bodies.  Rather than have a list of rigid/bias fixes like temp/rigid/com, temp/rigid/partial, temp/rigid/deform, it would probably be best to have a single compute temp/rigid with extra optional keywords that specified the bias of interest.

 But that compute has not yet been written ...

 Steve