search

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
226 stars 65 forks source link

[Bug-gmx_MMPBSA]: function `get_corrstd` fails in rare case #401

Closed jsjyhzy closed 1 year ago

jsjyhzy commented 1 year ago

Bug summary

In a very very rare case, the square root seems fail for the negative number, ends up with the ValueError: math domain error. This particular formula looks like the absolute value of the deviation between val1 and val2, I think it would be more robust to use abs(val1 - val2) instead of this, and using abs(..) will be little bit faster

>>> timeit.timeit('a=random(); b=random(); abs(a-b)', setup='from random import random; from math import sqrt')
0.05487921601161361
>>> timeit.timeit('a=random(); b=random(); sqrt(a**2 + b**2 -2*a*b)', setup='from random import random; from math import sqrt')
0.11629197793081403

Terminal output

[INFO   ]   calculating complex contribution...

              0%|          | 0/2 [elapsed: 00:00 remaining: ?]
            100%|##########| 2/2 [elapsed: 00:10 remaining: 00:00]

            100%|##########| 2/2 [elapsed: 00:10 remaining: 00:00]
[INFO   ]   calculating receptor contribution...

              0%|          | 0/2 [elapsed: 00:00 remaining: ?]
            100%|##########| 2/2 [elapsed: 00:06 remaining: 00:00]

            100%|##########| 2/2 [elapsed: 00:06 remaining: 00:00]
[INFO   ]   calculating ligand contribution...

              0%|          | 0/2 [elapsed: 00:00 remaining: ?]
            100%|##########| 2/2 [elapsed: 00:01 remaining: 00:00]

            100%|##########| 2/2 [elapsed: 00:01 remaining: 00:00]
[INFO   ] Parsing results to output files...

  File "/public/home/neotrident/apps/mmpbsa_1.6.1/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 109, in gmxmmpbsa
    app.parse_output_files()
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 1226, in parse_output_files
    self.calc_types.normal[key]['delta'] = BindingStatistics(self.calc_types.normal[key]['complex'],
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/amber_outputs.py", line 977, in __init__
    self._delta()
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/amber_outputs.py", line 1004, in _delta
    self[key] = temp.corr_sub(self.lig[key])
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/utils.py", line 114, in corr_sub
    comp_std = get_corrstd(self_std, other_std)
  File "/public/home/neotrident/apps/mmpbsa_1.6.1/lib/python3.9/site-packages/GMXMMPBSA/utils.py", line 158, in get_corrstd
    return sqrt(val1 ** 2 + val2 ** 2 - 2 * val1 * val2)
ValueError: math domain error
Error occurred on rank 0.
Exiting. All files have been retained.

gmx_MMPBSA.log

[INFO   ] Starting gmx_MMPBSA v1.6.1
[DEBUG  ] WDIR          : /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] AMBERHOME     : /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] PYTHON EXE    : /public/home/neotrident/apps/mmpbsa_1.6.1/bin/python
[DEBUG  ] PYTHON VERSION: 3.9.15 | packaged by conda-forge | (main, Nov 22 2022, 08:45:29) [GCC 10.4.0]
[DEBUG  ] MPI           : /public/home/neotrident/apps/mmpbsa_1.6.1/bin/mpirun
[DEBUG  ] ParmEd        : 3.4.3+11.g41cc9ab
[DEBUG  ] OS PLATFORM   : Linux-3.10.0-957.el7.x86_64-x86_64-with-glibc2.17
[DEBUG  ] OS SYSTEM     : Linux
[DEBUG  ] OS VERSION    : #1 SMP Fri Oct 25 22:34:08 UTC 2019
[DEBUG  ] OS PROCESSOR  : x86_64

[INFO   ] Command-line
  mpirun -np 2 gmx_MMPBSA MPI -O -i mmpbsa.in -ci index.ndx -cg 0 1 -ct charmm-gui_input_bio.xtc -cs charmm-gui_input_bio.tpr -cp charmm-gui_input_bio.top -cr reference.pdb -nogui -eo mmpbsa_ana.csv -deo mmpbsa_dec.csv

