PB calculations error

[Susana Cruz]

Hello, 

I've simulated a membrane protein complex system using gromacs, with inputs generated in CHARMM-GUI with charmm36m force field, and now I want to calcute the free energy of binding between the two proteins but I'm getting the following error:

"[INFO   ] 1000 frames were processed by cpptraj for use in calculation.
[INFO   ] Running calculations on normal system...
[INFO   ] Beginning PB calculations with /home/mbcstudent/.conda/envs/gmxMMPBSA/bin/sander
[INFO   ]   calculating complex contribution...
[ERROR  ] CalcError /home/mbcstudent/.conda/envs/gmxMMPBSA/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 "/home/mbcstudent/.conda/envs/gmxMMPBSA/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/home/mbcstudent/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 104, in gmxmmpbsa
    app.run_mmpbsa()
  File "/home/mbcstudent/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 203, in run_mmpbsa
    self.calc_list.run(rank, self.stdout)
  File "/home/mbcstudent/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/calculation.py", line 85, in run
    calc.run(rank, stdout=stdout, stderr=stderr)
  File "/home/mbcstudent/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/calculation.py", line 479, in run
    GMXMMPBSA_ERROR('%s failed with prmtop %s!\n\t' % (self.program, self.prmtop) +
  File "/home/mbcstudent/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/exceptions.py", line 169, in __init__
    raise exc(msg + '\nCheck the gmx_MMPBSA.log file to report the problem.')
CalcError: /home/mbcstudent/.conda/envs/gmxMMPBSA/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.
Error occurred on rank 8.
Exiting. All files have been retained."

I removed pbc and fitted my trajectory, so I don't know what the problem is.

This is the We Transfer link with my input files: https://we.tl/t-tFWudM9Iiu

The comand line used was: mpirun -np 20 gmx_MMPBSA MPI -O -i mmpbsa.in -cs prod.tpr -ci prod2.ndx -cg 6 5 -ct fskip10_nosolv.xtc -cp topol.top -o prod.dat -do prod_decomp.dat -eo prod_energy.csv -deo prod_energy_decomp.csv

Thank you in advance.

[Mario Ernesto]

could you send the files as a google drive attachment? I am pretty sure the biggest file is the trajectory but, we only need a few frames of the trajectory 

--
You received this message because you are subscribed to the Google Groups "gmx_MMPBSA" group.
To unsubscribe from this group and stop receiving emails from it, send an email to gmx_mmpbsa+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/gmx_mmpbsa/69a5798b-f64a-494f-a721-b19afeaaa3f3n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Susana Cruz]

Here's the link to my google drive folder with the files:

https://drive.google.com/drive/folders/1w0aQWZ86vFg8_XZSgiun_Jat0DwziJeB?usp=sharing

[Susana Cruz]

I have since changed some of my files and parameters in the mmpbsa.in file, as well as change command line to only do pb calculations (and not energy decomposition) and now I'm getting this error in my mdout.file :

 "error at /home/conda/feedstock_root/build_artifacts/ambertools_1645002103866/work/ambertools/src/pbsa/pb_bldsys.f90 , line 1370"

I can't find anything about this error, so I don't know what I'm doing wrong as I'm pretty new to this.

Thank you for helping me, I'll await your response.

On Sunday, July 10, 2022 at 8:01:20 PM UTC+1 marioe911116 wrote:

[Mario Sergio Valdes]

Hi Susana, can you attach in Google drive the topology file along with the folder of the itp and prm files (toppar or charmm36*)?

[Susana Cruz]

Yes, I've attached them now.

[Mario Sergio Valdes]

I have reproduced the error. Sander crashes because the configuration you are using consumes a significant amount of RAM. In my case, it crashes after consuming more than 18GB in a single thread, so using 20 threads should consume more than 400GB. However, the good news is that this configuration is not explicitly necessary. The parameter that generates this excessive consumption is maxarcdot. In the example we have in the documentation it appears as 15000, which is certainly a bit large. In your case, 5000 works correctly (for membrane,  in other cases you can use the default value 1500) and decreases considerably the RAM required for the calculation. Additionally, the fillratio value generates a considerable RAM consumption so that you can decrease it for your system as well (1.25 is commonly used for membrane, but you can probe how much RAM consumes).
Notes:

1. In this case, it is possible to dispense with the membrane calculation, since it is apparently distant (>15 A) from the interaction site. To do this, you must simply generate a group in the index that excludes the region across the membrane, say up to the residue SER241. This would allow you to use the most common parameters to estimate the interaction energy of the complex without sacrificing accuracy due to the high consumption of computational resources used to compute the membrane.
2. To compute the membrane you are using some incorrect parameters. mctrdz determines where the membrane will be located in the z-coordinate, so you should set it correctly in the center of the membrane of your system, according to my estimation at approximately 102.190. mthick on the other hand determines the thickness, so it should also be set correctly, according to my estimation it would be 47.0 approximately. 
3. Note that computing the membrane is extremely expensive, so the configuration must be modified and several parameters used for the same purpose are already modified by default (e.g. mprob = 2.7 and prbrad = 1.4). If you decide to use it, you can start from the attached configuration (mmpbsa_mem.in), otherwise, you can start from the default configuration of the file that is generated with gmx_MMPBSA --create_input pb
   
4. Keep in mind that PB is computationally expensive, so I suggest that you estimate how much it consumes doing the calculation with only one thread and for 2 or 3 frames. With the attached configuration, you can consume a maximum of 10GB per thread approximately. This means that you must take into account how much RAM you have before sending the number of CPUs. Let's say, if the node you are assigned is 128 GB, only, you will be able to send a job of maximum 11 (or 12?) CPUs.
5. Changing the fillratio variable from 1.25 to 2 generates an additional 5.5GB consumption.
   

HTH!

Mario S.

[Susana Cruz]

Hi again, 

I took your advice and created index groups for the protein and complex without the membrane part. I tried to change the trajetory in order to only include the new complex group but it caused an error, so I used my original trajectory and it seemed to be working fine. But while it hasn't crashed, there hasn't been any progress in the last 5 hours. Is this normal? I changed almost all parameters in the mmpbsa.in file to the default values, using the same command line and trajectory file, so how long should I expect it to take? 

Thank you for your help

[Mario Sergio Valdes]

It is difficult for me to know how long it will take. However, you can open the *.mdout.0 file and count the number of frames that have been processed (minimizing coord set #) and estimate the time as follow:

                            # processed frames

Time 1 frame = --------------------------------   

                               calculation time

                                                                    # Total frames

Estimated time = Time 1 frame * ---------------------------------

                                                                       # CPUs

Example: 

Time 1 frame = 10 min

# CPUs = 20

# Total frames = 1000

Estimated time = 10 min * 1000/20 = 500 min = 8,333333333 hours

Please note that this estimate is for the complex only. The total time would be this time x2 approximately. In this example, the total time would be approximately 16 or 17 hours.

We are planning to implement a system for tracking and estimating computation time, but it could still take some time.

Mario S.

[Mario Sergio Valdes]

Hi Susan. 

Recently, we have implemented a progress bar that shows the elapsed time, the number of frames processed, and the estimated time remaining. To use it, simply update from Github (python -m pip install git+https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA.git -U). It is in the beta phase but seems completely stable (at least in most of the calculations we have tested). Let me know if you find it useful

Mario S.