[lammps-users] fix rigid, thermostatting and barostatting

hi all,

My apologies on the previous incomplete post on the same topic, that email was accidentally sent.

I am trying to do molecular dynamics on a collection of rigid molecules, using 7 July 2009 version of LAMMPS. The molecules are composed of charged point masses and the intra-molecular interactions are turned off. Inter-molecular interactions are computed as a sum of Buckingham potential and long range Coloumbic interactions. I have a few concerns which i am listing below.

  1. I am using fix rigid to update the positions, velocities and orientation of the rigid bodies in my simulation. What algorithm does LAMMPS use to update the orientations? (Evans quaternion or similar?)

  2. According to the doc files, the way to control the temperature of such a system is through a Langevin thermostat. But the doc page for the fix_langevin explicitly states that this fix modifies only translational velocities. I understand that one can modify the temperature compute used by this fix to account for both rotational and translational KE but is it possible to modify this fix to scale both the translational and rotational velocities?

  3. Trying to barostat this sytem using a fix_press/berendsen results in undesirable results (blows up). I understand that LAMMPS uses a default value for bulk- modulus and this might not be good for my system. I am still to investigate it. Looking in the mailing lists, the latest post on this topic ( http://sourceforge.net/mailarchive/message.php?msg_id=33ea1b0a0909140717k1f06104ah93f6ce54cb5861dd%40mail.gmail.com) suggests new capabilities being added.

  4. Fix_rigid is not invoked during minimization. So are the atoms being individually moved, if i try to minimise the system?

Any comments on these issues will be greatly appreciated. Thanks in advance.

Answers below.

Steve

hi all,

My apologies on the previous incomplete post on the same topic, that email
was accidentally sent.

I am trying to do molecular dynamics on a collection of rigid molecules,
using 7 July 2009 version of LAMMPS. The molecules are composed of charged
point masses and the intra-molecular interactions are turned off.
Inter-molecular interactions are computed as a sum of Buckingham potential
and long range Coloumbic interactions. I have a few concerns which i am
listing below.

1) I am using fix rigid to update the positions, velocities and orientation
of the rigid bodies in my simulation. What algorithm does LAMMPS use to
update the orientations? (Evans quaternion or similar?)

It represents the rigid bodies with quaternions. I don't know what "Evans"
means. The integration is done with a Richardson iteration (see the
fix_rigid.cpp
richardson() and intial_integrate() routines).

2) According to the doc files, the way to control the temperature of such a
system is through a Langevin thermostat. But the doc page for the
fix_langevin explicitly states that this fix modifies only translational
velocities. I understand that one can modify the temperature compute used by
this fix to account for both rotational and translational KE but is it
possible to modify this fix to scale both the translational and rotational
velocities?

In this context, fix langevin is being applied not to the rigid bodies,
but to each atom in the rigid body. This effectively randomizes both
the translation and rotational DOF, and does a reasonable job of
equilibrating systems I've worked with.

3) Trying to barostat this sytem using a fix_press/berendsen results in
undesirable results (blows up). I understand that LAMMPS uses a default
value for bulk- modulus and this might not be good for my system. I am still
to investigate it. Looking in the mailing lists, the latest post on this
topic (
http://sourceforge.net/mailarchive/message.php?msg_id=33ea1b0a0909140717k1f06104ah93f6ce54cb5861dd%40mail.gmail.com)
suggests new capabilities being added.

Fix rigid corrrectly computes the virial contribution due to the rigid bodies,
so I don't know why it would blow up, unless you try and compress the
system too much. I've successfully run NPT to shrink a box around an
initially very dilute system.

(hopefully) soon, we will release some code that does NVT and NPT directly
on the rigid body degrees of freedom - this is due to some code provided
by U Michigan folks.

4) Fix_rigid is not invoked during minimization. So are the atoms being
individually moved, if i try to minimise the system?

yes, you can't run minimize with rigid bodies.