validating OpenMM energy terms

[luke czapla]

Hi,

I sort of moved some of the private functions into public space to get the energy terms one by one, was there an easier way to do it? After moving "Context::getImpl" and "*Impl::getImplInContext(Context *)" into public space I do something like this after every 1000 MD steps or so to print out my numbers:

HarmonicBondForceImpl& bondimpl = (HarmonicBondForceImpl&)bonded->getImplInContext(*context);
HarmonicAngleForceImpl& angleimpl = (HarmonicAngleForceImpl&)angle->getImplInContext(*context);
PeriodicTorsionForceImpl& torsionimpl = (PeriodicTorsionForceImpl&)torsion[0]->getImplInContext(*context);
CustomTorsionForceImpl& impropertorsionimpl = (CustomTorsionForceImpl&)impropertorsion->getImplInContext(*context);
NonbondedForceImpl& nbimpl = (NonbondedForceImpl&)nonbonded->getImplInContext(*context);
ContextImpl& cimpl = context->getImpl();
if (usegbsa) {
GBSAOBCForceImpl& gbsaimpl = (GBSAOBCForceImpl&)gbsa->getImplInContext(*context);
cout << "GBSA ENERGY = " << KcalPerKJ*gbsaimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
}
cout << "BOND ENERGY = " << KcalPerKJ*bondimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
cout << "ANGLE ENERGY = " << KcalPerKJ*angleimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
cout << "DIHEDRAL ENERGY = " << KcalPerKJ*torsionimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
cout << "IMPROPER DIHEDRAL ENERGY = " << KcalPerKJ*impropertorsionimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
if (hascmap) {
CMAPTorsionForceImpl& cmaptorsionimpl = (CMAPTorsionForceImpl&)cmaptorsion->getImplInContext(*context);
cout << "CMAP DIHEDRAL ENERGY = " << KcalPerKJ*cmaptorsionimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;
}
cout << "NONBONDED ENERGY = " << KcalPerKJ*nbimpl.calcForcesAndEnergy(cimpl, false, true, 0xff) << endl;


I'm trying to rigorously validate the OpenMM numbers against CHARMM and AMBER and it looks very good so far. The GBSAOBCForce and NonbondedForce are the tough ones where I need to figure out the decomposition but the rest agree with "Sander" in AMBER (single point calculation) and with the NAMD and CHARMM simulation using explicit solvent.

Luke


