SimTk/openmm — pull request #157 by rmcgibbo (2013-10-02 02:05:04+00:00)
Title: [WIP] Desmond System Test

Included in this PR are two system tests that compare OpenMM and Desmond in two ways. This isn't actually meant to be merged though.

**(1a)** Use desmond's tools to create a "desmond molecular system" (DMS) file that includes all of the parameters mapped onto all of the pairs, triplets, etc of atoms for ala2 in amber99sb-ildn in tip3p. Load this DMS file in OpenMM and calculate the forces on every atom in the system.
**(1b)** Take the same geometry of ala2/tip3p, but setup the system in OpenMM from a pdb via OpenMM's `amber99sbildn.xml`

**These forces agree within 0.1 kcal/mol/nm**. I have done this test also on a larger system, and I know there are two sources of error. The first is that I think the templating of certain improper torsions that are declared in the amber forcefields with wildcards, like *, *, CT, NH is done slightly differently between OpenMM and `viparr`, which is desmond's system setup tool. The second source of error, I think, which is much smaller and the one that is active for this ala2 system, comes from some numerical issues with the nonbonded exceptions. This could be in 2 parts. One is that the setup of nonbonded exceptions requires converting coefficients for LJ forces specified with `a` and `b` as in `a/r^12 + b^6` to `sigma` and `epsilon` in the more "standard" way of writing the LJ force: this could introduce some small numerical issues. The second part of this error is that desmond lists the amber fudgeQQ parameter as 0.8333, whereas OpenMM uses two extras threes (0.833333). So there's some going to be some small difference there.

For the second test, I setup the system with desmond's templating engine (as in 1a), and then
**(2a)** Calculate the forces with desmond, by calling the desmond binary with the dms file.
**(2b)** Calculate the forces in OpenMM, from the dms file.

**Here, the forces do not agree**. But furthermore, the forces calculated *with desmond* using differtent PME cutoffs don't agree with one another. For example, here are the forces on the first three atoms in the system calculated with different PME cutoffs in OpenMM and Desmond. If you state at these numbers a little bit, you'll see that in the desmond numbers, there is a strong trend in which each of the forces move monotonically (by quite a bit) as you increase the cutoff. So clearly something is wrong with my desmond config file. If anyone has any ideas (cc: @leeping @kyleabeauchamp @jchodera) I'm all ears.

```
OpenMM nonbondedCutoff = 8 A
[[ 126.87218047 136.04450917 24.66186398]
[ 39.37468144 478.242654 68.87983042]
[ 32.79352221 -57.76740471 -53.52405496]]
OpenMM nonbondedCutoff = 9 A
[[ 127.02834782 136.24668641 24.6223779 ]
[ 38.54905691 478.211485 69.27124718]
[ 32.81578078 -57.57349163 -53.109029 ]]
OpenMM nonbondedCutoff = 10 A
[[ 126.9447592 136.19952351 24.71190588]
[ 38.64393074 478.064535 69.1117237 ]
[ 32.84153748 -57.85112691 -53.31015571]]
OpenMM nonbondedCutoff = 11 A
[[ 126.95114064 136.11638617 24.67819236]
[ 38.39931538 478.19056683 69.13630179]
[ 32.76650786 -57.72095102 -53.2609079 ]]
OpenMM nonbondedCutoff = 12 A
[[ 126.97312378 136.18475034 24.84544781]
[ 38.33276344 478.436878 68.93574197]
[ 32.82778515 -57.77204697 -53.30481262]]

Desmond nonbondedCutoff = 8 A
[[ 128.29093933 123.14955139 30.95908928]
[ -13.34107208 504.97229004 38.50177002]
[ 41.41397095 -12.21356869 -48.0520668 ]]
Desmond nonbondedCutoff = 9 A
[[ 125.7900238 127.39439392 31.84701729]
[ -7.42175341 501.32952881 36.75513458]
[ 45.25564957 -5.83741856 -47.85729599]]
Desmond nonbondedCutoff = 10 A
[[ 123.88157654 130.65797424 33.0093956 ]
[ -3.13715291 498.11022949 34.66875458]
[ 48.03310013 -1.72448182 -48.11376953]]
Desmond nonbondedCutoff = 11 A
[[ 122.67459106 133.14651489 34.34104919]
[ -0.69900793 495.22064209 32.1003685 ]
[ 50.16262054 1.20843768 -48.27210617]]
Desmond nonbondedCutoff = 12 A
[[ 121.91383362 135.07867432 35.75953674]
[ 0.50310475 492.65002441 29.18287659]
[ 51.81376266 3.28319073 -48.2746048 ]]
```