# 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.
forcecap_equilibrator = ForceCappedEquilibration()
equilibrated_structure = forcecap_equilibrator.runEquilibration(configuration)
nlsave('ForceCapped_Configuration.hdf5', equilibrated_structure)

# Create a polymer equilibration object and run the equilibration process
polymer_equilibrator = PolymerEquilibrator(
    configuration=equilibrated_structure,
    filename='PolymerEquilibrationRestart.hdf5',
    object_id='PolymerEquilibration'
)
polymer_equilibrator.update()

# Return the final equilibrated configuration and save it to disk
final_configuration = polymer_equilibrator.equilibratedConfiguration()
nlsave('Final_Configuration.hdf5', final_configuration)