# -*- coding: utf-8 -*-
setVerbosity(MinimalLog)

# Load the initial configuration from the HDF5 file
bulk_configuration = nlread('PVC_Input.hdf5')[-1]

# Create and set a default OPLS calculator on the configuration
potential_builder = OPLSPotentialBuilder()
calculator = potential_builder.createCalculator(bulk_configuration)
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

# Set-up and run NVT molecular dynamics with a 2fs time step
initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=True,
    enforce_temperature=True,
)

method = NVTNoseHoover(
    time_step=2*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    chain_length=3,
    initial_velocity=initial_velocity,
)

# Set bond constraint to constrain all C-H bonds
bond_constraint = BondConstraint([[Carbon, Hydrogen]])
com_constraint = FixCenterOfMass()

# With 2fs time steps 5,000 MD steps corresponds to 10ps of simulation time.
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[bond_constraint, com_constraint],
    trajectory_filename='PVC_Bond_Constraint.hdf5',
    steps=5000,
    log_interval=500,
    method=method
)

# Save the resulting configuration
bulk_configuration = md_trajectory.lastImage()
nlsave('PVC_Bond_Constraint.hdf5', bulk_configuration)