Discovered a critical bug in gmx_mmpbsa! Inconsistent VDW and EEL Energies between FINAL_DECOMP_MMPBSA.csv and FINAL_RESULTS_MMPBSA.csv in gmx_mmpbsa

[Bingxu Qi]

I have discovered a potential issue in the energy decomposition of gmx_MMPBSA. When calculating the van der Waals (vdW) and electrostatic (EEL) energy terms, the results in FINAL_DECOMP_MMPBSA.csv and FINAL_RESULTS_MMPBSA.csv are inconsistent.

For example, for the overall decomposition, the reported values are:

vdW = 2A EEL = B

However, using the pairwise decomposition method, the values become:

vdW = 2A (consistent) EEL = 2B

Taking frame 2500 as an example: the pairwise method shows two entries for the interaction between R:A:ARG:68 and L:C:TYR:8, each with vdW = -0.63 and EEL = -7.22. This yields a total of -1.26 for vdW (which matches the overall value) but -14.44 for EEL, while FINAL_RESULTS_MMPBSA.csv reports only -7.22 for EEL.

Since the receptor and ligand in my system are represented by a single amino acid each, the overall energy should equal the sum of the pairwise contributions. This discrepancy in the EEL term suggests there might be a double-counting or another calculation error in the pairwise method.

Could this be a bug, or is there an explanation for this difference? Any insights or suggestions would be greatly appreciated.

Terminal output

For FINAL_DECOMP_MMPBSA.csv DELTAS:
Total Decomposition Contribution (TDC)
Frame # Resid 1 Resid 2 Internal van der Waals Electrostatic Polar Solvation Non-Polar Solv.
TOTAL DELTAS: Total Decomposition Contribution (TDC)
Frame # Resid 1 Resid 2 Internal van der Waals Electrostatic Polar Solvation Non-Polar Solv. TOTAL
2500 R:A:ARG:68 R:A:ARG:68 0 0 0 1.48 0 1.48
2500 R:A:ARG:68 L:C:TYR:8 0 -0.63 -7.22 1.78 0 -6.07
2500 L:C:TYR:8 R:A:ARG:68 0 -0.63 -7.22 1.78 0 -6.07
2500 L:C:TYR:8 L:C:TYR:8 0 0 0 1.07 0 1.07
2501 R:A:ARG:68 R:A:ARG:68 0 0 0 1.29 0 1.29
2501 R:A:ARG:68 L:C:TYR:8 0 -0.96 -6.3 1.89 0 -5.37
2501 L:C:TYR:8 R:A:ARG:68 0 -0.96 -6.3 1.89 0 -5.37
2501 L:C:TYR:8 L:C:TYR:8 0 0 0 0.94 0 0.94
2502 R:A:ARG:68 R:A:ARG:68 0 0 0 1.47 0 1.47
2502 R:A:ARG:68 L:C:TYR:8 0 -0.7 -6.72 1.79 0 -5.63
2502 L:C:TYR:8 R:A:ARG:68 0 -0.7 -6.72 1.79 0 -5.63
2502 L:C:TYR:8 L:C:TYR:8 0 0 0 0.91 0 0.91
For FINAL_RESULTS_MMPBSA.csv
Delta Energy Terms
Frame # BOND ANGLE DIHED VDWAALS EEL 1-4 VDW 1-4 EEL EPB ENPOLAR EDISPER GGAS GSOLV TOTAL
2500 0 0 0 -1.26 -7.22 0 0 37.74 -0.57 0 -8.48 37.17 28.69
2501 0 0 0 -1.92 -6.3 0 0 37.9 -0.51 0 -8.22 37.39 29.17
2502 0 0 0 -1.39 -6.72 0 0 37.72 -0.55 0 -8.11 37.17 29.06

[mariosergi...@gmail.com]

Sorry for the delay. We've been swamped. We'll look into this carefully, although the problem would come from Sander directly if it were the case. The per-wise formula considers interaction duplication, but it could be a problem we haven't noticed. In any case, we'll be looking into the issue in the next few days.
Thanks for reporting it and for using our tool.

All the best

Mario S.

[Bingxu Qi]

Dear Mario,

