Details of analysis

[baba-hashimoto]

When I prepared the sample file myself and looked at the calculation results, I found that the following large numbers were detected.

Energy Average Err.Est. RMSD Tot-Drift

Potential -27534.7 6.7 786.594 20.1302 (kJ/mol)

Statistics over 2500001 steps [ 0.0000 through 5000.0000 ps ], 1 data sets
All statistics are over 2501 points (frames)

Energy Average Err.Est. RMSD Tot-Drift

Potential 4.57488e+09 2.7e+09 1.57499e+11 -1.40655e+10 (kJ/mol)

Statistics over 2500001 steps [ 0.0000 through 5000.0000 ps ], 1 data sets
All statistics are over 2501 points (frames)

What kind of command is outputting the above information?
I would appreciate your help.

[luancarvalhomartins]

Is this the output from the MD run, right? It seems to be GROMACS averages over the simulation. Which window is it? Is this large potential due to VdW or other term? Are you using REST? Did any of the MD crashed?

Potential energies could be large due to swapping replicas with different coulombic and VdW terms leading to a large potential in the immediately after the swapping. Also, the total potential energy could be very high in a few frames and this could affect the averages. Maybe the minimization and equilibration were not long enough. One thing you might consider is comparing to the potential energy during the equilibration (in case REST is being applied) and to the potential energy of the first and last windows.

I guess I would be more concerned about the ddG convergence, overlap matrices, forward and revers ddG, cycle closure errors, and other metrics related to the ddG estimation than to the potential energy during the MD.

[baba-hashimoto]

I'm so sorry for the late reply. I took your advice and did some research." I didn't know that the line containing "Potential" was the data from "gmx energy".

Fatal error:
Step 0: The total potential energy is 8.11303e+09, which is extremely high.
The LJ and electrostatic contributions to the energy are 8.11306e+09 and
-31622, respectively. A very high potential energy can be caused by
overlapping interactions in bonded interactions or very large coordinate
values. Usually this is caused by a badly- or non-equilibrated initial
configuration, incorrect interactions or parameters in the topology.

When I checked the rerun folder, I found that MD calculation was crashing in the coord0 files such as "rerun_struct5_coord0" and "rerun_struct6_coord0". All files except coord0 are working fine.

In such a case, where should I start recalculation? I would appreciate your comments.

[luancarvalhomartins]

I vaguely remember seeing this issue before, but IIRC my inputs were broken at that time. This is for step 0 and all rerun_structX_coord0.edr are problematic, but rerun_struct0_coord0.edr, which suggests there is an atom between an atom which is fully decoupled at lamba=0 and something else in the system. This is not something wrong per si, so I would not restart the calculation because of it. There should be a way to tell GROMACS to stop complaining and just go on, let me take a look at the gmx energy source and try figure out.

[baba-hashimoto]

Thanks for your reply. I have checked the details regarding the calculation and found the following.
It seems that numerical anomalies are output, so I added the "-pforce" command and checked the results.

Command line:
gmx_mpi mdrun -rerun lambda0/lambda.trr -s lambda5/lambda_nocont.tpr -e rerun/rerun_struct5_coord0.edr -g rerun/rerun_struct5_coord0.log -pforce 100000

starting md rerun 'Protein in water', reading coordinates from input trajectory 'lambda0/lambda.trr'

trr version: GMX_trn_file (single precision)
Reading frame 0 time 0.000 step 0 atom 31 x 2.983 0.364 1.359 force 4.20103e+06
step 0 atom 32 x 2.927 0.456 1.344 force 3.20596e+12
step 0 atom 43 x 0.242 2.651 1.477 force 1.80614e+08
step 0 atom 524 x 2.913 0.431 1.334 force 3.20596e+12
step 0 atom 884 x 0.272 2.591 1.519 force 1.80614e+08

---

Program: gmx mdrun, version 2020.6-MODIFIED
Source file: src/gromacs/mdlib/sim_util.cpp (line 465)

Fatal error:
Step 0: The total potential energy is 8.11303e+09, which is extremely high.
The LJ and electrostatic contributions to the energy are 8.11306e+09 and
-31622, respectively. A very high potential energy can be caused by
overlapping interactions in bonded interactions or very large coordinate
values. Usually this is caused by a badly- or non-equilibrated initial
configuration, incorrect interactions or parameters in the topology.

The data that worked fine before lambda5 is also attached.

Command line:
gmx mdrun -rerun lambda0/lambda.trr -s lambda0/lambda_nocont.tpr -e rerun/rerun_struct0_coord0.edr -g rerun/rerun_struct0_coord0.log

       Step           Time
          0        0.00000

Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih. LJ-14
5.36337e+01 1.12068e+02 6.13360e+01 9.92114e+00 4.92113e+01
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
-1.19071e+03 4.15771e+03 -2.43270e+02 -3.13429e+04 1.81366e+02
Potential
-2.81517e+04

