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:
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.
In QuantumATK the Lennard-Jones function is defined as:
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.