CrosslinkBuilder¶
- class CrosslinkBuilder(crosslink_connector, restart_trajectory=None, skip_initial_crosslink=None, nvt_steps_per_crosslink=None, npt_steps_per_crosslink=None, temperature=None, pressure=None, isotropic_pressure=None, cycles_per_md_step=None, cycles_to_terminate=None, reaction_percentage=None, bond_search_increment=None, max_bond_search_radii=None, use_fbmc_relax=None, fbmc_steps_per_crosslink=None, h_passivate_types=None, potential_builder=None, opls_type_list=None, **kwargs)¶
Creates a polymer crosslink builder object.
- Parameters:
crosslink_connector (
CrosslinkConnector
) – The polymer cross-linker that defines the possible reactions.restart_trajectory (
Trajectory
| None) – A trajectory from which to continue the cross-linking reaction. Default: None.skip_initial_crosslink (bool) – Whether or not to skip the initial cross-linking step when starting the simulation. Default: False.
nvt_steps_per_crosslink (int) – The number of NVT molecular dynamics steps to perform per cycle. Default: 10000.
npt_steps_per_crosslink (int) – The number of NPT molecular dynamics steps to perform per cycle. Default: 10000.
temperature (PhysicalQuantity of type temperature) – The temperature that the cross-linking is performed at. Default: 300K.
pressure (PhysicalQuantity of type pressure) – The pressure that the cross-linking is performed at. Default: 1 bar.
isotropic_pressure (bool) – Whether or not pressure is applied isotropically. Default: True.
cycles_per_md_step (int) – The number of cross-linking cycles performed between NPT-MD cycles. Default: 1.
cycles_to_terminate (int) – The number of cycles the simulation will continue without any cross-linking reactions. None does not apply this termination condition. Default: None.
reaction_percentage (float) – The target reaction percentage at the completion of the reaction. Default: 100.
bond_search_increment (PhysicalQuantity of type length) – The increment to add to the bond search radius when no reactions are found in a cross-linking step. Default: 0.25 Angstrom.
max_bond_search_radii (PhysicalQuantity of type length) – The maximum bond search radius to be reached during the cross-linking process. Default: 10 Angstrom.
use_fbmc_relax (bool) – Whether or not to use force-bias Monte Carlo to relax the system after cross-linking. By default NVT molecular dynamics with a gradual relaxation of the reaction coordinates is used. Setting this variable to True changes the relaxation step to a combination of force-bias Monte Carlo and NVT molecular dynamics. Default: False.
fbmc_steps_per_crosslink (int) – The number of force-bias Monte Carlo steps to perform per cycle. This step is only performed if force-bias Monte Carlo relaxation is selected with the input use_fbmc_relax. Default: 10000.
h_passivate_types (str | list | None) – List of tags that specify atoms to be hydrogen passivated at the end of the cross-linking reaction. Default: None.
potential_builder (
DreidingPotentialBuilder
|UFFPotentialBuilder
|OPLSPot entialBuilder
) – The potential builder used to create the potential for the cross-linking system. The default is to use the Dreiding potential without electrostatic interactions.opls_type_list (list of str | None) – Optional list of all of the OPLS atom types that can arise during cross-linking. Works only if the OPLSPotentialBuilder is used, and can increase the performance in this case. Default: None.
- runCrosslinking(trajectory_filename=None)¶
Run the cross-linking simulation.
- Parameters:
trajectory_filename (str | None) – The file where the result trajectory is written.
- Returns:
The trajectory of the system as it is being cross-linked.
- Return type:
Usage Examples¶
Build a network polymer system based on an amine-epoxy reaction.
# Load the configuration to be reacted.
configuration = nlread('Amine_Epoxy_Crosslink.hdf5')[-1]
# Specify the reactions for primary amines with an epoxy using the two different hydrogens.
reaction_a = CrosslinkReaction(
add_bonds=[('REACT_A_C', 'REACT_B_N')],
remove_bonds=[('REACT_A_C', 'REACT_A_O')],
transfer_groups=[('REACT_B_H1', 'REACT_B_N', 'REACT_A_O')],
secondary_reaction_map={'REACT_B_N': 'REACT_C_N', 'REACT_B_H2': 'REACT_C_H', },
exclude_rings_below_size=0,
bond_cutoff=4*Angstrom
)
reaction_b = CrosslinkReaction(
add_bonds=[('REACT_A_C', 'REACT_B_N')],
remove_bonds=[('REACT_A_C', 'REACT_A_O')],
transfer_groups=[('REACT_B_H2', 'REACT_B_N', 'REACT_A_O')],
secondary_reaction_map={'REACT_B_N': 'REACT_C_N', 'REACT_B_H1': 'REACT_C_H', },
exclude_rings_below_size=0,
bond_cutoff=4*Angstrom
)
# Specify the reaction of a secondary amine with an epoxy.
reaction_c = CrosslinkReaction(
add_bonds=[('REACT_A_C', 'REACT_C_N')],
remove_bonds=[('REACT_A_C', 'REACT_A_O')],
transfer_groups=[('REACT_C_H', 'REACT_C_N', 'REACT_A_O')],
exclude_rings_below_size=0,
bond_cutoff=4*Angstrom
)
# Define the CrosslinkConnector, using the defined reactions and functional groups to find the
# reactive atom groups in the configuration.
crosslink_connector = CrosslinkConnector(
configuration,
reactions=[reaction_a, reaction_b, reaction_c],
functional_groups=[['REACT_A_C', 'REACT_A_O'], ['REACT_B_N', 'REACT_B_H1', 'REACT_B_H2']],
max_number_of_reactions=100
)
# Define the CrosslinkBuilder using the default Dreiding potential
polymer_builder = CrosslinkBuilder(
crosslink_connector,
cycles_to_terminate=None,
reaction_percentage=100,
nvt_steps_per_crosslink=1000,
npt_steps_per_crosslink=1000,
temperature=480*Kelvin,
pressure=1*bar,
bond_search_increment=0.5*Angstrom,
max_bond_search_radii=6*Angstrom
)
# Run the polymer building simulation.
polymer_builder.runCrosslinking(trajectory_filename='Amine_Epoxy_Thermoset.hdf5')
Amine_Epoxy_Crosslink.py
Amine_Epoxy_Crosslink.hdf5
Build a system of linear polymers of poly(ethylene oxide).
# Load the configuration to be reacted.
configuration = nlread('PEO_Crosslink.hdf5')[-1]
# Specify the reactions of ethylene glycol using both possible atom arrangements.
reaction_a = CrosslinkReaction(
add_bonds=[('REACT_C_A', 'REACT_O_B')],
remove_atoms=['REACT_O_A', 'REACT_H_A', 'REACT_H_B'],
bond_cutoff=5*Angstrom
)
reaction_b = CrosslinkReaction(
add_bonds=[('REACT_O_A', 'REACT_C_B')],
remove_atoms=['REACT_O_B', 'REACT_H_A', 'REACT_H_B'],
bond_cutoff=5*Angstrom
)
# Explicitly define the reactive atom groups in the configuration.
reactive_atom_groups = [{'REACT_C_A': 0, 'REACT_O_A': 2, 'REACT_H_A': 8},
{'REACT_C_A': 1, 'REACT_O_A': 5, 'REACT_H_A': 9},
{'REACT_C_A': 10, 'REACT_O_A': 12, 'REACT_H_A': 18},
{'REACT_C_A': 11, 'REACT_O_A': 15, 'REACT_H_A': 19},
{'REACT_C_A': 20, 'REACT_O_A': 22, 'REACT_H_A': 28},
{'REACT_C_A': 21, 'REACT_O_A': 25, 'REACT_H_A': 29},
{'REACT_C_A': 30, 'REACT_O_A': 32, 'REACT_H_A': 38},
{'REACT_C_A': 31, 'REACT_O_A': 35, 'REACT_H_A': 39},
{'REACT_C_A': 40, 'REACT_O_A': 42, 'REACT_H_A': 48},
{'REACT_C_A': 41, 'REACT_O_A': 45, 'REACT_H_A': 49},
{'REACT_C_A': 50, 'REACT_O_A': 52, 'REACT_H_A': 58},
{'REACT_C_A': 51, 'REACT_O_A': 55, 'REACT_H_A': 59},
{'REACT_C_A': 60, 'REACT_O_A': 62, 'REACT_H_A': 68},
{'REACT_C_A': 61, 'REACT_O_A': 65, 'REACT_H_A': 69},
{'REACT_C_A': 70, 'REACT_O_A': 72, 'REACT_H_A': 78},
{'REACT_C_A': 71, 'REACT_O_A': 75, 'REACT_H_A': 79},
{'REACT_C_A': 80, 'REACT_O_A': 82, 'REACT_H_A': 88},
{'REACT_C_A': 81, 'REACT_O_A': 85, 'REACT_H_A': 89},
{'REACT_C_A': 90, 'REACT_O_A': 92, 'REACT_H_A': 98},
{'REACT_C_A': 91, 'REACT_O_A': 95, 'REACT_H_A': 99},
{'REACT_C_B': 0, 'REACT_O_B': 2, 'REACT_H_B': 8},
{'REACT_C_B': 1, 'REACT_O_B': 5, 'REACT_H_B': 9},
{'REACT_C_B': 10, 'REACT_O_B': 12, 'REACT_H_B': 18},
{'REACT_C_B': 11, 'REACT_O_B': 15, 'REACT_H_B': 19},
{'REACT_C_B': 20, 'REACT_O_B': 22, 'REACT_H_B': 28},
{'REACT_C_B': 21, 'REACT_O_B': 25, 'REACT_H_B': 29},
{'REACT_C_B': 30, 'REACT_O_B': 32, 'REACT_H_B': 38},
{'REACT_C_B': 31, 'REACT_O_B': 35, 'REACT_H_B': 39},
{'REACT_C_B': 40, 'REACT_O_B': 42, 'REACT_H_B': 48},
{'REACT_C_B': 41, 'REACT_O_B': 45, 'REACT_H_B': 49},
{'REACT_C_B': 50, 'REACT_O_B': 52, 'REACT_H_B': 58},
{'REACT_C_B': 51, 'REACT_O_B': 55, 'REACT_H_B': 59},
{'REACT_C_B': 60, 'REACT_O_B': 62, 'REACT_H_B': 68},
{'REACT_C_B': 61, 'REACT_O_B': 65, 'REACT_H_B': 69},
{'REACT_C_B': 70, 'REACT_O_B': 72, 'REACT_H_B': 78},
{'REACT_C_B': 71, 'REACT_O_B': 75, 'REACT_H_B': 79},
{'REACT_C_B': 80, 'REACT_O_B': 82, 'REACT_H_B': 88},
{'REACT_C_B': 81, 'REACT_O_B': 85, 'REACT_H_B': 89},
{'REACT_C_B': 90, 'REACT_O_B': 92, 'REACT_H_B': 98}]
# Define the CrosslinkConnector, using the defined reactions and reactive atom groups.
crosslink_connector = CrosslinkConnector(
configuration,
reactions=[reaction_a, reaction_b],
reactive_atom_groups=reactive_atom_groups,
max_number_of_reactions=5
)
# Define the CrosslinkBuilder object, using termination after 3 steps without reactions.
polymer_builder = CrosslinkBuilder(
crosslink_connector,
nvt_steps_per_crosslink=1000,
npt_steps_per_crosslink=1000,
temperature=350*Kelvin,
pressure=90*bar,
cycles_per_md_step=2,
cycles_to_terminate=3,
bond_search_increment=0.0*Angstrom,
)
# Run the polymer building simulation.
polymer_builder.runCrosslinking(trajectory_filename='PEO_Polymer.hdf5')
PEO_Crosslink.py
PEO_Crosslink.hdf5
Note that in both these cases the examples have been created to be illustrative of the process of building the material. For high-quality simulations it may be necessary to increase the number of molecular dynamics steps or, as in the case of PEO, the size of the system.
Notes¶
The CrosslinkBuilder
class manages the building and cross-linking of polymer
systems by iteratively performing both cross-linking reactions and also controlling the motion of
atoms between each cross-linking step [1][2].
In each cycle, relaxation of both the atomic positions and simulation cell can be performed. There
are two different relaxation protocols implemented. The default relaxation protocol is based on
molecular dynamics. First the position of the transfer atoms are relaxed with a constrained geometry
optimization. After this the atomic positions are relaxed using NVT molecular dynamics. To avoid
large forces on the reacted atoms, the reaction coordinates are also gradually scaled from reactant
to product. The density of the system can then be further relaxed in an NPT molecular dynamics
simulation. Using molecular dynamics in both stages allows for the preservation of velocities
between cross-linking steps. An alternative relaxation method using Force-Bias Monte Carlo (FBMC)
is also available. In this protocol, all atomic positions are initially optimized with a geometry
optimization before being relaxed in an FBMC step. Because FBMC is not sensitive to large
forces, no scaling of the reaction coordinates is necessary in this step. An NVT molecular dynamics
step is then run before relaxing the density with NPT molecular dynamics. In this case, equilibrated
velocities are lost during the FBMC step, and so additional NVT is required before NPT molecular
dynamics to first equilibrate the velocities.
On initialization, the CrosslinkBuilder
requires a CrosslinkConnector
that contains the details of the possible reactions. There are also a number of optional arguments
that control how the configuration is evolved between cross-linking steps. The number of NVT and
NPT steps in each molecular dynamics cycle is specified with the arguments
nvt_steps_per_crosslink
and npt_steps_per_crosslink
. If FMBC based relaxation is desired,
this can be accessed by setting the argument use_fbmc_relax
to True
. The number of FBMC
steps can then be set with fbmc_steps_per_crosslink
. The temperature and pressure used in the
dynamics simulations can be set with the arguments temperature
and pressure
respectively.
The argument isotropic_pressure
also sets whether pressure is applied isotropically, or
independently for the three cell lengths. This can be useful in cases where the material is
anisotropic, such as in building thermoset materials at a solid surface. It is also possible to
skip some NPT molecular dynamics steps at each cycle. This is done with the argument
cycles_per_md_step
. This takes an integer that specifies how many reaction and NVT molecular
dynamics steps are done for every NPT molecular dynamics step. Skipping NPT molecular dynamics
steps can be appropriate in cases where the density of the system does not significantly change with
the progress of the reaction.
There are three separate termination criteria that can be set to determine when the cross-linking
reaction process is complete. First, the argument cycles_to_terminate
specifies a number of
steps in which the simulation will terminate if no reactions are performed. If a reaction is
performed at any step the counter is returned to zero. The default value of None
allows for
the simulation to continue through any number of unreactive steps. Secondly, the simulation can
also terminate when a defined reaction extent is reached, which is specified through the argument
reaction_percentage
. The default value allows for the reaction to go to completion. Finally the
reaction can be terminated when a specific bond search radius is reached. Using the argument
bond_search_increment
, an increase in the bond search radius used to determine possible
reactions can be added if no reactions are found during a reaction step. Increasing the bond search
radius during the simulation prioritizes reactions between closer reactive groups
[3]. When the bond search radius reaches the limit specified by
max_bond_search_radii
, the simulation terminates.
To set the potential used to determine the motion of the atoms during the simulation, a
PotentialBuilder
object can be passed in that determines the forcefield used at each step. As
the configuration is reacting, the object used must be able to automatically assign atomic
forcefield types for the configuration. By default, a Dreiding potential ignoring electrostatic
interactions is used. It should be noted that in both the UFFPotentialBuilder
and
DreidingPotentialBuilder
, electrostatic interactions are calculated in each cycle using the QEq
method. This may significantly slow down simulations in network materials which cannot be easily
decomposed into a number of molecules. For the OPLS potential it is also possible to specify a list
of types for the atoms in different stages of the reaction. This is intended to limit the search
range that the OPLSPotentialBuilder
must use to determine the types and speed up type
assignment. If no list is given, the OPLSPotentialBuilder
will search though all available types
to determine the best fits of types to the structure.
Partially complete simulations can be restarted using the argument restart_trajectory
. This
takes a Trajectory
from a previous calculation created by a CrosslinkBuilder
.
Restarting the calculation from the trajectory resets the count of the extent of the reaction and
the bond search radius. Note that if a restart trajectory is given, the last frame of the
trajectory should also be used to define the CrosslinkConnector
object used to create the
CrosslinkBuilder
. When restarting a calculation it is also possible to skip the
initial cross-linking step. This is useful in cases where the simulation terminated because of a
problem with a reaction. Skipping the initial cross-linking step in the simulation allows for more
movement of the atoms and can resolve problems with the reaction. Skipping the initial
cross-linking step is done by setting the argument skip_initial_crosslink
to True
.
After the reaction is completed the final configuration can be passivated with hydrogen. The
argument h_passivate_types
takes a list of tags of atoms that are to be passivated. This
requires that the tags be present in the final configuration. This is useful in cases such as in
the reaction of acrylates, where radicals may be left at the conclusion of the simulation. After
hydrogen atoms are added, a geometry optimization is then performed on the resulting configuration.
Once the CrosslinkBuilder
object is defined, the actual cross-linking simulation
is performed by calling the runCrosslinking
method. This returns a trajectory of the
resulting configuration after every geometry relaxation cycle, as well as the final structure after
possible passivation. The argument trajectory_filename
allows the specification of a file
location where the trajectory can be written. If this is not given, the trajectory is stored in
memory before being returned.