SoftMatterDynamicsSimulation

Included in QATK.Dynamics

class SoftMatterDynamicsSimulation(configuration, method=None, atomic_constraints=None, bond_constraints=None)

Create a class for running soft matter dynamics simulations.

Parameters:

  • configuration (MoleculeConfiguration | BulkConfiguration) – The starting configuration for the simulation

  • method (BaseMDMethod) – The dynamics method.
    Default: NVEDynamics()

  • atomic_constraints (None | FixAtomConstraints | FixCenterOfMass | list of (FixAtomConstraints | FixCenterOfMass)) – The constraints to add to atoms

  • bond_constraints (None | HydrogenBonds | AllBonds | HydrogenBondAngles | AllBondHydrogenAngles) – The bonds to constrain in the configuration. Values can be None, HydrogenBonds to constrain bonds to hydrogen atoms, AllBonds to constrain all bonds, HydrogenBondAngles to constrain all hydrogen bonds and H-X-H angles and AllBondHydrogenAngles to constrain all bonds and H-X-H and H-O-X angles. The latter two options also result in rigid water molecules.

currentConfiguration(velocities=True)

Get the current configuration from the dynamics simulation.

Parameters:

velocities (bool) – Include the velocities in the configuration

Returns:

The current configuration

Return type:

MoleculeConfiguration | BulkConfiguration

relax(max_force=PhysicalQuantity(0.05, eV / Ang), max_steps=200, report_frequency=1)

Relax the dynamics system to the given steps and tolerance.

Parameters:

  • max_force (PhysicalQuantity of type force) – Stop relaxation when maximum atomic force is below this value.
    Default: 0.05 eV/Ang

  • max_steps (int | None) – Maximum relaxation cycles. None means no iteration limit and the optimization continues until the force limit is reached.
    Default: 200

  • report_frequency (int | None) – Step interval between reporting the relaxation progress. None means that the relaxation runs without writing to the log.
    Default: 1

Returns:

Whether the relaxation is converged.

Return type:

bool

setLogging(log_interval)

Set how often simulation information is written to the log.

Parameters:

log_interval (int) – The number of steps between writing log information.

setTrajectory(trajectory_interval, trajectory_filename=None, trajectory_object_id=None, write_velocities=None, write_forces=None, write_stress=None)

Set writing the simulation to a trajectory.

Parameters:

  • trajectory_interval (int) – The step frequency of writing trajectory images.

  • trajectory_filename (str) – The name of the trajectory file to write to.

  • trajectory_object_id (str | None) – The object id of the trajectory in the file.

  • write_velocities (bool) – If True, velocities will be written to the trajectory.
    Default: True

  • write_forces (bool) – If True, forces will be written to the trajectory.
    Default: False

  • write_stress (bool) – If True, stress will be written to the trajectory.
    Default: False

simulate(total_steps, profiles=None, hook_functions=None)

Perform a given number of simulation steps.

Parameters:

  • total_steps (int) – The total number of steps to take

  • profiles (SimulationQuantityProfile | list of SimulationQuantityProfile | None) – Profiles that change variables during the simulation.

  • hook_functions (SoftMatterDynamicsHook | list of SoftMatterDynamicsHook | None) – Hook functions that are called during the simulation.

trajectory()

Return the MDTrajectory used for storing the trajectory.

Returns:

The trajectory object, or None if no trajectory is being written.

Return type:

MDTrajectory | None

Usage Examples

Calculate the pressure in an NVT simulation of poly(vinyl chloride), measuring fluctuations in the dipole moment, This simulation uses bond constraints on hydrogen bonds to increase the timestep to 2fs.

# Create simulation
method = NVTNoseHoover(
    time_step=2 * fs,
    initial_velocity=ConfigurationVelocities(),
)
simulation = SoftMatterDynamicsSimulation(
    configuration, method, bond_constraints=HydrogenBonds
)
simulation.setLogging(1000)
simulation.setTrajectory(1000, "pvc_trajectory.hdf5")

# Create a measurement of the pressure
hook = SoftMatterDynamicsMeasurements(dipole=All, call_interval=100)

# Run the simulation
simulation.simulate(100000, hook_functions=[hook])

equilibrated_pvc.hdf5 pvc_simulation.py

Notes

The SoftMatterDynamicsSimulation class implements molecular dynamics and relaxation using the SoftMatterCalculator. This calculator implements bonded forcefields such as OPLS and DREIDING with algorithms optimized for running on GPUs. This gives these dynamics improved performance over those performed with the MolecularDynamics() function and the TremoloXCalculator and the cost of some reduced flexibility in how the calculation is performed. In particular, the SoftMatterDynamicsSimulation only supports hook functions that run in-between molecular dynamics steps, instead of during as in the MolecularDynamics() method. For calculations that fit within the assumptions of the SoftMatterDynamicsSimulation it offers significantly improved GPU performance over the TremoloXCalculator

