Get forces#

Get forces.

import molsysmt as msm

import openmm as mm
from openmm import app
from openmm import unit

molsys = msm.convert(msm.systems['Trp-Cage']['1l2y.h5msm'], structure_indices=0)
openmm_topology = msm.convert(molsys, to_form='openmm.Topology')
positions = msm.get(molsys, element='atom', coordinates=True)
positions = msm.pyunitwizard.convert(positions[0], to_form='openmm.unit')

forcefield = app.ForceField("amber14-all.xml", "amber14/tip3p.xml")
system = forcefield.createSystem(openmm_topology, nonbondedMethod=app.NoCutoff,
                                 constraints=app.HBonds)

integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 1*unit.femtosecond)
platform = mm.Platform.getPlatform('CPU')

context = mm.Context(system, integrator, platform)
context.setPositions(positions)

forces = msm.molecular_mechanics.get_forces(context)

forces[100]
Magnitude
[-457.15502669008725 19.59674699306821 287.9360605937297]
Unitskilojoule/(mole nanometer)
forces = msm.molecular_mechanics.get_forces(context, magnitude=True)

view = msm.view(molsys)
view.clear()
msm.thirds.nglview.set_color_by_value(view, values=forces, min_value=0.0,
                                      representation='licorice', cmap='Oranges')
view