ForceCappedEquilibration

class ForceCappedEquilibration(temperature=None, md_steps_per_force_capped_simulation=None, md_time_step=None, force_capped_simulations=None, starting_factor=None, ending_factor=None, fixed_indices=None, lennard_jones_cutoff=None)

Force-capped equilibration for polymer melt systems.

Parameters:
  • temperature (PhysicalQuantity of type temperature) – The temperature at which the force-capped equilibration will be run.
    Default: 300 Kelvin

  • md_steps_per_force_capped_simulation (int) – The number of MD steps that are run at each level of force capping.
    Default: 40000

  • md_time_step (PhysicalQuantity of type time) – The timestep used in the molecular dynamics simulation.
    Default: 0.5 fs

  • force_capped_simulations (int) – The number of force-capped simulations. Each simulation is run at a different force-capping factor.
    Default: 4

  • starting_factor (float) – The starting position for the inclusion of force capping, expressed as a fraction of the sigma value in each potential.
    Default: 1.04

  • ending_factor (float) – The ending position for the inclusion of force capping, expressed as a fraction of the sigma value in each potential.
    Default: 0.8

  • fixed_indices (list of int | None) – A list of indices that should be fixed during the force-capped equilibration.

  • lennard_jones_cutoff (PhysicalQuantity of type length) – The cutoff beyond which non-bonding interactions are ignored.
    Default: 8 Angstrom

runEquilibration(configuration, initial_bond_only_step=False, trajectory_filename=None, log_interval=None)

Run the force-capped equilibration protocol on the given configuration.

Parameters:
  • configuration (BulkConfiguration) – The starting configuration for the equilibration.

  • initial_bond_only_step (bool) – Run an initial equilibration cycle without any non-bonding interactions.
    Default: False

  • trajectory_filename (str | None) – The name of the output trajectory file.

  • log_interval (int | None) – The frequency of saving configurations to output trajectory.
    Default: Number of MD steps

Returns:

The equilibrated configuration.

Return type:

BulkConfiguration

Usage Examples

Create a polymer model of poly(methyl methacrylate) and use the ForceCappedEquilibration object to remove atomic overlaps in the structure created by the PolymerMonteCarloBuilder.


# Load the monomer for the polymer
monomer = nlread('Methyl_Methacrylate.hdf5')[-1]

# Define the polymer sequence as containing 10 chains of 100 monomers each, joined atactically.
sequence = PolymerSequence(
    [monomer],
    number_of_monomers=100,
    number_of_chains=10,
)

# Build the initial polymer structure to a correct density for the given polymer.
builder = PolymerMonteCarloBuilder(
    sequence,
    density=1.18*gram*cm**-3,
)
configuration = builder.buildPolymerConfiguration()

# Place an OPLS potential on the configuration as the basis of the force capped equilibration.
potential = OPLSPotentialBuilder()
calculator = potential.createCalculator(configuration)
configuration.setCalculator(calculator)

# Equilibrate the structure and save the resulting configuration to disk.
equilibrator = ForceCappedEquilibration()
equilibrated_structure = equilibrator.runEquilibration(configuration)
nlsave('Final_Configuration.hdf5', equilibrated_structure)


ForceCappedEquilibration_Example.py Methyl_Methacrylate.hdf5

Notes

When creating realistic models of polymer systems, it is important to ensure that the internal degrees of freedom in the polymer chain are properly randomized. In QuantumATK, this is achieved with the PolymerMonteCarloBuilder. Once the chains are generated, it is then necessary to pack them correctly into a unit cell so that none of the atoms overlap. Overlapping atoms generally result in very large interatomic forces that make molecular dynamics and other types of simulations impossible. To ensure that the polymer chains are well packed, one approach is to use the ForceCappedEquilibration class [1].

The ForceCappedEquilibration class contains an equilibration protocol that combines optimization and molecular dynamics. This separates overlapping polymer chains and returns a structure that is ready for further simulation. In order to do this the input configuration needs to contain a TremoloX calculator that contains LennardJonesSplinePotential potentials. Both the DreidingPotentialBuilder and the OPLSPotentialBuilder create potentials that contain LennardJonesSplinePotential potentials to model van der Waals non-bonding interactions.

The force-capped equilibration protocol avoids the large forces associated with overlapping or close atoms by assuming a constant non-bonding interaction force at particle distances below a certain switching value. An example of such as potential is shown in the figure below, along with the reference Lennard-Jones potential. Here force-capping has been applied where the potential crosses the x-axis. Beyond this point both potentials are identical at larger separation.

../../../_images/Force_Capped_Graph.png

In QuantumATK the Lennard-Jones function is defined as:

\[V(r) = 4\epsilon \left(\left[\frac{\sigma}{r}\right]^{12} - \left[\frac{\sigma}{r}\right]^6\right)\]

Here \(\sigma\) is the position at which the potential is zero and \(\epsilon\) is the depth of the potential minimum. The potential minimum is located at \(2^{1/6}\sigma\). In a molecule the values for \(\sigma\) and \(\epsilon\) usually differ between pairs of atomic types to reflect differing chemical environments. Since \(\sigma\) varies between potentials, the position where force-capping is applied is defined relative to the \(\sigma\) value for that potential. During equilibration, this position where force-capping is applied is reduced over a number of cycles. This makes the applied potential increasingly similar to the full Lennard-Jones potential as the equilibration continues and the atomic overlaps are resolved. As electrostatic interactions are difficult to treat in this manner, any defined electrostatic interactions in the input configuration are ignored during force-capped equilibration and can be switched on after the force-capped equilibration.

The number of force-capped equilibration cycles is controlled by the parameter force_capped_simulations. A new force-capped potential is used in each cycle of equilibration. The range of force-capping factors is controlled by the input parameters starting_factor and ending_factor, which select the starting and final ratios of sigma. If force_capped_simulations is larger than two, the remaining factors are evenly spaced between the starting and ending factors.

During each equilibration cycle a geometry optimization, followed by an NVT molecular dynamics run, are performed. These can be modified using some additional parameters. The parameter temperature sets the temperature at which the molecular dynamics are carried out. This allows the simulation temperature to be set to the desired temperature of the final polymer model. The number of molecular dynamics steps in each cycle is controlled with the parameter md_steps_per_force_capped_simulation. The longer each MD simulation is run, the slower the non-bonded interactions are introduced, which can reduce deformation of the chain structure for long polymer melts [1].

Once an instance of the ForceCappedEquilibration class is defined, the equilibration is carried out using the runEquilibration function. As well as an input configuration, this takes an optional input argument initial_bond_only_step. Setting this argument to True causes an additional initial equilibration step to be carried out in which no non-bonding interactions are taken into account. This can be useful in cases where the polymer system has some highly strained internal degrees of freedom, such as torsions or bond angles, that need to be relaxed before packing of the polymer chains can be considered. By default this initial equilibration step is not performed. It is also possible to output trajectories of the molecular dynamics steps using the parameters trajectory_filename and log_interval to set the file name and the number of steps between images respectively.