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:

Trajectory

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.