Oscillational period warning during tripeptide A2I forward simulation

pmx
1

Hello, I’m trying to perform forward and reverse simulations for Mutation free energy calculations. I have successfully calculated the total ddG of several tripeptides. They output similar ddG values to pmx webserver. However, there is some tripeptide that cannot be successfully simulated for either forward or reverse simulation. In this tripeptide A2I case, reverse simulation has succeeded, but the forward gave me these two warnings during the NPT equilibration:

WARNING 1 [file posre_tripep_A2I.itp, line 10]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B

WARNING 2 [file tripep_A2I.top, line 17]:
The bond in molecule-type Protein_chain_X between atoms 29 DCD1 and 30
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.

I’m using GROMACS v2022 and repeating the forward and reverse simulation 5 times for each protein. But all of the outputs of the forward simulations have the same error. I also get the same error for other several tripeptides. Is this related to the preprocessing or GROMACS version? I took the predetermined itp and topology files for the tripeptide from the pmx webserver.

Thank you

2

You seem to be using the older python2 based pmx version. There you will need to constrain all-bonds. If you want to constrain only hydrogen containing bonds, use the python3 based pmx version which you can find in the ‘develop’ branch. This will resolve the oscillation period warning.

3

Hello, thank for your reply. I do use python3 based pmx version from the ‘develop’ branch. I also have tried to reinstall the pmx and did the same protocol but still get the oscillation warning. Is it ok to just ignore the warning?

4

No, you shouldn’t ignore this warning.
The A2I mutation in the ‘amber99sb-star-ildn-mut.ff’ and ‘charmm36m-mut.ff’ forcefields in the ‘develop’ branch doesn’t have DHD1 atom, thus the warning that you see will not be triggered. Could you try ‘develop’ branch with one of these forcefields?

5

Thank you for your replies, it helps me a lot.

Now I’m performing double system-single box simulation for I2M mutation. The system consisted of the folded_wt (I2M) and unfolded_mutant (tripeptide M2I) separated at least 3 nm which solvated and neutralized in a dodecahedron box.
I generated the ACE-GLY-MET-GLY-NME capped tripeptide using tleap program and amber99sb-star-ildn force field.
I do use the python 3 pmx ‘develop’ branch to mutate both the folded_wt and unfolded_mutant pairs. However, I always face these warnings and notes during nvt, npt, and equilibirum simulations:

WARNING 1 [~/dssb_mdp/eqA/nvt-A.mdp, line 62]:
You are using soft-core interactions while the Van der Waals interactions
are not decoupled (note that the sc-coul option is only active when using
lambda states). Although this will not lead to errors, you will need much
more sampling than without soft-core interactions. Consider using
sc-alpha=0.

NOTE 1 [file ~/dssb_mdp/eqA/nvt-A.mdp]:
With PME there is a minor soft core effect present at the cut-off,
proportional to (LJsigma/rcoulomb)^6. This could have a minor effect on
energy conservation, but usually other effects dominate. With a common
sigma value of 0.34 nm the fraction of the particle-particle potential at
the cut-off at lambda=0.5 is around 3.9e-05, while ewald-rtol is 1.0e-05.

Generating 1-4 interactions: fudge = 0.5

WARNING 2 [file posre_Protein_chain_A.itp, line 233]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B

NOTE 2 [file newtop.top, line 25]:
The bond in molecule-type Protein_chain_X between atoms 18 CB and 31 DCG2
has an estimated oscillational period of 1.9e-02 ps, which is less than
10 times the time step of 2.0e-03 ps.
Maybe you forgot to change the constraints mdp option.

Number of degrees of freedom in T-Coupling group Protein is 6647.00
Number of degrees of freedom in T-Coupling group non-Protein is 126471.00

NOTE 3 [file ~/dssb_mdp/eqA/nvt-A.mdp]:
Removing center of mass motion in the presence of position restraints
might cause artifacts. When you are using position restraints to
equilibrate a macro-molecule, the artifacts are usually negligible.

There were 3 notes

There were 2 warnings

I think Warning_1 and Note_1 can be ignored as described by the book. But how about the Warning_2 and (especially) Note_2?

  • Warning_2, posre_Protein_chain_A.itp, line 233:

453 1 1000 1000 1000

  • Note 2: I think DCG2 atom for M2I exists in ‘develop’ branch and not in ‘master’ but why I still get oscillational period messages?

And also I have a question about the soft-core function. Most of the existing examples are using the sc parameters for sc-function=beutler. As now the GROMACS 2022 version supports sc-function=gapsys. Is it better to use this function instead?

6

Yes, I see now what is causing it. The default value for “–scale_mass” flag is set too low (currently it is 0.25). This flag tells to scale down dummy atom masses by multiplying the mass of the real atom in another state by the scaling factor. Such scaling helps with simulation stability for the perturbations where many dummies need to be introduced. However, scaling also cannot be too low, because it will trigger the Notes about too frequent vibration period that you see now. Ideally, I would need to check for the oscillation period within pmx, but as a quick fix for now I increased this scaling factor to 0.33. I pushed this change to the ‘develop’ branch.

