[gmx-users] Inconsistent behavior for bonds involving dummy atoms for alchemical perturbation
Ryan Muraglia
rmuraglia at gmail.com
Thu Mar 17 23:04:52 CET 2016
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.
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.
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!
--
Ryan Muraglia
rmuraglia at gmail.com
More information about the gromacs.org_gmx-users
mailing list