Valdes-Tresanco-MS / gmx_MMPBSA

gmx_MMPBSA is a new tool based on AMBER's MMPBSA.py aiming to perform end-state free energy calculations with GROMACS files.
https://valdes-tresanco-ms.github.io/gmx_MMPBSA/
GNU General Public License v3.0
229 stars 66 forks source link

[Bug-gmx_MMPBSA]: #528

Closed JuuelHerza closed 3 months ago

JuuelHerza commented 3 months ago

Bug summary

For this project I have used the program for around ~8 compounds, three replicates for each compound, everything is fine except for one compound. But I don't know why it is impossible to get a calculation for one ligand. The programs just gets stuck at "[INFO ] calculating ligand contribution..." No error messages to follow. I checked the trajectory it is fine, without jumps, PBC already removed, and I cannot get a number out of any of the three trajectories involving such compound.

Terminal output

[INFO   ] Running calculations on normal system...
[INFO   ] Beginning PB calculations with /home/juuelherza/anaconda3/envs/gmxMMPBSA/bin/sander
[INFO   ]   calculating complex contribution...
              0%|          | 0/100 [elapsed: 00:00 remaining: ?              6%|6         | 6/100 [elapsed: 00:36 remaining: 0             12%|#2        | 12/100 [elapsed: 00:44 remaining:              20%|##        | 20/100 [elapsed: 01:11 remaining:              24%|##4       | 24/100 [elapsed: 01:21 remaining:              26%|##6       | 26/100 [elapsed: 01:43 remaining:              36%|###6      | 36/100 [elapsed: 01:54 remaining:              38%|###8      | 38/100 [elapsed: 02:17 remaining:              48%|####8     | 48/100 [elapsed: 02:29 remaining:              52%|#####2    | 52/100 [elapsed: 02:55 remaining:              60%|######    | 60/100 [elapsed: 03:08 remaining:              61%|######1   | 61/100 [elapsed: 03:21 remaining:              69%|######9   | 69/100 [elapsed: 03:34 remaining:              72%|#######2  | 72/100 [elapsed: 03:47 remaining:              75%|#######5  | 75/100 [elapsed: 04:01 remaining:              84%|########4 | 84/100 [elapsed: 04:14 remaining:              85%|########5 | 85/100 [elapsed: 04:27 remaining:              91%|#########1| 91/100 [elapsed: 04:40 remaining:              96%|#########6| 96/100 [elapsed: 04:54 remaining:              98%|#########8| 98/100 [elapsed: 05:07 remaining:             100%|##########| 100/100 [elapsed: 05:20 remaining:                                                                           100%|##########| 100/100 [elapsed: 05:20 remaining: 00:00]
[INFO   ]   calculating receptor contribution...
              0%|          | 0/100 [elapsed: 00:00 remaining: ?              8%|8         | 8/100 [elapsed: 00:03 remaining: 0              9%|9         | 9/100 [elapsed: 00:12 remaining: 0             11%|#1        | 11/100 [elapsed: 00:27 remaining:              17%|#7        | 17/100 [elapsed: 00:33 remaining:              20%|##        | 20/100 [elapsed: 00:39 remaining:              21%|##1       | 21/100 [elapsed: 00:45 remaining:              23%|##3       | 23/100 [elapsed: 00:59 remaining:              27%|##7       | 27/100 [elapsed: 01:06 remaining:              32%|###2      | 32/100 [elapsed: 01:13 remaining:              33%|###3      | 33/100 [elapsed: 01:20 remaining:              34%|###4      | 34/100 [elapsed: 01:27 remaining:              36%|###6      | 36/100 [elapsed: 01:34 remaining:              41%|####1     | 41/100 [elapsed: 01:42 remaining:              44%|####4     | 44/100 [elapsed: 01:49 remaining:              45%|####5     | 45/100 [elapsed: 01:56 remaining:              47%|####6     | 47/100 [elapsed: 02:03 remaining:              49%|####9     | 49/100 [elapsed: 02:10 remaining:              53%|#####3    | 53/100 [elapsed: 02:17 remaining:              57%|#####6    | 57/100 [elapsed: 02:25 remaining:              59%|#####8    | 59/100 [elapsed: 02:40 remaining:              64%|######4   | 64/100 [elapsed: 02:48 remaining:              68%|######8   | 68/100 [elapsed: 02:57 remaining:              69%|######9   | 69/100 [elapsed: 03:05 remaining:              71%|#######1  | 71/100 [elapsed: 03:13 remaining:              76%|#######6  | 76/100 [elapsed: 03:21 remaining:              77%|#######7  | 77/100 [elapsed: 03:29 remaining:              81%|########1 | 81/100 [elapsed: 03:38 remaining:              83%|########2 | 83/100 [elapsed: 03:46 remaining:              87%|########7 | 87/100 [elapsed: 03:54 remaining: 00:34]Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
             89%|########9 | 89/100 [elapsed: 04:02 remaining: 00:32]Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
             93%|#########3| 93/100 [elapsed: 04:11 remaining:              94%|#########3| 94/100 [elapsed: 04:19 remaining:              95%|#########5| 95/100 [elapsed: 04:27 remaining:              97%|#########7| 97/100 [elapsed: 04:45 remaining:              98%|#########8| 98/100 [elapsed: 04:54 remaining:              99%|#########9| 99/100 [elapsed: 05:03 remaining: 00:06]Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
            100%|##########| 100/100 [elapsed: 05:12 remaining:                                                                           100%|##########| 100/100 [elapsed: 05:12 remaining: 00:00]
[INFO   ]   calculating ligand contribution...
              0%|          | 0/100 [elapsed: 00:00 remaining: ?]Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL

gmx_MMPBSA.log

gmx_MMPBSA.log

Operating system

Zorin Os 17.1 (Ubuntu)

gmx_MMPBSA Version

gmx_MMPBSA v1.6.3

Python version

Python 3.10.14

Installation

conda AmberTools + pip

Valdes-Tresanco-MS commented 3 months ago

Looks like your ligand has a problem and Sander can't do the calculation. Please, check the mdout file.

JuuelHerza commented 3 months ago

Hello, thanks for the fast response, I see, so it is at the Sander level. Is there any command for gmx_MMPBSA to run only that step or to continue the gmx_MMPBSA calculation from that point?

I checked the file "_GMXMMPBSA_ligand_pb.mdout.0".

This is the last line after the "Atomic solvent accessible surface area":

PBSA BOMB: in MG: ncyc, itn, gid, norm, onorm 10 40 1 0.0253 0.0253

Valdes-Tresanco-MS commented 3 months ago

Sander fails here, so we must check what is going on. Please, send me the file to review them. Unfortunately, the calculation must be restarted from the beginning. Remember you can accelerate the calculation using MPI parallelization

JuuelHerza commented 3 months ago

Thanks, here's the file:

_GMXMMPBSA_ligand_pb.mdout.0.log

Valdes-Tresanco-MS commented 3 months ago

According to the gmx_MMPBSA.log file, you are performing a calculation for a protein-ligand system embedded in a membrane. Is that correct? Does it make sense to include the membrane in the calculation? In several cases, you can omit the membrane since it does not interact with the ligand. Could you compare all your systems with this one and try to identify the differences? Also, check the PBC

JuuelHerza commented 3 months ago

Hello, all the other complexes tested have the same protein (a porin) with a "pore blocker" and basically the same PBC treatment/MMPBSA input parameters have been used. I wanted to test later whether considering an implicit membrane would improve or decrease the correlation with experimental data. I'll try running a calculation now without the membrane.

First, I treat all the trajectories with this line:

echo -e "4\n0\n" | gmx trjconv -s rep${i}/step8_100ns.$ligand.rep${i}.tpr -f rep${i}/step8_100ns.$ligand.rep${i}.xtc -pbc mol -center -o rep${i}/step8_100ns.$ligand.rep${i}.PBCmol.xtc

Second, I remove the membrane since Parmed struggles with Lipid21 parameter conversion:

echo -e "$NoMEMB\n" | gmx trjconv -s rep${i}/step8_100ns.$ligand.rep${i}.tpr -f rep${i}/step8_100ns.$ligand.rep${i}.PBCmol.xtc -n index.ndx -o rep${i}/step8_100ns.$ligand.rep${i}.NoMEMB.xtc

Now that I notice, I don't know how the program knows the location of the membrane. Bu the orientation of the membrane is the same with all the other systems.

Valdes-Tresanco-MS commented 3 months ago

oh, I thought you had defined it. These parameters control the membrane center and thickness respectively mctrdz=94.1, and mthick=36.100.

JuuelHerza commented 3 months ago

Oh, yeah I did, a few months ago. I forgot I did some testing on a slightly different starting configuration. How embarrassing. I'm sorry, I think that's the cause, I just corrected the value mctrdz and it worked !

Thank you so much for your assistance.