Recalculate forces on atoms

[Marcelo D Poleto]

Hi,

I am trying to recalculate the forces acting on each atom of my system. I am combining OpenMM and MDTraj to 1) build the System, 2) loop over the trajectory frames and 3) extract the array of Forces and slice it to get only a few atoms of interest. Here is my current code for a single frame:

Code: Select all

print("\n> Setting the system:\n")
print("\t- Reading force field directory...")
charmm_params = (toppar_files_here)

# reading and parsing topology and coordinates
psf = CharmmPsfFile(my_psf_file)
crd = CharmmCrdFile(my_crd_file)

print("\t- Setting box (using user's information)...")
box = 67.160261399015955*angstrom
psf.setBox(box, box, box)

print("\t- Creating system and setting parameters...")
system = psf.createSystem(charmm_params, nonbondedMethod=PME, nonbondedCutoff=1.2*nanometers, switchDistance=1.0*nanometer, ewaldErrorTolerance = 0.0005, constraints=HBonds)

# making sure each ForceGroup has its own index
forcegroups = {}
for i, force in enumerate(system.getForces()):
        force.setForceGroup(i)
        forcegroups[force.getName()] = i

integrator = LangevinIntegrator(temperature, 1/picosecond, dt)
simulation = Simulation(psf.topology, system, integrator)
simulation.context.setPositions(crd.positions)

######################################
# Calculating Energies for entire trajectory
traj = md.load(my_pdb_file, top=my_psf_file)
top = traj.topology
ligand = top.select("resname LYS and residue 1")

#print(str(molecule))
print("\nAtom IDs of selection:", ligand)

forces_dict = {}

print("\n>>> Reading and Calculating:")

# update simulation context with new positions
#simulation.context.setPositions(traj.xyz[frame])

# update simulation context with new box vectors
#a, b, c = traj.unitcell_vectors[frame]
#simulation.context.setPeriodicBoxVectors(*(a,b,c))

# get updated energies
energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
energy = energy.value_in_unit(kilocalories_per_mole) # remove unit
print(energy)


forces = simulation.context.getState(getForces=True).getForces(asNumpy=True)
forces = forces.value_in_unit(kilocalorie/(nanometer*mole))
for id in ligand:
    atomname = top.atom(id).name
    resname = top.atom(id).residue.name
    globalresid = top.atom(id).residue.index + 1

    identifier = resname + "_" + str(globalresid) + "_" + atomname

    print(identifier, forces[id])

The thing is, the forces are different from what CHARMM provides me:

Code: Select all

               Fx              Fy              Fz
OpenMM     327.441942      174.08579      -1.07871011
CHARMM     -32.745453      -17.41029       0.10701801

OpenMM    -103.377675      -37.89938      54.70902991
CHARMM      10.332535        3.79285      -5.47057189

I can understand the negative/positive force convention, but I don't understand the difference in a factor of 10. Does anyone have any clue? Am I missing something in my code?

[Peter Eastman]

I think it's just different units. You're printing the units in kcal/mol/nm. If you instead print them in kcal/mol/A they'll be 10 times smaller. That's probably what CHARMM is using.

[Marcelo D Poleto]

Thanks Peter. I think I stared at the code for so long that I forgot that CHARMM uses Angstrom as default distance unit. That solves the case.

One quick question: if by any chance I decide to decompose the forces by ForceGroup, would it make sense to use the code below?

Code: Select all

forces = simulation.context.getState(getForces=True, groups={forceid}).getForces(asNumpy=True)

And should the sum of the contributions from all ForceGroups be the same as the total force as I was obtaining it previously?

[Peter Eastman]

Correct on both.