Command line:
gmx mdrun -rerun lambda0/lambda.trr -s lambda1/lambda_nocont.tpr -e rerun/rerun_struct1_coord0.edr -g rerun/rerun_struct1_coord0.log

       Step           Time
          0        0.00000

Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih. LJ-14
5.36337e+01 1.12068e+02 6.13360e+01 9.92114e+00 4.91667e+01
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
-8.04359e+02 3.12002e+08 -2.43389e+02 -3.14089e+04 1.80636e+02
Potential
3.11970e+08

Command line:
gmx mdrun -rerun lambda0/lambda.trr -s lambda4/lambda_nocont.tpr -e rerun/rerun_struct4_coord0.edr -g rerun/rerun_struct4_coord0.log

Rerun does not report kinetic energy, total energy, temperature, virial and pressure.
Step Time
0 0.00000

Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih. LJ-14
5.36337e+01 1.12068e+02 6.13360e+01 9.92114e+00 4.90117e+01
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
0.00000e+00 2.68205e+09 -2.43828e+02 -3.16220e+04 1.83914e+02
Potential
2.68202e+09

I looked for a way to force it to work with "gmx mdrun -rerun", but I couldn't find one.
If you know of any information, I would appreciate your comments.

[luancarvalhomartins]

It seems gmx energy will always test for large potentials and die in case of finding it. However:

    /* We only check for large potential energy at the initial step,
     * because that is by far the most likely step for this too occur
     * and because computing the average kinetic energy is not free.
     * Note: nstcalcenergy >> 1 often does not allow to catch large energies
     * before they become NaN.
     */

So I guess maybe it's possible to trick gmx rerun by skipping the first frame. Would you mind sending me a example of a problematic .xtc/.trr and .tpr pair? It would be much simpler should I be able to test.

Thanks

[baba-hashimoto]

I will send e-mail(gmail) soon.

[luancarvalhomartins]

Thanks for your data. I think this is indeed an issue that PyAutoFEP should somehow address. Ideally, gmx mdrun -rerun should have a some flag to ignore this error and keep going, but it doesn't. I filled an issue (#30) for further discussing this.

[baba-hashimoto]

Thank you for your comment and support, I will consider changing the Window. Incidentally, if I change windows, will I have to write the input file from scratch?

[luancarvalhomartins]

Nopes, you only need to add lambda_input to the input. If you are referring to the lambda input, the answer is also no, but you can if you rather do so. You can (a) usa a custom lambda scheme by passing a text file or dict to lambda_input; (b) use one of the presets: lambdas12, lambdas24, lambdas16wang, and lambdas23wang; (c) mix several lambda schemes in the same perturbation map, passing a file to pertubation_map. See manual for further info. In the last case, the perturbation map file must be formatted as follows:

# Use the default lambdas12 for these pairs
LIG1 LIG2
LIG4 LIG5

# Use a custom value for this one
LIG3 LIG4 lambdas24

PS: just out of curiosity: how many heavy atoms are being perturbed? I am considering adding a auto selector for the number of lambdas windows, something like: more than X1 heavy atoms -> use Y2 lambda windows, more than X2 torions -> use Y2 lambda windows, charge hopping -> use Y3 lambda windows, and so on.

[baba-hashimoto]

Thank you for your comment. I have found the cause of the problem." How many heavy atoms are perturbed? I checked and found that 24/27 of LIG14 are perturbed atoms.

In lambda0 and lambda11, the ligand, with both vdw and charge set to 0, is free to move. In this case, there were only three atoms in the common site, resulting in a collision with the protein.　 I put the two picture

Snap shot after min03

Snap shot after npt

It is rotating around the center.

I think it would be useful if the number or percentage of perturbed atoms is calculated in the first step.

[luancarvalhomartins]

(I am assuming you are obfuscating the structures on purpose, which is OK)

I guess the rings differs in heteroatoms, so the MCS won't match. Currently, PyAutoFEP wont perturb bonds, so it would not be possible to obtain a physical representation for the common core, should heteroatoms be perturbed. I hope I will be able to add bond perturbing in a future release, but it is not something I am working on. BTW, I am currently working on charged perturbations. I hope to start working on covalent FEP after that, than, finally, I hope to put together the GUI.

Indeed, adding some warning if less than, say, 10 atoms are common or if more than 10 are being perturbed is a great idea. I am adding that do my to do.

[baba-hashimoto]

Thank you for your comment.
I think an automatic warning system would be beneficial. I look forward to improvements.