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 rmurag...@gmail.com -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.