Querying 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.
[1]:
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.
[2]:
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.
[3]:
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.
[4]:
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)