To create a SoftMatterDynamicsSimulation both an initial configuration and a molecular dynamics method are required. The molecular dynamics method details the specific parameters of the dynamics. This includes things like initial velocities, time steps and how temperature and pressure might be controlled. Note that not all methods that work with MolecularDynamics() are also compatible with the SoftMatterDynamicsSimulation. Differences in the two forcefield implementations can provide some impediments to providing full compatibility and support for all methods. The table below lists each molecular dynamics method and whether it is available using the TremoloXCalculator or the SoftMatterCalculator

Table 34 Molecular dynamics methods and which calculators they are compatible with.

MD Method

Ensemble

TremoloX

Soft Matter

NVEDynamics

NVE

Yes

Yes

Langevin

NVT

Yes

Yes

NVTBerendsen

NVT

Yes

No

NVTBussiDonadioParrinello

NVT

Yes

No

NVTNoseHoover

NVT

Yes

Yes

NVTAndersen

NVT

No

Yes

NPTBerendsen

NPT

Yes

No

NPTMartynaTobiasKlein

NPT

Yes

No

NPTBernettiBussi

NPT

Yes

No

NPTLangevinMonteCarlo

NPT

No

Yes

NPTNoseHooverMonteCarlo

NPT

No

Yes

NPTAndersenMonteCarlo

NPT

No

Yes

One notable difference between the TremoloXCalculator and SoftMatterDynamicsSimulation is in how pressure and stress are calculated. While the TremoloXCalculator calculates these analytically, the SoftMatterCalculator calculates them numerically using a finite difference method. This makes barostats that calculate the stress and every step, such as the NPTBernettiBussi barostat, difficult to implement. The SoftMatterCalculator instead supports a Monte Carlo type barostat. Here random simulation cell changes are attempted and accepted or rejected based on the overall change in enthalpy. By default these attempted cell changes are done every 25 steps. Over long simulation runs the isothermal-isobaric ensemble is effectively sampled. The barostat can be implemented independently of the thermostat, and so different combinations of barostat and thermostat are available.

In addition to the molecular dynamics method, atomic and bond constraints can also be provided. These are considered as part of the dynamics system, and are used in all dynamics simulations. To change constraints a new simulation object must be created. Atom constraints are constraints that apply to individual atoms. These are the FixAtomConstraints and FixCenterOfMass constraints. Note that when fixing atomic positions all cartesian directions must be constrained. Bond constraints apply to specific bonds within the configuration. Here these bonds are fixed to their equilibrium bond length set by the forcefield. One common use of this constraint is to freeze high frequency motions in the dynamics, allowing for the use of longer time steps. As an exmaple, when all bonds involving hydrogen atoms are constrained the dynamics time step can be increased from 1fs to 2fs. There are 4 levels of bond constraints available, each selected with a flag. The flag HydrogenBonds constrains all bonds involving hydrogen. The flag AllBonds increases the constraints to all bonds in the configuration. Angle constraints can be added with the flag HydrogenBondAngles which constrains all hydrogen bonds and H-X-H angles. Finally the flag AllBondHydrogenAngles constrains all bonds and H-X-H and H-O-X (i.e. alcohol) bond angles. Using one of the latter two flags means e.g. all water molecules are kept rigid.

Once the simulation object is created the logging and trajectory writing during the simulation are specified with the setLogging and setTrajectory methods. By default the simulation does not write anything. Images are written to the trajectory during the simulation, so if a trajectory file is given this is written to during the simulation.

The simulation is then carried out using the simulate method. This takes as a required argument the number of steps to simulate. It also optionally takes a list of profiles and hook functions. Profiles are used to continuously change a property of the simulation. Values such as reservoir temperature, reservoir pressure, simulation cell length and NPT Monte Carlo sampling frequency can be all specified with profiles. Hook functions are used to either change the simulation or to extract information from it in a more complicated way. Examples of hook functions include the SoftMatterDynamicsNonEquilibrium and SoftMatterDynamicsMeasurements classes. At the end of a simulate call the state of the simulation is preserved in the SoftMatterDynamicsSimulation object. Subsequent calls to simulate will then continue the simulation from this point allowing the initial simulation to be extended. Both profiles and hook functions are considered properties of the simulation call. This means they are only active during the simulate calls they are passed into. When using multiple simulate calls the hook functions and profiles can be changed between simulations.

In addition to molecular dynamics, atomic relaxation can be performed with the SoftMatterDynamicsSimulation object. This is primarily designed to remove large atomic forces before running a dynamics simulation. Atomic relaxation is performed using the relax method. This optionally takes a few arguments. The max_force argument specifies the convergence limit for the relaxation. Likewise the max_steps argument gives the number of steps before non-convergence is reported. The frequency of writing log information is controlled with the report_frequency argument. Note that if None is specified for this the relaxation will not produce any log output. The relax method returns a bool indicating if the optimization converged.

The current configuration in the simulation can be returned using the currentConfiguration method. This is typically used after a simulation to get the final configuration. The current trajectory can also be returned using the trajectory method.