How do we introduce constraints?

[g-bauer]

There are a lot of force fields that represent molecules as semi-flexible or even rigid. To have that functionality in a simulation, constraints have to be introduced.

The complexity that comes along with constraints is very different for MC and MD. While it is trivial to implement e.g. constant bonding lengths in MC, we need to have specific algorithms to do that in MD.

We should collect ideas how we can implement constraints so that they can be used conveniently for both propagators.

Some questions that pop into my mind:

• Where/how to keep track of degrees of freedom?
• Where would we implement information about constraints? Could they be part of Interactions?
• Which functionalities do we want/need to provide? Constant bonding lengths, -angles, completely rigid molecules?
• Where would constraint algorithms be implemented for MD?

[Luthaf]

Where/how to keep track of degrees of freedom?
Where would we implement information about constraints? Could they be part of Interactions?

We could store both of them in the Molecule struct. Something like:

struct Molecule {
    bonds: Vec<Bond>,
    angles: Vec<Angle>,
    rigid_bonds: Vec<Bond>,
    rigid_angles: Vec<Angle>,
}

And then it would be up to the algorithms to obey the constraints. A rigid molecule would have only rigid bonds. I still have to check the MD constraints algorithms to see if this could work.

Which functionalities do we want/need to provide? Constant bonding lengths, -angles, completely rigid molecules?

All of them?

Where would constraint algorithms be implemented for MD?

RATTLE and SHAKE algorithm can be implemented as part of a new Integrator

[g-bauer]

We could store both of them in the Molecule struct.

It is not clear to me: Where would we get e.g the bonding length from? Would we check the bonding potential's equilibrium distance? So for MC, we need to add a harmonic potential to describe the bond and then - if it is fixed - use the bonding length and ignore the force constant?

[Luthaf]

This could be a solution. Or we could store them aside with the Bond, or in a RigidBond struct:

struct Molecule {
    bonds: Vec<Bond>,
    angles: Vec<Angle>,
    rigid_bonds: Vec<Bond>,
    rigid_angles: Vec<Angle>,
}

////////////////////////////////////////////////////////

struct RigidBond {
    usize: i,
    usize: i,
    length: f64
}

struct Molecule {
    bonds: Vec<Bond>,
    angles: Vec<Angle>,
    rigid_bonds: Vec<RigidBond>,
    rigid_angles: Vec<RigidAngle>,
}

And then we can get the bong length from:

• the initial configuration. This might be bad as we might not get the same bond length for every bond;
• a parameter in the forcefield; This might be harder to parse and a bit ackward: why an Harmonic potential specifically?
• a parameter in the [molecules] section.

(I am mainly throwing ideas as they come, most of them could work here 😃).

[g-bauer]

the initial configuration. This might be bad as we might not get the same bond length for every bond;

Yepp. Also, if we want to insert a molecule, we would need a separate configuration.

a parameter in the forcefield; This might be harder to parse and a bit ackward: why an Harmonic potential specifically?

I'll have to read a bit more about the MD constraint algorithms, but I think they need the potential function to compute the Lagrange multipliers that have to be introduced. For MC the potential function does not matter; we only need the equilibrium distance. In MC codes you'll often find the option to apply a harmonic potential that either takes an force constant or a flag fixed as input.

In any case, we need both information since we can only get the underlying potentials of bonds, angles, etc from the mapping via Molecule. It's hard to decide where to put the functionality before we sketch the constraint algorithms.