(POSSIBLE BUG): GBn and GBn2 exhibit unexpected behaviors unlike HCT, OBC1, OBC2 in openmm and igb=8 in amber 18

[Iñigo Marcos Alcalde]

Hello all.

I am facing the following problem and I don't know if I am missing something really obvious or if it is an actual bug.

When simulating dsDNA with implicit solvent models, GBn and GBn2 exhibit unexpected behaviors unlike HCT, OBC1, OBC2 in openmm and igb=8 in amber 18.

Running implicit solvent simulations of dsDNA using GBn and GBn2 rapidly disrupts the interactions between both strands and even clusters the phosphate groups together. This behavior cannot be reproduced using HCT, OBC1, OBC2 nor running the simulation with amber 18 using igb=8 (the equivalent flag to GBn2).

This behavior was first detected on a run using the " linux-ppc64le/openmm-7.4.0-py37_cuda101" dev version with power9 processors and a Tesla V100 GPU, but has also been reproduced on and old 48 core cluster running the standard x64 release of openmm-7.4.0 (only the first 40 ps were computed given the limited performance, as it was enough to reproduce the behavior).

On contrast, the same system simulated via the pmemd.cuda_SPFP.MPI executable of an amber 18 ppc64le compilation, using a Tesla V100 and under igb=8 conditions, did not present this problem.

Here a tar.gz with the results and scripts used in these tests can be downloaded:
https://drive.google.com/open?id=1rdtb1 ... mJSLf-L1r-

[Peter Eastman]

That does sound strange. The first thing I would try is loading up a single conformation in both Amber and OpenMM and writing out the force on every atom. How well do they agree? If there are significant differences, the next thing would be to write out the GB parameters for every atom and see if any parameters are getting assigned differently in the two programs.

I'm not an Amber user, so I don't know exactly how to do those, but I assume there's a way? I can provide instructions on how to do them in OpenMM if that would help.

Note that Amber and OpenMM use different approximations for the surface area, so small differences in forces are expected. But the same approximation is used for all GB models, so that probably isn't the cause.

[Iñigo Marcos Alcalde]

