[gmx-users] Inconsistent behavior for bonds involving dummy atoms for alchemical perturbation
Justin Lemkul
jalemkul at vt.edu
Fri Mar 18 12:57:41 CET 2016
On 3/17/16 6:04 PM, Ryan Muraglia wrote:
> Dear GROMACS users,
>
> I've been attempting to calculate delta G values for GXG tripeptide
> mutations as a test case to validate and benchmark my set up, but I have
> run into a host of issues and curiosities, one of which I describe below.
>
> Depending on the direction in which I carry out a mutation, a given bond
> may or may not throw a (fatal) warning as I execute the grompp call to set
> up the NVT equilibration (and subsequently mdruns).
>
> To clarify, here is an example: when running a GAG to GLG mutation, during
> the grompp call I receive the following warning:
>
> "WARNING 1 [file topol.top, line 388]:
> The bond in molecule-type Protein between atoms 21 DCD1 and 22 DHD1 has
> an estimated oscillational period of 8.3e-03 ps, which is less than 5
> times the time step of 2.0e-03 ps.
> Maybe you forgot to change the constraints mdp option."
>
> If I choose to ignore this warning by setting maxwarn to 1, the NVT
> equilibration blows up.
>
> In the other direction, when running a GLG to GAG mutation, for the same
> mdp file, the grompp step before the NVT throws no warning, and the NVT
> equilibration runs to completion without a hitch.
>
> I can reduce the time step for the GAG to GLG mutation to have it run to
> completion, but I do not understand why I only see this warning in one
> direction, when this bond exists in both topologies. I also do not
> understand why there is only one bond for which this warning comes up, when
> there are several other identical bonds (in terms of atom types, etc) that
> fly under the radar. Lastly, I have set "constraints = h-bonds" so it
> puzzles me that a bond involving a hydrogen is the culprit here.
>
If the atoms are named "D" then grompp won't detect them as hydrogen. The
checking is done with respect to the first character in the atom name. So
that's why you're blowing up; a bond that should be constrained isn't.
> I use the pmx package to generate topologies for amino acid substitutions,
> and my system uses the amber99sb forcefield with tip3p water. The [atoms]
> directive of my A2L topology is reproduced below:
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass typeB
> chargeB massB
> 1 N3 1 GLY N 1 0.294300 14.0100
> 2 H 1 GLY H1 2 0.164200 1.0080
> 3 H 1 GLY H2 3 0.164200 1.0080
> 4 H 1 GLY H3 4 0.164200 1.0080
> 5 CT 1 GLY CA 5 -0.010000 12.0100
> 6 HP 1 GLY HA1 6 0.089500 1.0080
> 7 HP 1 GLY HA2 7 0.089500 1.0080
> 8 C 1 GLY C 8 0.616300 12.0100
> 9 O 1 GLY O 9 -0.572200 16.0000
> 10 N 2 A2L N 10 -0.415700 14.0100
> 11 H 2 A2L H 11 0.271900 1.0080
> 12 CT 2 A2L CA 12 0.033700 12.0100
> CT -0.051800 12.0100
> 13 H1 2 A2L HA 13 0.082300 1.0080
> H1 0.092200 1.0080
> 14 CT 2 A2L CB 14 -0.182500 12.0100
> CT -0.110200 12.0100
> 15 HC 2 A2L HB1 15 0.060300 1.0080
> CT 0.353100 12.0100
> 16 HC 2 A2L HB2 16 0.060300 1.0080
> HC 0.045700 1.0080
> 17 HC 2 A2L HB3 17 0.060300 1.0080
> HC 0.045700 1.0080
> 18 C 2 A2L C 18 0.597300 12.0100
> 19 O 2 A2L O 19 -0.567900 16.0000
> 20 DUM_HC 2 A2L DHG 20 0.000000 1.0000
> HC -0.036100 1.0080
> 21 DUM_CT 2 A2L DCD1 21 0.000000 1.0000
> CT -0.412100 12.0100
> 22 DUM_HC 2 A2L DHD1 22 0.000000 1.0000
> HC 0.100000 1.0080
> 23 DUM_HC 2 A2L DHD2 23 0.000000 1.0000
> HC 0.100000 1.0080
> 24 DUM_HC 2 A2L DHD3 24 0.000000 1.0000
> HC 0.100000 1.0080
> 25 DUM_CT 2 A2L DCD2 25 0.000000 1.0000
> CT -0.412100 12.0100
> 26 DUM_HC 2 A2L DHD4 26 0.000000 1.0000
> HC 0.100000 1.0080
> 27 DUM_HC 2 A2L DHD5 27 0.000000 1.0000
> HC 0.100000 1.0080
> 28 DUM_HC 2 A2L DHD6 28 0.000000 1.0000
> HC 0.100000 1.0080
> 29 N 3 GLY N 29 -0.382100 14.0100
> 30 H 3 GLY H 30 0.268100 1.0080
> 31 CT 3 GLY CA 31 -0.249300 12.0100
> 32 H1 3 GLY HA1 32 0.105600 1.0080
> 33 H1 3 GLY HA2 33 0.105600 1.0080
> 34 C 3 GLY C 34 0.723100 12.0100
> 35 O2 3 GLY OC1 35 -0.785500 16.0000
> 36 O2 3 GLY OC2 36 -0.785500 16.0000
>
> The L2A topology is similar. The second residue is as follows:
>
> 10 N 2 L2A N 10 -0.415700 14.0100
> 11 H 2 L2A H 11 0.271900 1.0080
> 12 CT 2 L2A CA 12 -0.051800 12.0100
> CT 0.033700 12.0100
> 13 H1 2 L2A HA 13 0.092200 1.0080
> H1 0.082300 1.0080
> 14 CT 2 L2A CB 14 -0.110200 12.0100
> CT -0.182500 12.0100
> 15 HC 2 L2A HB1 15 0.045700 1.0080
> HC 0.060300 1.0080
> 16 HC 2 L2A HB2 16 0.045700 1.0080
> HC 0.060300 1.0080
> 17 CT 2 L2A CG 17 0.353100 12.0100
> HC 0.060300 1.0080
> 18 HC 2 L2A HG 18 -0.036100 1.0080
> DUM_HC 0.000000 1.0000
> 19 CT 2 L2A CD1 19 -0.412100 12.0100
> DUM_CT 0.000000 1.0000
> 20 HC 2 L2A HD11 20 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 21 HC 2 L2A HD12 21 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 22 HC 2 L2A HD13 22 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 23 CT 2 L2A CD2 23 -0.412100 12.0100
> DUM_CT 0.000000 1.0000
> 24 HC 2 L2A HD21 24 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 25 HC 2 L2A HD22 25 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 26 HC 2 L2A HD23 26 0.100000 1.0080
> DUM_HC 0.000000 1.0000
> 27 C 2 L2A C 27 0.597300 12.0100
> 28 O 2 L2A O 28 -0.567900 16.0000
>
> An additional quirk is that in the A2L direction (the one that is
> problematic), mdrun reports that I have 16 constraints in my system,
> whereas in the L2A direction, mdrun reports 19 constraints.
>
Use gmx dump or gmx check on the .tpr files to see what the differences are.
-Justin
> Other mutations are a mixed bag. A2V and V2A, despite having the same type
> of bond that causes the error in A2L, both go through with no warnings
> (although their constraints number do differ, but they follow the pattern
> described in the pmx paper). Q2A goes through with no warning, but A2Q
> prints a similar warning as above, except this one is between atoms 22 DCD
> and 23 DOE1 with an estimated oscillational period of 6.4e-03 ps,
> indicating that the problem can arise for bonds that do not involve
> hydrogen as well.
>
> Pertinent portions of my mdp file are reproduced below:
>
> integrator = sd ; Langevin dynamics
> dt = 0.002 ; 2 fs
> constraints = h-bonds
> constraint-algorithm = LINCS
> continuation = no
> lincs-order = 6
> lincs-iter = 1
>
> ; Free energy calculations
> free_energy = yes
> delta_lambda = 0 ; no Jarzynski non-eq
> calc_lambda_neighbors = 1 ; only calculate energy to immediate neighbors
> (suitable for BAR, but MBAR needs all)
> sc-alpha = 0.5
> sc-coul = no
> sc-power = 1.0
> sc-sigma = 0.3
> couple-moltype = Protein ; name of moleculetype to decouple
> couple-lambda0 = vdw-q ; all interactions
> couple-lambda1 = vdw ; remove electrostatics, only vdW
> couple-intramol = no
> nstdhdl = 100
>
> ; lambda vectors ; decharging only. Keep as A. Next file deals with
> uncharged A to uncharged L to charged L
> ; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10
> init_lambda_state = 05
> coul_lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
> vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
> bonded_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; match
> vdw
> mass_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; match
> vdw
> temperature_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ; not
> doing simulated tempering
>
> I considered that perhaps the issue was due to which atoms had "appeared"
> or "disappeared," so I also ran a test using an mdp file where the vdw,
> bonded and mass lambdas were set to 1 for the L2A mutation, which I believe
> means I should be simulating the same state as the A2L which causes an
> error, but I got no warning from grompp.
>
> Any information or advice as to why I am seeing unexpectedly different
> behaviors depending on the direction of my mutation would be appreciated.
>
> Thank you!
>
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list