Goals

Placed in a rough guess as to order of importance. Please comment on order.

1. reduce the space explored during optimization
   • can be achieved by anything that forces the particles to stay rigid
2. accelerate non-bonded list computations
   • requires some sort of hierarchical proximity detection scheme
3. reduce the work done by the optimizer
   • for general purpose optimizers, need to allow the optimizer to work directly on rigid bodies by either computing derivatives or implementing rigid body mechanic
   • for MD, this requires a way for the optimizer to figure out what particles represents a rigid body

Desired Functionality

• ability to change the composition of the rigid body during optimization
  • Manually
  • Automatically, with some criteria ( rmsd between last stored rigid bodies is more than X-- DSR: I don't understand)

Existing Modeller Approach

• A simple interface providing a list of atoms in each body is available.
• Calculation of a non-bonded restraint: All intra-body atom pairs are removed from the nonbonded list, since those distances cannot change. Thus these atom pairs will no longer contribute to any nonbonded interactions, such as Coulomb or Lennard-Jones interactions.
• Optimizers:
  • The MD optimizer is "aware" of the rigid bodies, and does quaternion-based integration.
  • All other optimizers have no knowledge of rigid bodies. They optimize an N-vector (where N is the number of degrees of freedom: 3 per atom and 6 per rigid body) of numbers (system vector state).

  • Certain algorithms like CG work best when all of numbers in the system state verctor have a similar range. Since in the rigid case those numbers are not all xyz coordinates (atom coordinates plus coordinates of the center of mass of each rigid body are, but the 3 Euler angles for each body are not), various schemes have been proposed to force the angles into a similar range to the xyz coordinates. One could imagine scaling the angles by the radius of gyration of each body. Other labs do a dynamic rescaling using the second derivatives. In Modeller we obtained reasonable results using a simple user-tunable scaling parameter (the scale_factor option).

  • An option to tune the scale_factor member , which is used to scale the system state vector (used by CG and quasi_newton optimizers) to rigid body orientation Euler angles (in radians). (Note that no scaling is done for the position of the rigid body; thus the units of this factor are effectively radians/Å.) This can improve optimization convergence in some cases. By default the scaling factor is 1.0; values larger than 1 increase the rotational sampling, while values less than 1 will decrease it.

Current/Proposed IMP Solution

• RigidBodyDecorator and associated helpers allows one to define a rigid body

• in optimization rigid bodies can be handled two ways
  1. as a post optimization step, compute a rigid alignment between the pre-opt and post-opt configuration of a rigid body and fix deviations of the particle from the rigid_body "constraint".
  2. propagate the derivatives from the member particles to the position and orientation of the rigid body. Allow the optimizer to optimize the position and orientation and then copy things back.

Informal experiments so far indicate that convergence is quite slow in general and that things work better for method 1 (the tests fail less often)

• There is a wide space to explore with regards to speeding up collision detection
  • ignore the rigid bodies

  • create a PairScore which uses rasterized versions of the rigid body to directly compute some score

  • use a hierarchical collision detection scheme (in a PairScore) to perform collision detection.

  • make the ClosePairsScoreState handle certain pairs differently and create a CloseBipartitePairsScoreState for them and add the results to the list

• is it worth having a scaling factor (or more intuitively, a range) for each type of optimized attribute so that they can be scaled internally in optimizers which care? That is Model::set_attribute_range(FloatKey, Float max_value).

• it is currently hard to tell if a particle is part of a rigid body as you need to know the rigid body traits used before you check (with is_instance_of). This is not easily fixed without removing the ability for a particle to be part of multiple rigid bodies.