I think I can try to do what you say. Printing forces atom-by-atom seems reasonably trivial in amber (http://archive.ambermd.org/201502/0174.html), and, as far as I know, gb radii are present in plain text in the prmtop file, so I should be able to extract those as well. Regarding OpenMM, it is pretty much new to me, so those instructions to extract forces and radii would be extremely helpful.

[Peter Eastman]

You can retrieve a list of the force on every atom like this:

Code: Select all

forces = simulation.context.getState(getForces=True).getForces()

GBn and GBn2 both get implemented with CustomGBForce. To look up the particle parameters, you first look up the force, then loop over all the particles and print out the parameters for each one.

Code: Select all

for force in system.getForces():
    if isinstance(force, CustomGBForce):
        for i in range(force.getNumParticles()):
            print(force.getParticleParameters(i))

The one thing that complicates this is that the parameters stored in the CustomGBForce have been preprocessed slightly to save work when computing the forces. The second element in the list of parameters is radius-0.009. The third element is screen*(radius-0.009). So to print out the original radius and screen values, change the above to

Code: Select all

params = force.getParticleParameters(i)
print(params[1]+0.009, params[2]/params[1])

One other important note. When using GBn or GBn2, it ignores the screen values given in the prmtop file and instead uses the values specified at https://github.com/openmm/openmm/blob/m ... #L392-L402. My understanding is that Amber does the same thing, but I don't know much about the details.

[Iñigo Marcos Alcalde]

I am afraid there might be a problem with the current parameter specification. In amber, nucleic acids have separate parameters of igb=8 (i.e. GBn2). In pages 62 and 63 of amber 2019 reference manual (http://ambermd.org/doc12/Amber19.pdf) the parameters indicated are (I have edited out the alpha, beta and gamma symbols so that it would allow me to post it here):

The following are the default parameters sander uses with igb=8:

Sh=1.425952, Sc=1.058554, Sn=0.733599,
So=1.061039, Ss=-0.703469, Sp=0.5,
offset=0.195141, gbneckscale=0.826836,
gbalphaH=0.788440, gbbetaH=0.798699, gbgammaH=0.437334,
gbalphaC=0.733756, gbbetaC=0.506378, gbgammaC=0.205844,
gbalphaN=0.503364, gbbetaN=0.316828, gbgammaN=0.192915,
gbalphaOS=0.867814, gbbetaOS=0.876635, gbgammaOS=0.387882,
gbalphaP=1.0, gbbetaP=0.8, gbgammaP=4.85 screen_hnu=1.69654,
screen_cnu=1.26890, screen_nnu=1.425974,
screen_onu=0.18401, screen_pnu=1.54506,
gb_alpha_hnu=0.53705, gb_beta_hnu=0.36286, gb_gamma_hnu=0.11670,
gb_alpha_cnu=0.33167, gb_beta_cnu=0.19684, gb_gamma_cnu=0.09342,
gb_alpha_nnu=0.68631, gb_beta_nnu=0.46319, gb_gamma_nnu=0.13872,
gb_alpha_onu=0.60634, gb_beta_onu=0.46301, gb_gamma_onu=0.14226,
gb_alpha_pnu=0.41836, gb_beta_pnu=0.29005, gb_gamma_pnu=0.10642

Parameters for proteins and for nucleic acids were optimized separately and can be independently specified. Protein parameters: Sh, Sc, Sn, So, Ss and Sp are scaling parameters, gbalphaX, gbbetaX, gbgammaX are the alpha, beta, gamma set for element X. gbalphaOS, gbbetaOS, gbgammaOS is the alpha, beta, gamma set applied to both O and S. The phosphorus parameters (in proteins) were not optimized and are simply taken as the parameters used in the OBC-2 model (igb=5). Nucleic acid parameters (end with "nu"): screen_Xnu (X=h, c, n, o, p) are scaling parameters, gb_alpha_Xnu (X=h, c, n, o, p) are the alpha, beta, gamma set for element X.

I have obtained the forces applied during the first step both in openmm and amber 18. To ensure the same initial positions for both, no energy minimization was performed prior MD this time in any of the cases. If I have made no mistakes during force extraction, after converting amber forces to (kilojoule/(nanometer*mole)))) and sorting by difference magnitude, it seems clear that the atoms participating in the phosphate groups present the highest differences. I attach both the sorted and unsorted force values in csv format.

Is there any way to force a GB parameter set for each atom? That way I could try to reproduce amber parameters and see if that solves the weird phosphate behavior.

[Peter Eastman]

Is there any way to force a GB parameter set for each atom?

Yes, it's very similar to the code above but call setParticleParameters() to give it the new parameters to use for each particle. You need to do that before you create the context/simulation. After you create it, it will ignore any further changes to the CustomGBForce.

We may want to consult the Amber developers about this. I'm not sure who would be most appropriate. This code was originally written by Jason Swails, but he mostly doesn't work on Amber anymore.

[Iñigo Marcos Alcalde]

Thank you very much for your help!

To be honest, I have no idea who we could ask. The amber reflector (amber@ambermd.org) is publicly available but it specifies that should be used for "scientific questions about installing or using Amber". On the other hand, the corresponding author of the paper used to reference the parameter optimization for nucleic acids in the amber manual is Carlos Simmerling:

Nguyen, H.; Perez, A.; Bermeo, S.; Simmerling, C. Refinement of Generalized Born Implicit Solvation Parameters for Nucleic Acids and Their Complexes with Proteins. J. Chem. Theory Comput., 2015, 11, 3714-3728.


In the meantime, I will try to set the parameters manually to fit those indicated in the amber manual and I'll report if I have any degree of success.

[Peter Eastman]

Let me see if I can find someone who can help us.

[Peter Eastman]

Ok, I think I know the problem. I checked with David Case and Carlos Simmerling, and it turns out Amber uses different parameters for nucleic acids than for proteins. The OpenMM code only looks at the element, not the type of residue, and always uses the protein parameters. Possibly it was based on an older version of Amber that didn't have the nucleic acid support? Anyway, it should be easy to fix. Given me a day or two.

[Iñigo Marcos Alcalde]

Thank you very much!