BondConstraint¶

class BondConstraint(fixed_bonds=None, fixed_angles=None, inversion_order=4, rotation_order=1)¶

This constraint fixes the bond lengths of the system in a MolecularDynamics() simulation or OptimizeGeometry() calculation based on the LINCS method.

Parameters:

fixed_bonds (list) – The bonds to constrain. These are specified as a list of pairs. If the pair contains two integers, these are taken as the indices of the bond. If the pair contains two strings, these are taken as tags describing two sets of atoms between which the bonds are to be constrained. If the pair contains two elements, these are taken as the elements between which all bonds are to be constrained.
fixed_angles (list) – The angles to constrain. These are specified as a list of triples. The angle is constrained by constraining the two bonds and the distance between the end atoms in the angle. In this way small molecules such as water can be constrained to be rigid. If the triple contains to integers, these are taken as the indices of the angle. If the triple contains two strings, these are taken as tags describing three sets of atoms between which the angles are to be constrained. If the triple contains three elements, these are taken as the elements between which all angles are to be constrained.
inversion_order (int) – The number of terms used to approximate the matrix inversion in the LINCS equations.
Default: 4
rotation_order (int) – The number of cycles used in the correction for rotational bond lengthening.
Default: 1

bondIndices()¶

Returns:: List of indices of the bonds that are constrained
Return type:: ndarray

bondLengths()¶

Returns:: List of the constrained bond lengths
Return type:: PhysicalQuantity of type length

frozenDegreesOfFreedom(local_atoms=None)¶

The number of constrained degrees of freedom.

Parameters:: local_atoms (list of int | None) – The group of atoms from which the frozen degrees of freedom should be calculated, e.g. a thermalized group of atoms.
Default: All atoms.
Returns:: The number of degrees of freedom that are frozen by this constraint object.
Return type:: int

uniqueString()¶: Return a unique string representing the state of the object.

Usage Examples¶

Run a molecular dynamics calculation on a PVC polymer with constrained bonds for all hydrogen atoms. This allows using a longer time-step of 2 femtoseconds.

method = NVTNoseHoover(
    time_step=2*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    chain_length=3,
    initial_velocity=initial_velocity,
)

# Set bond constraint to constrain all C-H bonds
bond_constraint = BondConstraint([[Carbon, Hydrogen]])
com_constraint = FixCenterOfMass()

# With 2fs time steps 5,000 MD steps corresponds to 10ps of simulation time.
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[bond_constraint, com_constraint],
    trajectory_filename='PVC_Bond_Constraint.hdf5',
    steps=5000,
    log_interval=500,
    method=method
)

PVC_Bond_Constraint.py PVC_Input.hdf5

Notes¶

The BondConstraint class implements the LINCS bond constraint [1] for both molecular dynamics and geometry optimization simulations. In the LINCS algorithm bonds are held rigid by introducing bond distance constraints \(g\):

\[g_i(\mathbf{r}_n) = \left| \mathbf{r}_{i_1} - \mathbf{r}_{i_2} \right| - d_i = 0 \; \; \; i=1,.....K\]

These are then included in the equations of motion using Lagrange multipliers \(\lambda_i(t)\) such that:

\[-\mathbf{M} \frac{d^2\mathbf{r}}{dt^2} = \frac{\partial}{\partial \mathbf{r}} \left( \mathbf{V} - \mathbf{\lambda} \cdot \mathbf{g} \right)\]

These equations of motion are solved at each step of the simulation, allowing the changes in the selected bond lengths to be projected out of the atomic coordinates. Bond length changes are also projected out of the velocities, forces and stresses calculated during the simulation.

These bond constraints can be used with atoms at relative distances to each other during a simulation. This can aid in locating transition states or other intermediate configurations. The bond constraints can also be used to remove high frequency motion during molecular dynamics. In these kinds of simulations, the necessary time step of the simulation is often determined by the high frequency bond vibrations. This is especially true for bonds containing light elements such as hydrogen. By constraining these bonds it can be possible to increase the time step used for the molecular dynamics simulation, as those high frequency modes no longer have to be represented in the dynamics.

In BondConstraint, the bonds to be constrained are specified in the argument fixed_bonds. This takes a list of pairs of either elements, tags or atom indices as integers. For elements and tags, all static bonds between those elements or pairs will be constrained. If atom indices are given, the distance between these atoms will be constrained whether or not they are bonded. The constrained bond length will be taken from the starting configuration in the simulation.

It is also possible to constrain bond angles using the argument fixed_angles. This argument takes a list of triples of elements, tags or bond indices. The atoms specified in this triple are then constrained in a rigid triangle, approximating a rigid angle. Regardless of the type in input used, a distance constraint is added between the two end atoms whether or not they are bonded. This type of constraint can be used to create rigid small molecules, such as rigid water.

Warning

The solution to the constraint equations is performed using an approximate matrix inverse. Including angle constraints makes this approximation less valid, and increases the number of terms that need to be included to get a good solution. If angle constraints are included the parameter inversion_order may need to be significantly increased to account for this.

In the LINCS algorithm there are two iterative parts to the solution of the constraint equation. The first is the solution of a matrix inverse, which is approximated as a series. The number of terms in this series is modified using the argument inversion_order. Secondly the LINCS algorithm corrects changes in bond length due to rotation using an iterative solution. The number of iterations is modified using the argument rotation_order. By default these have values of 4 and 1 respectively. This is suitable for most calculations with just constrained bonds. Including constrained angles or large number of bonds may require increasing these values to obtain stable constraints.