Pin atoms

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) # masa del átomo de argón
system.addParticle(atom.element.mass) # masa del átomo de argón

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.pin_atoms(simulation, selection='chain_name=="B"', 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: 3.28 kJ/mol ± 2.43 kJ/mol

Execution time: 0 days, 0 hours, 0 minutes, and 22.12 seconds (3905.109 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.show()