P.S. - OpenMM seems to crash with just one water molecule in the System. I print out these energies and they are all extremely tiny (indeed there's no interaction because of constraints and exclusions and its just very small bonds and angle energies, the Urey Bradley is in there and is a constraint) but I get "Segmentation fault: 11" even if I skip minimization and go right to "integrator->step(Nsteps)" in my code. I can write an exception for one water molecule and it's for water GCMC, but perhaps this is a bug?

[Lee-Ping Wang]

Hi Luke,

It is possible to do an energy decomposition analysis using force groups. Basically, set each Force in your system to a different force group number and then calculate the energy with only that force group activated. The drawback is that this will disable any other use of force groups (in MTS integrators for example).

I think it would benefit OpenMM to have a more "official" way to do energy component analysis, and one that is able to decompose the nonbonded interactions into VdW and electrostatic (perhaps real space and reciprocal space) contributions.

Thanks,

- Lee-Ping

[Peter Eastman]

Also, the code you've written will not work correctly on all platforms. There are multiple ways a force can return its contribution to the energy. There are good reasons those methods are not part of the public API.

Peter

[luke czapla]

Hi Luke,

It is possible to do an energy decomposition analysis using force groups. Basically, set each Force in your system to a different force group number and then calculate the energy with only that force group activated. The drawback is that this will disable any other use of force groups (in MTS integrators for example).

I think it would benefit OpenMM to have a more "official" way to do energy component analysis, and one that is able to decompose the nonbonded interactions into VdW and electrostatic (perhaps real space and reciprocal space) contributions.

Thanks,

- Lee-Ping

Thanks Lee-Ping,

Peter is right, it just returns all zeros when I run on OpenCL, but it did the job to look at the energies on CPU. I will try it with the force groups instead. The problem is mainly the NonbondedForce with GBSAOBCForce (or GBVIForce), I'm trying to do comparison with CHARMM and AMBER and see if I can adapt these forces or CustomGBForce to reproduce some of the models. NonbondedForce just returns one energy so that's where the decomposition might help. I have setReactionDielectricField(1.0) as suggested by the manual with GB and CUDA seems to work with GBSAOBC and GBVI and NonbondedForce is perfect for the explicit solvent and I can compare it against any of the programs there.

It seems GBSW (a CHARMM model) is something more like what I need but I'm also working on using the OpenMM with a ligand for scoring with MM/GBSA. The protein doesn't move here and I just shuffle states from an existing simulation to estimate the ensemble free energy, so I might be able to get away with a less accurate treatment but it's still a matter of adjusting the GB terms and validating it against some existing model. It looks like I'm getting closer and so I will keep looking at it...

AMBER doesn't have the ACE approximation for the SASA that's used in GBSAOBCForce but it has everything else, and CHARMM has this ACE approximation so I may need to validate against a combination or find out how to make CHARMM behave like GBSAOBCForce and GBVIForce. I found "customgbforces.py" to reproduce AMBER models but it seems quite a bit slower for some reason.

Thanks,
Luke

[Peter Eastman]

On the CPU platform, CustomGBForce will be very very slow. That's because it has to run an interpreter to evaluate the expressions. On OpenCL or CUDA, it generates GPU kernels to evaluate them efficiently, so it's much faster there.

Peter

[luke czapla]

Oh I see, I was told by Jason Swails in Dave Case's lab that there is a GBOBC2Force in the works and available for reference and CPU? This one would be a bit different than the GBSAOBCForce and I am just interested in trying it out for validation of the modeling against AMBER Sander [seems the GBSW is not yet ready for testing but I have CHARMM c35b5 to validate these models too].

I moved some of this stuff into force groups for the energy evaluation and found ParmEd, and next I'm trying to implement the dummy particles to take up space in the context for insertion and deletion of molecules. The current version is more of a hack than anything because I am deleting stuff off the end of the list and calling reinitialize(). It's for the bit with the GCMC using MD to integrate the system for sampling.

[Peter Eastman]

I don't know what you mean by "GBOBC2Force". Perhaps you're thinking of the GBn2 model Jason added? You'll find that in customgbforces.py.

Peter

[luke czapla]

I don't know what you mean by "GBOBC2Force". Perhaps you're thinking of the GBn2 model Jason added? You'll find that in customgbforces.py.

Peter

Thanks, he was actually saying that there was no need for customgbforces.py for this one because there was OBC2 somewhere with CPU and Reference platform implemented and it was something you did, but it is not important at this point - I figured I will not run anything until I fully understand the theory behind it and the implementation so I'm just studying these problems.

It looks like I can work on some of these other GB models that haven't been done in OpenMM yet and learn what really goes into them and there's a couple really good ones for CHARMM, so I've been studying what goes into these methods because people have parameterized the models very accurately for proteins but it doesn't seem to have the acceleration yet - the custom expressions seems like an interesting solution to do the easiest problems to the hardest ones and so I will just start with this customgbforce concept where it can be written into these expressions and I have a ton of CustomNonbonded's and things that are generally called "pair potentials" in other programs.

I'm sort of listening to what people need both around here and around the community because I think OpenMM would be great to implement everything and my calculations (besides the ones I do here for my group) are less important than the big picture - I can't even apply for funding for a lot of this stuff because I've been working as a postdoc too long and took a break from everything for almost a year.

Luke

[Peter Eastman]

Code: Select all

because there was OBC2 somewhere with CPU and Reference platform

Oh, he's referring to GBSAOBCForce. It's uses the OBC-II parameters. And it's implemented on all platforms, not just CPU and Reference. But for CUDA and OpenCL, the speed penalty of using CustomGBForce instead of the hand-coded version is much much smaller - maybe 2x rather than 100x for CPU.

Peter

[luke czapla]

Code: Select all

because there was OBC2 somewhere with CPU and Reference platform

Oh, he's referring to GBSAOBCForce. It's uses the OBC-II parameters. And it's implemented on all platforms, not just CPU and Reference. But for CUDA and OpenCL, the speed penalty of using CustomGBForce instead of the hand-coded version is much much smaller - maybe 2x rather than 100x for CPU.

Peter

Great, yes I'm perfectly willing to throw the factor of two to make things work in here.. they're using LAMMPS to run a lot of stuff and I think with the Custom forces and expressions and Python it could be made much easier in OpenMM and there's nothing special in these models here besides different functions that are available in the expressions and I can validate energy values against everything we have.

I just have to port some code bits over that were done by hand in C++, because some of this code was done very ad hoc, but I have people who want to use it and it looks like putting in a lot of the stuff that it really needed [distance-dependent dielectrics was one simple thing for electrostatics] will be a lot easier with OpenMM. Turns out there are grad students over at Rutgers who would like to try these old models, so I just put them on my XSEDE because I would call this educational use of resources.



Luke