# Load the configuration
configuration = nlread("equilibrated_urethane.hdf5")[0]

# Create and set the calculator
builder = OPLSPotentialBuilder()
builder.addCustomOplsTorsion(
    "HC", "C", "N", "CA", 2.3*kiloCaloriePerMol, 6.09*kiloCaloriePerMol
)

calculator = SoftMatterCalculator(
    builder, use_lennard_jones_pme=True, assign_atom_types=True
)
configuration.setCalculator(calculator)

# Set the simulation
atomic_constraints = [FixCenterOfMass()]
method = NPTAndersenMonteCarlo(
    time_step=2*fs,
    initial_velocity=ConfigurationVelocities(),
    reservoir_temperature=300*Kelvin,
    reservoir_pressure=[10, 1, 1] * bar,
    pressure_coupling=DiagonalOnly,
)
simulation = SoftMatterDynamicsSimulation(
    configuration,
    method,
    atomic_constraints=atomic_constraints,
    bond_constraints=HydrogenBonds
)
simulation.setLogging(1000)
simulation.setTrajectory(10000, "urethane_trajectory.hdf5")

# Run the simulation
simulation.simulate(100000)