[DEBUG  ] |Input file:
[DEBUG  ] |--------------------------------------------------------------
[DEBUG  ] |&general
[DEBUG  ] | startframe=0
[DEBUG  ] | endframe=1000000000
[DEBUG  ] | interval=10
[DEBUG  ] | temperature=298.15
[DEBUG  ] | interaction_entropy=0
[DEBUG  ] | ie_segment=25
[DEBUG  ] | c2_entropy=0
[DEBUG  ] | qh_entropy=0
[DEBUG  ] | exp_ki=0
[DEBUG  ] |
[DEBUG  ] |/
[DEBUG  ] |
[DEBUG  ] |&gb
[DEBUG  ] | igb=5
[DEBUG  ] | saltcon=0.15
[DEBUG  ] | surfoff=0
[DEBUG  ] | surften=0.0072
[DEBUG  ] | intdiel=1
[DEBUG  ] | extdiel=78.5
[DEBUG  ] |
[DEBUG  ] |/
[DEBUG  ] |--------------------------------------------------------------
[DEBUG  ] 
[INFO   ] Checking mmpbsa.in input file...
[WARNING] &decomp namelist has not been defined in the input file. Ignoring '-deo' flag... 
[WARNING] The startframe variable must be >= 1. Changing startframe from 0 to 1
[INFO   ] Checking mmpbsa.in input file...Done.