Thank you for your prompt response and for looking into this issue. I really appreciate the effort your team is putting into investigating it.

If needed, I’d be happy to provide my test files to assist with debugging. Please let me know if that would be helpful.

From my perspective, when using GMX_MMPBSA, the consistency of the van der Waals (vdW) term—where vdW = 2A—makes sense. However, for electrostatics (EEL), I personally find that EEL = 2B would be more logical than EEL = B, maintaining the same approach as vdW. If there's an underlying reason for the current implementation, I’d love to understand it better.

I’d greatly appreciate any updates on this matter and look forward to your response at your earliest convenience. Thanks again for your work on this great tool!

Best regards,

Bingxu Qi

State Key Laboratory of Structural Chemistry
Fujian Institute of Research on the Structure of Matter, Chinese Academy of Sciences
Fuzhou, CHINA

[Bingxu Qi]

Dear Mario,

I hope this email finds you well. I truly appreciate the effort your team is dedicating to investigating the potential inconsistency between the FINAL_DECOMP_MMPBSA.csv and FINAL_RESULTS_MMPBSA.csv outputs.

However, since a few days have passed without further clarification, I kindly inquire whether any conclusions have been reached regarding which approach is correct. Specifically, I am interested in knowing if the pairwise method's treatment of the electrostatic (EEL) term—resulting in a doubled value—is intentional or indicative of a calculation error. Alternatively, if you recommend a temporary approach I should adopt for my current research, I would greatly appreciate your guidance.

Best regards,

Bingxu Qi
State Key Laboratory of Structural Chemistry
Fujian Institute of Research on the Structure of Matter, Chinese Academy of Sciences
Fuzhou, CHINA

On Wednesday, February 19, 2025 at 11:34:51 AM UTC+8 mariosergi...@gmail.com wrote:

[marioe911116]

HI there! Can you run a calculation with just 5 frames and send a folder with all the files you are using and obtaining after running gmx_MMPBSA?

[Bingxu Qi]

Thank you for your message. I am happy to run a calculation using just 5 frames as requested. I will compile a folder containing all the files used and generated during the gmx_MMPBSA run, and send it to you as follow.

I appreciate your continued support and look forward to your feedback.

[marioe911116]

Hi there! I received the files, thanks... Could you please, using these new files, pinpoint where the issue is exactly?

thanks in advance

[Bingxu Qi]

I’m happy to explain my files. In my analysis, “VDW” refers to van der Waals interactions and “EEL” to electrostatic interactions. When using the decomposition method, I obtain the files FINAL_DECOMP_MMPBSA.csv (and .dat) in addition to the general FINAL_RESULTS_MMPBSA.csv (and .dat) files, which are produced regardless of the method used.

Focusing on the “DELTAS:” section in both CSV files—for example, for frame 2501, FINAL_DECOMP_MMPBSA.csv shows:

Frame #    Resid 1      Resid 2      Internal   van der Waals   Electrostatic   Polar Solvation   Non-Polar Solv.   TOTAL

2501       R:A:ARG:68  R:A:ARG:68  0          0               0               1.26              0                 1.26
2501       R:A:ARG:68  L:C:TYR:8  0          -0.96           -6.3            1.9               0                 -5.36
2501       L:C:TYR:8  R:A:ARG:68  0          -0.96           -6.3            1.9               0                 -5.36
2501       L:C:TYR:8  L:C:TYR:8  0          0               0               0.93              0                 0.93

Here, the van der Waals energy between R:A:ARG:68 and L:C:TYR:8 is –0.96 (from one direction) plus –0.96 (the reciprocal), giving –1.92, and the electrostatic term is –6.3 – 6.3 = –12.6.

In contrast, FINAL_RESULTS_MMPBSA.csv for frame 2501 reports:

Frame #   BOND   ANGLE   DIHED   VDWAALS   EEL   1-4 VDW   1-4 EEL   EPB   ENPOLAR   EDISPER   GGAS   GSOLV   TOTAL

2501      0      0       0       -1.92     -6.3  0         0         37.79 -0.52     0         -8.22  37.27   29.05

Thus, while the VDW contributions are identical (–1.92), the EEL value in the decomposition file (–12.6) is exactly twice that reported in FINAL_RESULTS_MMPBSA.csv (–6.3).

