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
211 stars 64 forks source link

[Bug-gmx_MMPBSA]: pb_fdfrc(): Atom out of focusing box #476

Closed lachrymator closed 5 months ago

lachrymator commented 5 months ago

Bug summary

Conducting MMPBSA calculations fails on COM.prmtop for some reason. MMGBSA proceeds to completion without issue. Looking in the _GMX_complex_pb.out it seems there's something off with the focusing box?

Terminal output

[INFO   ] 30 frames were processed by cpptraj for use in calculation.
[INFO   ] Starting calculations in 1 CPUs...

[INFO   ] Running calculations on normal system...
[INFO   ] Beginning GB calculations with /shared/amber22/bin/sander
[INFO   ]   calculating complex contribution...
            100%|####################################################################################################################################################| 30/30 [elapsed: 01:16 remaining: 00:00]
[INFO   ]   calculating receptor contribution...
            100%|####################################################################################################################################################| 30/30 [elapsed: 01:05 remaining: 00:00]
[INFO   ]   calculating ligand contribution...
            100%|####################################################################################################################################################| 30/30 [elapsed: 00:15 remaining: 00:00]
[INFO   ] Beginning PB calculations with /shared/amber22/bin/sander
[INFO   ]   calculating complex contribution...
              0%|                                                                                                                                                         | 0/30 [elapsed: 00:00 remaining: ?][ERROR  ] CalcError

/shared/amber22/bin/sander failed with prmtop COM.prmtop!

If you are using sander and PB calculation, check the *.mdout files to get the sander error

