Energy decomposition of System using CHARMM36 FF

[Justin Chan]

I have been trying to do energy decomposition of a system running CHARMM36 FF on OpenMM.

I am using scripts provided by CHARMM-GUI to run the simulation because they have implemented CustomNonBondedForce to do force-based switching.

After creating the system and adding in the coordinates I ran the following script to get the energy decomposition:

Code: Select all

def getEnergyDecomposition(context, system):
    print(context.getState(getEnergy=True).getPotentialEnergy())

    forces = []

    for i in range(system.getNumForces()):
        current_force = system.getForce(i)
        forces.append(current_force)

        print(current_force, ':', context.getState(getEnergy=True, groups=2**i).getPotentialEnergy())

    return forces

The total potential energy is similar to the potential energy calculated by NAMD. But I am unable to reproduce the total potential energy from the energy decomposition. I get the following output by running the function (first line is the total potential energy):

-1542535.43777 kJ/mol
<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7f5c26948060> > : 76593.1529571 kJ/mol
<simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7f5c26948090> > : 137966.790553 kJ/mol
<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7f5c26948240> > : 150203.559384 kJ/mol
<simtk.openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x7f5c26948120> > : 45353.6148533 kJ/mol
<simtk.openmm.openmm.CustomTorsionForce; proxy of <Swig Object of type 'OpenMM::CustomTorsionForce *' at 0x7f5c269481b0> > : 8677.96254746 kJ/mol
<simtk.openmm.openmm.CMAPTorsionForce; proxy of <Swig Object of type 'OpenMM::CMAPTorsionForce *' at 0x7f5c26948180> > : 2961.66222989 kJ/mol
<simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7f5c269481e0> > : -1935585.15758 kJ/mol
<simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7f5c26948210> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f5c26948270> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7f5c269482a0> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7f5c269482d0> > : 4784.50361701 kJ/mol

I have managed identify each force type to the FF terms:
1. Bond (Excludes Hbonds, using rigid HBonds)
2. Angle
3. Bond (Urey-Bradley term)
4. Dihedral
5. Improper
6. Cross-terms
7. Non-bonded (electrostatics)
8. Custom-nonbonded (van der Waals)
9. Center of mass motion remover
10. Custom-nonbonded (van der Waals)
11. Custom-bonded (van der Waals)

I have attached the python script that implements the force switching for the van der Waals.

Another thing that I noticed is that the number of dihedral angles reported by PeriodicTorsionForce is a lot more than the number report by NAMD/VMD and in my PSF file. Weirdly, when I run this system on AMBER (using CHAMBER to convert PSF) I get the exact same number as in OpenMM.

Thanks in advance!

[Peter Eastman]

I think I know what's happening. It's a combination of two things.

First, I'm pretty sure your force groups aren't defined in quite the way you think they are (one group per Force object). When looping over forces, have it print out current_force.getForceGroup(). I think you'll find some of them are in the same group.

Second, you're running into this bug in 6.3 (but fixed in 7.0): https://github.com/pandegroup/openmm/issues/1400. The effect is that if 1) you're using PME, and 2) you request the energy of a set of force groups that does not include the NonbondedForce, then it will incorrectly include part of the Ewald energy in the value it returns.

There are two clues that led me to these conclusions. First, notice that it's outputting exactly the same value (4784.50361701) for all of the last four groups. That suggests there aren't really that many independent force groups. Furthermore, if you add up all the energies printed by your script, then subtract off the correct total energy (-1542535.43777), the difference is 47845.03: precisely 10 times the value reported for the last four groups. Notice that you query 11 groups, so 10 groups that do not include the NonbondedForce.

Try using 7.0 (just released earlier today), and that error should go away.

Peter

[Justin Chan]

Thank you so much, I managed to get OpenMM to report the energy decomposition and it matches the energy decomposition reported by NAMD!

But I still cannot figure out why I am getting a lot more dihedral angles in OpenMM compared to what is reported by NAMD and in the PSF file (the dihedral section reports the number of dihedrals)

[Jason Swails]

But I still cannot figure out why I am getting a lot more dihedral angles in OpenMM compared to what is reported by NAMD and in the PSF file (the dihedral section reports the number of dihedrals)

When you use getNumTorsions() on the PeriodicTorsionForce in OpenMM, you get a list of the total number of Fourier terms in all of the torsions combined. Remember that many torsions are comprised of several individual terms with different periodicities (i.e., the torsion term is a sum of cosine terms with different parameters). The PSF file (and seemingly NAMD as well) reports each individual torsion as "1" torsion, where each torsion can have many terms.