Looking at the .dat files confirms this discrepancy. In FINAL_DECOMP_MMPBSA.dat, for the R:A:ARG:68–L:C:TYR:8 interaction, each instance of the VDW contribution is –0.4394632 (summing to –0.88) and each EEL contribution is –6.9312608 (summing to –13.86). However, FINAL_RESULTS_MMPBSA.dat reports ΔVDWAALS as –0.88 and ΔEEL as –6.93.

In summary, regardless of the method, when calculating the binding energy between R:A:ARG:68 and L:C:TYR:8, the per-frame VDW and EEL energies should match exactly. Instead, while the VDW values are consistent, the EEL energy in the decomposition files appears doubled. This means that in FINAL_RESULTS_MMPBSA files, VDW = 2A (consistent) and EEL = B, whereas in FINAL_DECOMP_MMPBSA files, VDW = 2A and EEL = 2B. This factor-of-two discrepancy in the electrostatic term is unexpected and is the issue we need to pinpoint.

[Bingxu Qi]

I am confident that your development team is highly professional and possesses an in‐depth understanding of binding energy calculations. You should be fully aware of the significant discrepancy present here. Since we use the gmx_MMPBSA tool for our research, we do not merely rely on the overall conclusions presented in the FINAL_RESULTS_MMPBSA.csv (and .dat) files. Instead, we frequently need to analyze in detail—within the FINAL_DECOMP_MMPBSA.csv (and .dat) files—which energy components contribute to the total energies.

However, due to the fact that

“This means that in FINAL_RESULTS_MMPBSA files, VDW = 2A (consistent) and EEL = B, whereas in FINAL_DECOMP_MMPBSA files, VDW = 2A and EEL = 2B. This factor-of-two discrepancy in the electrostatic term is unexpected and is the issue we need to pinpoint.”

the energy values analyzed in FINAL_DECOMP_MMPBSA.csv (and .dat) do not correspond to those in FINAL_RESULTS_MMPBSA.csv (and .dat), which is quite astonishing. In other binding energy calculation tools I have used, decomposition energies have consistently matched the overall results. This issue is not specific to my system but arises in any protein–small molecule binding scenario.

I encourage you to conduct further tests. With your extensive experience as chemists, it should not be difficult to identify the key differences in these energy components. Thank you for addressing my inquiry and for your efforts—I look forward to your prompt response.

On Monday, February 24, 2025 at 10:15:16 AM UTC+8 marioe911116 wrote:

[mariosergi...@gmail.com]

We apologize for the delay in solving this problem. We are making an effort to provide a solution over other commitments. Thanks for the report and for taking the time to explain. Please send me the top, tpr, ndx, and 5 frames xtc files so I can test some config options and check what is happening.

[Bingxu Qi]

I believe I have explained my concerns very clearly, yet you still have not provided any helpful explanation regarding my questions. In the gmx_MMPBSA tool your team developed, there is content of this nature. Please refer to the following URL:
https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA/tree/master/examples/Decomposition_analysis
This directory includes the topology (.top), tpr, ndx, and other files.

[marioe911116]

Hello there! 

I believe I have explained my concerns very clearly,

Yes, you did! and we appreciate that.

yet you still have not provided any helpful explanation regarding my questions.

