Add MD support to rigid bodies

[benmwebb]

Note that we have an existing implementation in Modeller. This assumes that rigid bodies do not overlap, are not nested, and do not contain any non-rigid members, but it should be possible to use a similar approach.

[cgreenberg]

Those assumptions sound reasonable. Do you have to input torque for restraints or is that calculated automatically?

[benmwebb]

Well, these are the assumptions made by the Modeller code. IMP doesn't make the same assumptions, which is why things won't be exactly the same.

Torque is calculated automatically. It's basically the cross product of the force and the internal coordinate vector, IIRC. No need to mess with anything on the restraints side (Modeller automatically excludes intra-rigid body pairs from the nonbonded list, which is the only messing with restraints that it does IIRC, but this is already handled in IMP).

[cgreenberg]

How's the outlook on this issue?

[benmwebb]

Probably won't have time to make much progress on this until Feb-March.

[sethaxen]

@benmwebb:

This assumes that rigid bodies do not overlap, are not nested, and do not contain any non-rigid members

Nesting doesn't seem too difficult to handle, but do you have a sense for how non-rigid members change things?

[sethaxen]

Also, can rigid bodies overlap in IMP?

[benmwebb]

do you have a sense for how non-rigid members change things?

You'd just pretend for MD purposes that the member didn't belong to the rigid body (since non-rigid members can move independently). So if the member is a point, you'd handle it in the same way that points are currently handled. If the member is itself a rigid body, you'd use the new rigid body code to handle it. I think you'd just have to update the local coordinates of the non-rigid member after the move to reflect both its new global coordinates and the new reference frame of the rigid body - otherwise when you move the owning rigid body the member's global coordinates will be overwritten.

can rigid bodies overlap in IMP?

It certainly wouldn't work if you added a particle to two different rigid bodies. I'm not 100% sure that this is checked for, but probably you'd get an error on the second rigid body when it tried to create the rigid body attributes (since the first rigid body added them already).

[sethaxen]

@benmwebb do you have more details on Modeller's implementation of rigid body MD linked to above? The documentation is not explicit regarding how the equations of motion wrt the rotational component of the rigid body is integrated.

The reference given in the documentation (Rapaport, 1997) explains that leapfrog (Verlet) integration is inappropriate for this quaternion representation (due to error accumulation), and they recommended a Predictor-Corrector integrator. From scanning more recent literature on the topic, symplectic integrators like leapfrog are applied when considering the rotation matrix representation (also done in a later section of Rapaport), while other integrators are required when considering the quaternion representation.

[benmwebb]

do you have more details on Modeller's implementation of rigid body MD linked to above?

Sure, in the Sali lab you can look at the source code directly at https://svn.salilab.org/modeller/trunk/src/routines/ (warning: Fortran). I would start with md_routines.F90, then look at rigid_body_routines.F90 and maybe quaternion_routines.F90.

leapfrog (Verlet) integration is inappropriate for this quaternion representation (due to error accumulation)

If you're worried about integrator error in your MD code, IMP is probably the wrong tool for you... but I don't have any particular attachment to Verlet integration. Note that changing the integrator would probably break all of the thermostats though.

[sethaxen]

If you're worried about integrator error in your MD code, IMP is probably the wrong tool for you... but I don't have any particular attachment to Verlet integration. Note that changing the integrator would probably break all of the thermostats though.

To be honest, I'm only indirectly interested in MD; I'm more concerned about IMP.isd.HybridMonteCarlo. The Metropolis criterion will correct for errors in integration so long as the integrator maintains certain properties. For the Metropolis criterion in HMC to be valid, the integrator must be 1) volume-preserving (so that no Jacobian modification to the Hamiltonian is needed) and 2) reversible (so that detailed balance is ensured). Symplectic integrators by definition ensure (1); the leapfrog/velocity Verlet integrators satisfy (2), but I'm not certain that leapfrog applied directly to quaternions this way satisfies either. (Details can be found here)

The good news is that there are symplectic, reversible integration algorithms for quaternions that can be interleaved in a velocity Verlet scheme. Two I'm looking at are "DLM" 10.1063/1.474310 and NO_SQUISH 10.1063/1.1473654.

From examining the RigidBody code, I can't see any explicit calculation/updating of the center of mass outside of the RigidBody.get_initial_reference_frame method that is called when the RigidBody is setup with its members. Then the displacement of the reference frame is set to the center of mass, and all accumulation of member derivatives are done with respect to that center. If the alternate RigidBody setup with the ReferenceFrame3D was used, then the center is likely not at the center of mass at all. While I'm not an expert at rigid body dynamics, the simple form of the equations of motion into displacement and rotation assumes that the reference point is the center of mass, and they would need to change if the center was not indeed the COM. Am I missing something here?

[benmwebb]

While I'm not an expert at rigid body dynamics, the simple form of the equations of motion into displacement and rotation assumes that the reference point is the center of mass, and they would need to change if the center was not indeed the COM. Am I missing something here?

Yes, I believe everything relies on the COM being correct, and yes, rigid bodies can take an arbitrary reference frame (it need not be at the COM). So I guess you'd have to correct for this. You should certainly check with @barakr though, since he already uses rigid body torques in his BD code.

[sethaxen]

It turns out that the propagation of derivatives and torques requires no assumptions regarding the nature of the reference frame (the coordinates don't need to be the COM). The COM assumption is only needed to greatly simplify integrating the equations of motion. After speaking with @barakr, I agree that the most sensible thing is for IMP.atom.BrownianDynamics and IMP.atom.MolecularDynamics to simply assume that the coordinates are the COM. If the user wants to do MD or BD with rigid bodies, it's up to them to ensure that this is set up correctly. We just need to document that.

Quaternion derivatives and torques are now propagated correctly for nested rigid bodies. At the moment, derivatives are not handled properly for non-rigid members (this was the case before). We'll eventually want this supported, but it will probably require some special-casing during integration, as @benmwebb noted above. I believe BD also only works correctly for bodies with only rigid members.