Regarding sc-function, in many cases both functions will yield identical free energy estimates, so for now you could as well continue with the standard beutler variant. However, if you notice some erratic behavior of the calculated free energies (due to large jumps in the dhdl curves), I would recommend trying out another function with the default parameters. Such problematic cases usually occur for more complex perturbations, e.g. mutations involving charge change.

1 Like
7

Hi @vgapsys,
Should it be the same case of ligand hybrids where I let the tool to do the thing(default parameters) and do not explicitly provide a number to scale dummy masses or angle parameters or dihedral parameters.
Regards,
Pallav.

8

If you observe instabilities, e.g. simulations show LINCS warnings and crash, scaling dummy masses might help. Otherwise, you can continue with the default settings of non-scaled masses.

1 Like
9

Thanks a lot! now it worked smoothly without the oscillational period warning. The soft-core warning has gone because I changed the sc-function to ‘gapsys’. However, I still got the:

WARNING 1 [file posre_Protein_chain_A.itp, line 233]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B

Is it totally fine to ignore it?

10

What interaction is defined on line 233?

11

file posre_Protein_chain_A.itp, line 233:

453 1 1000 1000 1000

atom #453 in pmx_topol_Protein_chain_A.itp:

450 O 26 THR O 450 -0.567900 16.0000
451 N 27 I2M N 451 -0.415700 14.0100
452 H 27 I2M H 452 0.271900 1.0080
453 CT 27 I2M CA 453 -0.059700 12.0100 CT -0.023700 12.0100
454 H1 27 I2M HA 454 0.086900 1.0080 H1 0.088000 1.0080
455 CT 27 I2M CB 455 0.130300 12.0100 CT 0.034200 12.0100
456 HC 27 I2M HB 456 0.018700 1.0080 HC 0.024100 1.0080
457 CT 27 I2M CG1 457 -0.043000 12.0100 CT 0.001800 12.0100
458 HC 27 I2M HG11 458 0.023600 1.0080 H1 0.044000 1.0080
459 HC 27 I2M HG12 459 0.023600 1.0080 H1 0.044000 1.0080
460 CT 27 I2M CG2 460 -0.320400 12.0100 DUM_CT 0.000000 3.9633
461 HC 27 I2M HG21 461 0.088200 1.0080 DUM_HC 0.000000 1.0000
462 HC 27 I2M HG22 462 0.088200 1.0080 DUM_HC 0.000000 1.0000
463 HC 27 I2M HG23 463 0.088200 1.0080 DUM_HC 0.000000 1.0000
464 CT 27 I2M CD1 464 -0.066000 12.0100 S -0.273700 32.0600
465 HC 27 I2M HD11 465 0.018600 1.0080 DUM_HC 0.000000 1.0000
466 HC 27 I2M HD12 466 0.018600 1.0080 DUM_HC 0.000000 1.0000
467 HC 27 I2M HD13 467 0.018600 1.0080 DUM_HC 0.000000 1.0000
468 C 27 I2M C 468 0.597300 12.0100
469 O 27 I2M O 469 -0.567900 16.0000
470 DUM_HC 27 I2M HV1 470 0.000000 1.0000 HC 0.024100 1.0080
471 DUM_CT 27 I2M DCE 471 0.000000 3.9633 CT -0.053600 12.0100
472 DUM_H1 27 I2M HV2 472 0.000000 1.0000 H1 0.068400 1.0080
473 DUM_H1 27 I2M HV3 473 0.000000 1.0000 H1 0.068400 1.0080
474 DUM_H1 27 I2M HV4 474 0.000000 1.0000 H1 0.068400 1.0080

12

Right, so the perturbed Calpha atom has a position restraint. The posre is defined for stateA only, but since the atom is perturbed it is not clear what position restraint should be in the stateB. In this case Gromacs warns that it will simply use the same force constants for posre in both states. I think this is what you would like to have, so it is safe to ignore this warning.

1 Like
13

I noticed that some of the mutation types need larger values for the --scale_mass (i.e > 0.6) to avoid the oscillation warning. However, it seems that larger mass-scale values often cause energy minimization failure. On the other hand, when low/default value (0.33) is outputting an oscillation warning, the MD steps went smoothly through the production step. Then, when I use the default scale and ignore the oscillation warning, the ddG results from double-system single-box simulation on 5 different mutation cases (3 simulations each → 15 simulations in total) show a relatively accurate to experimental value with around <= 1 kcal/mol difference.

Is there any way to check the optimal scale value for each mutation type?
In my case, the results show somehow large scale value causes the EM step to fail and the oscillation warning from too low scale-mass value does not really affect the calculated ddG performance.

14

Masses should have no influence on energy minimization. Are you sure there is nothing else at play there? Maybe it could be related to some energy minimization bug, e.g. Energy minimisation won't go to the minimum point when GPU is used. (#4533) · Issues · GROMACS / GROMACS · GitLab

You are mentioning a warning about oscillation period with the scaling of 0.33: are you sure this is a warning and not a note? If it is a note, it is safe to ignore it

15

Yes, you are correct. Sorry it was my misunderstanding. Thanks a lot! The link you gave solved my EM problem. I tried to find solutions for the energy minimization failure problems previously, but only found generic solutions.

Since it was a Warning previous GROMACS version (before 2022), I thought it was harmful even if it became a note instead the latest version. I see, so it safe to ignore it. Many many thanks for your kind replies and solutions!

1 Like