In order to identify and solve the issue, we MUST understand what the issue is about first; otherwise, we cannot do our work properly. If you look carefully at other conversations of ours in this very group, you will notice that we usually message back and forth with the user in order to have a deep understanding of the issue first, and when it is an issue like the one you just identified, we usually request politely the files to reproduce the issue by ourselves. Sometimes issues arise when working with real-life examples and not the models used in the tutorials (it's understandable that we are unable to cover every single aspect with the models in the tutorials).

In the gmx_MMPBSA tool your team developed, there is content of this nature. Please refer to the following URL:
https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA/tree/master/examples/Decomposition_analysis
This directory includes the topology (.top), tpr, ndx, and other files.

We are aware of this; after all, we were the ones who designed the tutorials. However, and I want to make myself very clear here, you DID NOT identify the issue working with the files we provided in the tutorials; you did it using your files, hence the request from our part of the files you were using.

Given the circumstances, we see three clear paths forward:

1- Collaborate with us by providing the necessary files so we can troubleshoot efficiently. Please note that gmx_MMPBSA is maintained voluntarily and for free, and our time is limited.
2- Investigate and resolve the issue yourself, and, if you find a fix, contribute back via a pull request, as many others have done.
3- Consider alternative software if gmx_MMPBSA does not meet your requirements.

We hope to work together towards a solution. Let us know how you’d like to proceed.

Sincerely!

gmx_MMPBSA team

[Bingxu Qi]

Dear Team,

Thank you very much for your patient and detailed responses. I now understand which files need to be provided. I will be sending you the file named Decomposition_analysis-by_qbx.zip. This ZIP archive includes the topology (.top), tpr, ndx, and other files.

It is important to note that:

1. The topology (.top), tpr, ndx, and other files all originate from the following URL:
   https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA/tree/master/examples/Decomposition_analysis
2. I have only made some modifications to the mmpbsa-all_10frames.in file.

Despite these adjustments, the discrepancy still persists: in the FINAL_RESULTS_MMPBSA files, VDW = 2A (consistent) and EEL = B, whereas in the FINAL_DECOMP_MMPBSA files, VDW = 2A and EEL = 2B.

Thank you once again for your team's effort and rigorous approach. I look forward to your guidance on resolving this issue.

Sincerely,

Bingxu Qi

[marioe911116]

We have an idea now of what the issue might be... You are using dielectric constant = 2. The dielectric constant of 2 is being taken into account for the calculations reported in FINAL_RESULTS_MMPBSA.csv but not in FINAL_DECOMP_MMPBSA.csv, hence the discrepancy. This is an "issue" coming from Sander, although we could correct for that in the output file now that we know it happens.

We'll open an issue with this and fix it for the next version. We gotta perform some tests to see that everything works after we include the correction, and that could take some time.

Please let us know if you have any questions.

[Bingxu Qi]

Dear Team,

Thank you for your detailed explanation—it has been very helpful. I ran some tests on my end and confirmed that when the dielectric constant is set to 1, the VDW and EEL values in both FINAL_RESULTS_MMPBSA.csv and FINAL_DECOMP_MMPBSA.csv match perfectly. Since both files are generated using the same decomposition method, it makes complete sense that their results should be consistent.

I truly appreciate your prompt and thorough response.

Best regards,

Your BingxuQi

[Bingxu Qi]

Dear  Team,

Thank you for your previous response. I have some other doubts regarding the contents of the FINAL_RESULTS_MMPBSA.csv and FINAL_DECOMP_MMPBSA.csv files included in the Decomposition_analysis-by_qbx.zip I sent earlier.

In the FINAL_DECOMP_MMPBSA.csv file, the “Non-Polar Solv.” value is consistently 0, whereas in the FINAL_RESULTS_MMPBSA.csv file, the ENPOLAR value is non-zero. According to the information provided at https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/output/, it is mentioned that “For GSOLV, the polar and non-polar contributions are EGB (or EPB) and ESURF (or ENPOLAR + EDISPER), respectively for GB (or PB) calculations.”

Thus, one would expect that the “Non-Polar Solv.” in the FINAL_DECOMP_MMPBSA.csv file should be equal to the sum of ENPOLAR and EDISPER from the FINAL_RESULTS_MMPBSA.csv file. However, the values are clearly different. Could you please explain the reason for this discrepancy?

Thank you very much for your time and assistance.

Sincerely,
Bingxu Qi

On Wednesday, February 26, 2025 at 2:25:53 PM UTC+8 marioe911116 wrote:

[marioe911116]

According to this thread in the Amber official mainling list, it seems PB non-polar can't be decomposed, only GB...

http://archive.ambermd.org/201307/0733.html

[Bingxu Qi]

Thank you very much for your clarification. I have reviewed the thread on the Amber official mailing list (http://archive.ambermd.org/201307/0733.html) and now understand that PB non-polar contributions cannot be decomposed, whereas this decomposition is possible for GB calculations. I appreciate your detailed explanation and prompt assistance in resolving my doubts.