# -------------------------------------------------------------
# Initial Configuration
# -------------------------------------------------------------
bulk_configuration = nlread('PolyVinylChloride.hdf5')[-1]

# -------------------------------------------------------------
# OPLS-AA Calculator
# -------------------------------------------------------------
potential_builder = OPLSPotentialBuilder()
calculator = potential_builder.createCalculator(bulk_configuration)
bulk_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------
initial_velocity = ConfigurationVelocities(
    remove_center_of_mass_momentum=True
)

method = NPTMartynaTobiasKlein(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    reservoir_pressure=1*bar,
    thermostat_timescale=100*femtoSecond,
    barostat_timescale=500*femtoSecond,
    initial_velocity=initial_velocity,
    heating_rate=0*Kelvin/picoSecond,
)

constraints = [FixCenterOfMass()]

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=constraints,
    trajectory_filename='PVC_Trajectory.hdf5',
    steps=10000,
    log_interval=1000,
    method=method
)

bulk_configuration = md_trajectory.lastImage()

# -------------------------------------------------------------
# Calculate the Hildebrand Solubility
# -------------------------------------------------------------
analyzer = CohesiveEnergyDensity(md_trajectory)

# -------------------------------------------------------------
# Plot The Results
# -------------------------------------------------------------
model = Plot.PlotModel(ps, Joule**0.5/cm**1.5)
model.axis('x').setLabel('Time')
model.axis('y').setLabel('Hildebrand Solubility')
model.legend().setVisible(False)

line = Plot.Line(analyzer.times(), analyzer.solubilities())
line.setColor('red')
model.addItem(line)

model.setLimits('x', 0*ps, 0.1*ps)
model.setLimits('y', 10*Joule**0.5/cm**1.5, 20*Joule**0.5/cm**1.5)

Plot.save(model, 'CohesiveEnergyDensity_Plot.png')