# 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')