Download Notebook View in GitHub Open in Google Colab

Querying and Modifying Charges and Other Parameters

Sometimes you want to inspect the charges or other parameters of the particles or bonds in a System. Force field parameters are stored in the Force objects added to a System. As an example, let’s load a PDB file and model it using the Amber14 force field.

from openmm.app import *
from openmm import *

pdb = PDBFile('ala_ala_ala.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology)

Now suppose we want to retrieve the charge of every particle. These are stored in the NonbondedForce object. We can call getForces() to retrieve all forces, then find the one of the correct class.

nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]

Now that we have the NonbondedForce, we can call getParticleParameters() to get the nonbonded parameters (charge, sigma, and epsilon) of each particle. Let’s build a list of all the charges.

charges = []
for i in range(system.getNumParticles()):
    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
    charges.append(charge)

You can record the Lennard-Jones parameters (sigma and epsilon) in the same way, but be aware that some force fields compute Van der Waals interactions differently. If you model your system with CHARMM36, for example, you may find that epsilon is 0 for every atom. In that case, it is using a CustomNonbondedForce to implement the Lennard-Jones force.

You can retrieve bonded parameters in the same way. For example, let’s print out every bond involving particle 0 (the N-terminal nitrogen). This information is stored in the HarmonicBondForce.

bonded = [f for f in system.getForces() if isinstance(f, HarmonicBondForce)][0]
for i in range(bonded.getNumBonds()):
    particle1, particle2, length, k = bonded.getBondParameters(i)
    if particle1 == 0 or particle2 == 0:
        print(f'Particles ({particle1}, {particle2}), length = {length}, k = {k}')

Particles (4, 0), length = 0.1471 nm, k = 307105.5999999999 kJ/(nm**2 mol)
Particles (1, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)
Particles (2, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)
Particles (3, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)

In some cases you may want to modify those parameters to make a bond behave differently from what the ForceField assigned. We show you here that to modify these parameters are also easy.

Let’s say we would want to modify the force constant of the bond connecting particle 1 and 0 to a new number.

for i in range(bonded.getNumBonds()):
    particle1, particle2, length, k = bonded.getBondParameters(i)
    if particle1 == 1 and particle2 == 0:
        bonded.setBondParameters(i, particle1, particle2, length, 2666 * unit.kilojoules_per_mole/unit.nanometer**2)

Now if you query again you can see new parameters.

for i in range(bonded.getNumBonds()):
    particle1, particle2, length, k = bonded.getBondParameters(i)
    if particle1 == 0 or particle2 == 0:
        print(f'Particles ({particle1}, {particle2}), length = {length}, k = {k}')

Particles (4, 0), length = 0.1471 nm, k = 307105.5999999999 kJ/(nm**2 mol)
Particles (1, 0), length = 0.101 nm, k = 2666.0 kJ/(nm**2 mol)
Particles (2, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)
Particles (3, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)

Modifying Force objects will affect any new Simulations or Contexts you create, but it will have no effect on ones that already exist. If you want your modifications to apply to an existing Simulation, you can copy them over by calling bonded.updateParametersInContext(simulation.context).