Add plane harmonic restraint

import molsysmt as msm

import openmm as mm
from openmm import unit
from openmm import app
from tqdm import tqdm
import numpy as np
from matplotlib import pyplot as plt

topology = app.Topology()
chain = topology.addChain('A')
residue = topology.addResidue('Ar', chain)
atom = topology.addAtom(name='Ar', element= app.element.argon, residue=residue)
chain = topology.addChain('B')
residue = topology.addResidue('Ar', chain)
atom = topology.addAtom(name='Ar', element= app.element.argon, residue=residue)

system = mm.System()
system.addParticle(atom.element.mass)
system.addParticle(atom.element.mass)

1

temperature = 300*unit.kelvin

integration_timestep = 2.0*unit.femtoseconds
saving_interval = 1.00*unit.picoseconds
logging_interval = 100.00*unit.picoseconds
simulation_time = 1000.*unit.picoseconds

saving_steps = int(saving_interval/integration_timestep)
logging_steps = int(logging_interval/integration_timestep)
md_steps = int(simulation_time/integration_timestep)

friction   = 5.0/unit.picoseconds
integrator = mm.LangevinIntegrator(temperature, friction, integration_timestep)

platform = mm.Platform.getPlatformByName('OpenCL')

simulation = app.Simulation(topology, system, integrator, platform)

initial_positions  = [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]] * unit.nanometers
simulation.context.setPositions(initial_positions)

Lbox = 2.0
v1 = [Lbox,0,0] * unit.nanometers
v2 = [0,Lbox,0] * unit.nanometers
v3 = [0,0,Lbox] * unit.nanometers
simulation.context.setPeriodicBoxVectors(v1, v2, v3)

msm.thirds.openmm.forces.add_plane_harmonic_restraint(simulation, selection='all',
                                                      force_constant='5000 kilojoules/(mol*nanometers**2)',
                                                      vector=[0,0,1], pbc=True)

0

reporter_tqdm = msm.thirds.openmm.reporters.TQDMReporter(logging_steps, md_steps, temperature=False)
simulation.reporters.append(reporter_tqdm)

reporter_trajectory_dict = msm.thirds.openmm.reporters.TrajectoryDictReporter(saving_steps, time=True, coordinates=True)
simulation.reporters.append(reporter_trajectory_dict)

simulation.step(md_steps)

Potential energy: 2.22 kJ/mol ± 2.12 kJ/mol

Execution time: 0 days, 0 hours, 0 minutes, and 14.53 seconds (5945.952 ns/day).

trajectory_dict = reporter_trajectory_dict.finalize()

plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,0,2])
plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,1,2])
plt.axhline(y=0.0, color='red', linestyle='-')
plt.axhline(y=1.0, color='red', linestyle='-')
plt.show()

plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,0,0])
plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,0,1])
plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,1,0])
plt.plot(trajectory_dict['time'], trajectory_dict['coordinates'][:,1,1])
plt.show()