Check the gmx_MMPBSA.log file to report the problem.
  File "/shared/miniconda3/envs/mmpbsa/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/shared/miniconda3/envs/mmpbsa/lib/python3.10/site-packages/GMXMMPBSA/app.py", line 101, in gmxmmpbsa
    app.run_mmpbsa()
  File "/shared/miniconda3/envs/mmpbsa/lib/python3.10/site-packages/GMXMMPBSA/main.py", line 205, in run_mmpbsa
    self.calc_list.run(rank, self.stdout)
  File "/shared/miniconda3/envs/mmpbsa/lib/python3.10/site-packages/GMXMMPBSA/calculation.py", line 142, in run
    calc.run(rank, stdout=stdout, stderr=stderr)
  File "/shared/miniconda3/envs/mmpbsa/lib/python3.10/site-packages/GMXMMPBSA/calculation.py", line 625, in run
    GMXMMPBSA_ERROR('%s failed with prmtop %s!\n\t' % (self.program, self.prmtop) +
  File "/shared/miniconda3/envs/mmpbsa/lib/python3.10/site-packages/GMXMMPBSA/exceptions.py", line 171, in __init__
    raise exc('\n\n' + msg + '\n\nCheck the gmx_MMPBSA.log file to report the problem.')
CalcError:

/shared/amber22/bin/sander failed with prmtop COM.prmtop!

If you are using sander and PB calculation, check the *.mdout files to get the sander error

Check the gmx_MMPBSA.log file to report the problem.
Exiting. All files have been retained.

gmx_MMPBSA.log

[DEBUG  ] Using output precision of 0.001 (nm)
Last frame       3000 time 3000.000    ->  frame   2999 time 2999.000
[DEBUG  ] Last written: frame   3000 time 3000.000
[DEBUG  ]
[DEBUG  ]
[DEBUG  ] GROMACS reminds you: "In physics, you don't have to go around making trouble for yourself. Nature does it for you." (Frank Wilczek)
[DEBUG  ]
[DEBUG  ] Note that major changes are planned in future for trjconv, to improve usability and utility.
[DEBUG  ] Select group for output
[DEBUG  ] Selected 21: '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   ] 30 frames were processed by cpptraj for use in calculation.
[INFO   ] Starting calculations in 1 CPUs...
[INFO   ] Running calculations on normal system...
[INFO   ] Beginning GB calculations with /shared/amber22/bin/sander
[INFO   ]   calculating complex contribution...
[INFO   ]   calculating receptor contribution...
[INFO   ]   calculating ligand contribution...
[INFO   ] Beginning PB calculations with /shared/amber22/bin/sander
[INFO   ]   calculating complex contribution...
[ERROR  ] CalcError

/shared/amber22/bin/sander failed with prmtop COM.prmtop!

If you are using sander and PB calculation, check the *.mdout files to get the sander error

_GMX_complex_pb.mdout.0

          -------------------------------------------------------
          Amber 22 SANDER                              2022
          -------------------------------------------------------

| Run on 03/25/2024 at 01:19:56

|   Executable path: /shared/amber22/bin/sander
| Working directory: /shared/DEMO
|          Hostname: Unknown
  [-O]verwriting output

File Assignments:
|  MDIN: _GMXMMPBSA_pb.mdin
| MDOUT: _GMXMMPBSA_complex_pb.mdout.0
|INPCRD: _GMXMMPBSA_dummycomplex.inpcrd
|  PARM: COM.prmtop
|RESTRT: _GMXMMPBSA_restrt.0
|  REFC: refc
| MDVEL: mdvel
| MDFRC: mdfrc
|  MDEN: mden
| MDCRD: mdcrd
|MDINFO: mdinfo
|  MTMD: mtmd
|INPDIP: inpdip
|RSTDIP: rstdip
|INPTRA: _GMXMMPBSA_complex.mdcrd.0

 Here is the input file:

File generated by gmx_MMPBSA
&cntrl
 ntb=0, nsnb=99999, cut=999.0, imin=5,
 ipb=2, inp=1, dec_verbose=0,
/
&pb
 istrng=150.0, radiopt=0, maxitn=1000,
 fillratio=4.0, epsmem=4.0,
/

--------------------------------------------------------------------------------
   1.  RESOURCE   USE:
--------------------------------------------------------------------------------

| Flags:
 PBSA Warning in pb_read(): sprob=0.557 is optimized for inp=2 and  should not be used with inp=1. It has been reset to 1.4.
 PBSA Warning in pb_read(): cavity_surften=0.0378 is optimized for inp=2  and should not be used with inp=1. It has been reset to 0.005000.
 PBSA Warning in pb_read(): cavity_offset=-0.5692 is optimized for inp=2  and should not be used with inp=1. It has been reset to 0.000.
| New format PARM file being parsed.
| Version =    1.000 Date = 03/25/24 Time = 01:13:40
 NATOM  =    4406 NTYPES =    4406 NBONH =    2181 MBONA  =    2279
 NTHETH =    4961 MTHETA =    3076 NPHIH =   10970 MPHIA  =   10467
 NHPARM =       0 NPARM  =       0 NNB   =   24313 NRES   =     269
 NBONA  =    2279 NTHETA =    3076 NPHIA =   10467 NUMBND =      54
 NUMANG =      51 NPTRA  =     349 NATYP =       1 NPHB   =       0
 IFBOX  =       0 NMXRS  =      55 IFCAP =       0 NEXTRA =       0
 NCOPY  =       0

 Implicit solvent radii are H(N)-modified Bondi radii (mbondi2)

|     Memory Use     Allocated
|     Real            20494474
|     Hollerith          13489
|     Integer         19687698
|     Max Pairs              1
|     nblistReal             0
|     nblist Int             0
|       Total           237070 kbytes

| Note: 1-4 EEL scale factors are being read from the topology file.

| Note: 1-4 VDW scale factors are being read from the topology file.
| Duplicated    0 dihedrals
| Duplicated    0 dihedrals

--------------------------------------------------------------------------------
   2.  CONTROL  DATA  FOR  THE  RUN
--------------------------------------------------------------------------------

General flags:
     imin    =       5, nmropt  =       0

Nature and format of input:
     ntx     =       1, irest   =       0, ntrx    =       1

Nature and format of output:
     ntxo    =       2, ntpr    =      50, ntrx    =       1, ntwr    =       1
     iwrap   =       0, ntwx    =       0, ntwv    =       0, ntwe    =       0
     ioutfm  =       1, ntwprt  =       0, idecomp =       0, rbornstat=      0

Potential function:
     ntf     =       1, ntb     =       0, igb     =      10, nsnb    =   99999
     ipol    =       0, gbsa    =       0, iesp    =       0
     dielc   =   1.00000, cut     = 999.00000, intdiel =   1.00000

Frozen or restrained atoms:
     ibelly  =       0, ntr     =       0

Energy minimization:
     maxcyc  =       1, ncyc    =      10, ntmin   =       1
     dx0     =   0.01000, drms    =   0.00010
|  INFO: Old style inpcrd file read

--------------------------------------------------------------------------------
   3.  ATOMIC COORDINATES AND VELOCITIES
--------------------------------------------------------------------------------

Cpptraj Generated Restart
 begin time read from input coords =     4.000 ps

 Number of triangulated 3-point waters found:        0

--------------------------------------------------------------------------------
   4.  RESULTS
--------------------------------------------------------------------------------

 POST-PROCESSING OF TRAJECTORY ENERGIES
Cpptraj Generated trajectory
minimizing coord set #     1
     34.6300     50.2900     17.7800
pb_fdfrc(): Atom out of focusing box    11    15     7

Operating system

Ubuntu 22

gmx_MMPBSA Version

gmx_MMPBSA v1.6.3 based on MMPBSA version 16.0 and AmberTools 20

Python version

Python 3.10.13

Installation

AmberTools compilation + conda

marioernestovaldes commented 5 months ago

Could you please send the files you are using?

lachrymator commented 5 months ago

https://drive.google.com/file/d/1QznSylN24jKFXuhE6PzNInBpSMd4hpFr/view?usp=sharing

The command used to run is gmx_MMPBSA -O -i mmpbsa.in -cs md_0_1.pdb -ct md_0_1.xtc -ci index.ndx -cg 1 19 -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv -nogui

marioernestovaldes commented 5 months ago

this is extremely odd, I haven't found any solution to the issue so far, but I'll keep working on it to see if I can find something...

lachrymator commented 5 months ago

Much obliged! Will definitely sponsor!

Valdes-Tresanco-MS commented 5 months ago

What program do you use to generate the topology? Looks like the topology file was generated considering each atom as an independent atom. It means that every carbon (same for O, N, etc.) atom has a set of redundant parameters. Maybe, this is not related, but this consumes a lot of resources and generates a heavy topology file for gromacs and amber. image

I got the same error for the receptor and the ligand. I don't know why this happened, but I would like to perform the calculation using a more standard topology file. I feel that this is the problem. Maybe you can parameterize your system using CHARMM-GUI or Antechamber+tLeap and test if this solves the problem

lachrymator commented 5 months ago

Thank you for looking at that!

I used OpenFF to generate the combined topology in a manner similar to: https://docs.openforcefield.org/en/latest/examples/openforcefield/openff-toolkit/toolkit_showcase/toolkit_showcase.html

I don't have access to CHARMM FF or other commercial FFs

Valdes-Tresanco-MS commented 5 months ago

CHARMM-GUI is free for academic user. AmberTools is free for academic and commercial use. I will take a look to OpenFF because we don't support it out of the box

lachrymator commented 5 months ago

Allright, fair enough. Thank you very much.

lachrymator commented 5 months ago

Ok I can confirm using Antechamber / t-Leap makes everything work fine. Also, it appears to be much faster.

Valdes-Tresanco-MS commented 5 months ago

Reviewing OpenFF I see that it is not optimized yet. It generates a lot of redundant information. BTW, thank you for your support. You are the first and it represents a lot to us.