[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /public/home/neotrident/apps/mmpbsa_1.6.1/bin/cpptraj
[INFO   ] tleap found! Using /public/home/neotrident/apps/mmpbsa_1.6.1/bin/tleap
[INFO   ] parmchk2 found! Using /public/home/neotrident/apps/mmpbsa_1.6.1/bin/parmchk2
[INFO   ] sander found! Using /public/home/neotrident/apps/mmpbsa_1.6.1/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[INFO   ] Checking external programs...Done.

[INFO   ] Building AMBER topologies from GROMACS files...
[INFO   ] Get PDB files from GROMACS structures files...
[INFO   ] Making gmx_MMPBSA index for complex...
[DEBUG  ] Running command: echo -e "name 0 GMXMMPBSA_REC\n name 1 GMXMMPBSA_LIG\n  0 | 1\n q\n" | /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx make_ndx -n index.ndx -o _GMXMMPBSA_COM_index.ndx -f charmm-gui_input_bio.tpr
[DEBUG  ]                :-) GROMACS - gmx make_ndx, 2022.4-conda_forge (-:
[DEBUG  ] 
[DEBUG  ] Executable:   /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[DEBUG  ] Data prefix:  /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] Working dir:  /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] Command line:
[DEBUG  ]   gmx make_ndx -n index.ndx -o _GMXMMPBSA_COM_index.ndx -f charmm-gui_input_bio.tpr
[DEBUG  ] 
[DEBUG  ] 
[DEBUG  ] Reading structure file
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] 
[DEBUG  ] GROMACS reminds you: "Give a Man a Fish" (Arrested Development)
[DEBUG  ] 
[DEBUG  ] Going to read 1 old index file(s)
[DEBUG  ] 
[DEBUG  ]   0 new_group_0         :  6098 atoms
[DEBUG  ]   1 new_group_1         :   881 atoms
[DEBUG  ] 
[DEBUG  ]  nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
[DEBUG  ]  'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
[DEBUG  ]  't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
[DEBUG  ]  'r': residue              'res' nr         'chain' char
[DEBUG  ]  "name": group             'case': case sensitive           'q': save and quit
[DEBUG  ]  'ri': residue index
[DEBUG  ] 
[DEBUG  ] > 
[DEBUG  ] 
[DEBUG  ] > 
[DEBUG  ] 
[DEBUG  ] > 
[DEBUG  ] Copied index group 0 'GMXMMPBSA_REC'
[DEBUG  ] Copied index group 1 'GMXMMPBSA_LIG'
[DEBUG  ] Merged two groups with OR: 6098 881 -> 6979
[DEBUG  ] 
[DEBUG  ]   2 GMXMMPBSA_REC_GMXMMPBSA_LIG:  6979 atoms
[DEBUG  ] 
[DEBUG  ] > 
[INFO   ] Normal Complex: Saving group new_group_0_new_group_1 (0_1) in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb
[DEBUG  ] Running command: echo -e "GMXMMPBSA_REC_GMXMMPBSA_LIG"| /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_COM.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ]                :-) GROMACS - gmx trjconv, 2022.4-conda_forge (-:
[DEBUG  ] 
[DEBUG  ] Executable:   /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[DEBUG  ] Data prefix:  /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] Working dir:  /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] Command line:
[DEBUG  ]   gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_COM.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ] 
[DEBUG  ] Will write pdb: Protein data bank file
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Group     0 (  GMXMMPBSA_REC) has  6098 elements
[DEBUG  ] Group     1 (  GMXMMPBSA_LIG) has   881 elements
[DEBUG  ] Group     2 (GMXMMPBSA_REC_GMXMMPBSA_LIG) has  6979 elements
[DEBUG  ] Select a group: 
Reading frame       0 time    0.000   
[DEBUG  ] Precision of charmm-gui_input_bio.xtc is 0.001 (nm)
[DEBUG  ] 
Reading frame       1 time    1.000   
[DEBUG  ] Dumping frame at t= 0 ps
[DEBUG  ] Last written: frame      0 time    0.000
[DEBUG  ] 
[DEBUG  ] 
[DEBUG  ] GROMACS reminds you: "Give a Man a Fish" (Arrested Development)
[DEBUG  ] 
[DEBUG  ] Note that major changes are planned in future for trjconv, to improve usability and utility.
[DEBUG  ] Select group for output
[DEBUG  ] Selected 2: 'GMXMMPBSA_REC_GMXMMPBSA_LIG'
[INFO   ] No receptor structure file was defined. Using ST approach...
[INFO   ] Using receptor structure from complex to generate AMBER topology
[INFO   ] Normal Receptor: Saving group new_group_0 (0) in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_REC.pdb
[DEBUG  ] Running command: echo -e "0"| /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_REC.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ]                :-) GROMACS - gmx trjconv, 2022.4-conda_forge (-:
[DEBUG  ] 
[DEBUG  ] Executable:   /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[DEBUG  ] Data prefix:  /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] Working dir:  /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] Command line:
[DEBUG  ]   gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_REC.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ] 
[DEBUG  ] Will write pdb: Protein data bank file
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Group     0 (  GMXMMPBSA_REC) has  6098 elements
[DEBUG  ] Group     1 (  GMXMMPBSA_LIG) has   881 elements
[DEBUG  ] Group     2 (GMXMMPBSA_REC_GMXMMPBSA_LIG) has  6979 elements
[DEBUG  ] Select a group: 
Reading frame       0 time    0.000   
[DEBUG  ] Precision of charmm-gui_input_bio.xtc is 0.001 (nm)
[DEBUG  ] 
Reading frame       1 time    1.000   
[DEBUG  ] Dumping frame at t= 0 ps
[DEBUG  ] Last written: frame      0 time    0.000
[DEBUG  ] 
[DEBUG  ] 
[DEBUG  ] GROMACS reminds you: "Give a Man a Fish" (Arrested Development)
[DEBUG  ] 
[DEBUG  ] Note that major changes are planned in future for trjconv, to improve usability and utility.
[DEBUG  ] Select group for output
[DEBUG  ] Selected 0: 'GMXMMPBSA_REC'
[INFO   ] No ligand structure file was defined. Using ST approach...
[INFO   ] Using ligand structure from complex to generate AMBER topology
[INFO   ] Normal Ligand: Saving group new_group_1 (1) in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[DEBUG  ] Running command: echo -e "1"| /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_LIG.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ]                :-) GROMACS - gmx trjconv, 2022.4-conda_forge (-:
[DEBUG  ] 
[DEBUG  ] Executable:   /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[DEBUG  ] Data prefix:  /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] Working dir:  /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] Command line:
[DEBUG  ]   gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o _GMXMMPBSA_LIG.pdb -n _GMXMMPBSA_COM_index.ndx -dump 0
[DEBUG  ] 
[DEBUG  ] Will write pdb: Protein data bank file
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Group     0 (  GMXMMPBSA_REC) has  6098 elements
[DEBUG  ] Group     1 (  GMXMMPBSA_LIG) has   881 elements
[DEBUG  ] Group     2 (GMXMMPBSA_REC_GMXMMPBSA_LIG) has  6979 elements
[DEBUG  ] Select a group: 
Reading frame       0 time    0.000   
[DEBUG  ] Precision of charmm-gui_input_bio.xtc is 0.001 (nm)
[DEBUG  ] 
Reading frame       1 time    1.000   
[DEBUG  ] Dumping frame at t= 0 ps
[DEBUG  ] Last written: frame      0 time    0.000
[DEBUG  ] 
[DEBUG  ] 
[DEBUG  ] GROMACS reminds you: "Give a Man a Fish" (Arrested Development)
[DEBUG  ] 
[DEBUG  ] Note that major changes are planned in future for trjconv, to improve usability and utility.
[DEBUG  ] Select group for output
[DEBUG  ] Selected 1: 'GMXMMPBSA_LIG'
[INFO   ] Checking the structures consistency...
[INFO   ] Assigning chain ID to structures files according to the reference structure...
[INFO   ] 
[INFO   ] Using topology conversion. Setting radiopt = 0...
[INFO   ] Building Normal Complex Amber topology...
[INFO   ] Detected Amber/OPLS force field topology format...
[INFO   ] Assigning PBRadii mbondi2 to Complex...
[INFO   ] Writing Normal Complex AMBER topology...
[INFO   ] No Receptor topology file was defined. Using ST approach...
[INFO   ] Building AMBER Receptor topology from Complex...
[INFO   ] Assigning PBRadii mbondi2 to Receptor...
[INFO   ] Writing Normal Receptor AMBER topology...
[INFO   ] No Ligand topology file was defined. Using ST approach...
[INFO   ] Building AMBER Ligand topology from Complex...
[INFO   ] Assigning PBRadii mbondi2 to Ligand...
[INFO   ] Writing Normal Ligand AMBER topology...
[INFO   ] Cleaning normal complex trajectories...
[DEBUG  ] Running command: echo -e "GMXMMPBSA_REC_GMXMMPBSA_LIG"| /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o COM_traj_0.xtc -n _GMXMMPBSA_COM_index.ndx
[DEBUG  ]                :-) GROMACS - gmx trjconv, 2022.4-conda_forge (-:
[DEBUG  ] 
[DEBUG  ] Executable:   /public/home/neotrident/apps/mmpbsa_1.6.1/bin.AVX2_256/gmx
[DEBUG  ] Data prefix:  /public/home/neotrident/apps/mmpbsa_1.6.1
[DEBUG  ] Working dir:  /public/home/neotrident/job/271/9c0703ad-cc24-4d0c-b205-b253e2442c2b/1691156885252/16914212118381350131112400
[DEBUG  ] Command line:
[DEBUG  ]   gmx trjconv -f charmm-gui_input_bio.xtc -s charmm-gui_input_bio.tpr -o COM_traj_0.xtc -n _GMXMMPBSA_COM_index.ndx
[DEBUG  ] 
[DEBUG  ] Will write xtc: Compressed trajectory (portable xdr format): xtc
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Reading file charmm-gui_input_bio.tpr, VERSION 2021.5 (single precision)
[DEBUG  ] Group     0 (  GMXMMPBSA_REC) has  6098 elements
[DEBUG  ] Group     1 (  GMXMMPBSA_LIG) has   881 elements
[DEBUG  ] Group     2 (GMXMMPBSA_REC_GMXMMPBSA_LIG) has  6979 elements
[DEBUG  ] Select a group: 
Reading frame       0 time    0.000   
[DEBUG  ] Precision of charmm-gui_input_bio.xtc is 0.001 (nm)
[DEBUG  ] Using output precision of 0.001 (nm)
[DEBUG  ] 
Reading frame       1 time    1.000    ->  frame      0 time    0.000      

Reading frame       2 time    2.000    ->  frame      1 time    1.000      

Reading frame       3 time    3.000    ->  frame      2 time    2.000      

Reading frame       4 time    4.000    ->  frame      3 time    3.000      

Reading frame       5 time    5.000    ->  frame      4 time    4.000      

Reading frame       6 time    6.000    ->  frame      5 time    5.000      

Reading frame       7 time    7.000    ->  frame      6 time    6.000      

Reading frame       8 time    8.000    ->  frame      7 time    7.000      

Reading frame       9 time    9.000    ->  frame      8 time    8.000      

Reading frame      10 time   10.000    ->  frame      9 time    9.000      

Last frame         10 time   10.000   
[DEBUG  ] Last written: frame     10 time   10.000
[DEBUG  ] 
[DEBUG  ] 
[DEBUG  ] GROMACS reminds you: "The scientist is not a person who gives the right answers, he's one who asks the right questions." (Claude Levi-Strauss)
[DEBUG  ] 
[DEBUG  ] Note that major changes are planned in future for trjconv, to improve usability and utility.
[DEBUG  ] Select group for output
[DEBUG  ] Selected 2: 'GMXMMPBSA_REC_GMXMMPBSA_LIG'
[INFO   ] Building AMBER topologies from GROMACS files... Done.

[INFO   ] Loading and checking parameter files for compatibility...
[INFO   ] Preparing trajectories for simulation...

[INFO   ] 2 frames were processed by cpptraj for use in calculation.
[INFO   ] Starting calculations in 2 CPUs...
[INFO   ] Running calculations on normal system...
[INFO   ] Beginning GB calculations with /public/home/neotrident/apps/mmpbsa_1.6.1/bin/sander
[INFO   ]   calculating complex contribution...
[INFO   ]   calculating receptor contribution...
[INFO   ]   calculating ligand contribution...
[INFO   ] Parsing results to output files...

Operating system

No response

gmx_MMPBSA Version

No response

Python version

No response

Installation

None

Valdes-Tresanco-MS commented 1 year ago

2 frames are too few. Please define more (5-10 is a good number).