QM- calculation error

[maral afshinpour]

Hello, I am new to this program and when I run gmx_MMPBSA for the sample file provided on your website, the calculation was completed successfully without any error. 

but when I am running for my trajectory file, it got the error you can find in error.txt as attached.

[INFO   ] Starting
[WARNING] protein_forcefield and ligand_forcefield variables are deprecate since version 1.4.1 and will be remove in the next version. Please, use forcefield instead.
[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/cpptraj
[INFO   ] tleap found! Using /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/tleap
[INFO   ] parmchk2 found! Using /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/parmchk2
[INFO   ] sander found! Using /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /gpfs/shared/apps_local/gromacs/2019.6-nompi/bin/gmx
[INFO   ] Checking external programs...Done.

[INFO   ] Building AMBER Topologies from GROMACS files...
[INFO   ] Checking if supported force fields exists in Amber data...
[INFO   ] Get PDB files from GROMACS structures files...
[INFO   ] Making gmx_MMPBSA index for complex...
[INFO   ] Normal Complex: Saving group 10_11 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb
[INFO   ] No receptor structure file was defined. Using ST approach...
[INFO   ] Using receptor structure from complex to generate AMBER topology
[INFO   ] Normal Complex: Saving group 10 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_REC.pdb
[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 11 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[INFO   ] Building Normal Complex Amber Topology...
[INFO   ] Writing Normal Complex Amber Topology...
[INFO   ] No Receptor topology files was defined. Using ST approach...
[INFO   ] Building AMBER Receptor Topology from Complex...
[INFO   ] No Ligand Topology files was defined. Using ST approach...
[INFO   ] Building AMBER Ligand Topology from Complex...
[INFO   ] Cleaning normal complex trajectories...
[INFO   ] Building AMBER Topologies from GROMACS files...Done.

[INFO   ] Loading and checking parameter files for compatibility...

CHAMBER prmtops found. Forcing use of sander
Preparing trajectories for simulation...
11 frames were processed by cpptraj for use in calculation.

Running calculations on normal system...

Beginning GB calculations with /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/sander
  calculating complex contribution...
  calculating receptor contribution...
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/app.py", line 100, in gmxmmpbsa
    app.run_mmpbsa()
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/main.py", line 212, in run_mmpbsa
    self.calc_list.run(rank, self.stdout)
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/calculation.py", line 84, in run
    calc.run(rank, stdout=stdout, stderr=stderr)
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/calculation.py", line 158, in run
    raise CalcError('%s failed with prmtop %s!' % (self.program,
CalcError: /gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/bin/sander failed with prmtop REC.prmtop!

You can also find all the required files attached. and the commands I used is in run.slurm.  I want to calculate the binding affinity and qm between two chains of protein that I specified residues for qm in mmpbsa.in file. the version of gromacs 2019 and gmx_MMPBSA is v.1.4.1. and I used charmm36-jul2021 forcefield for md simulation.

[maral afshinpour]

I attached all files here.

[Mario Sergio Valdes]

Please attach the trajectory (10 frames), the charmm parameters folder (charmm36-jul2021), and the mmpbsa.in file.

Also, note that version v1.4.1 is not currently supported. However, the new versions reproduce the energetic results of the previous versions. Therefore, upgrading to the latest version (v1.5.6) would be advisable.

[maral afshinpour]

[maral afshinpour]

On Monday, July 18, 2022 at 1:48:19 PM UTC-5 maral afshinpour wrote:

[maral afshinpour]

[maral afshinpour]

I think I mistakenly send the wrong xtc file. please use the new one.

[Mario Sergio Valdes]

Probably, the error is due to the fact that your system is charged, so in version v1.4.1, you must explicitly define the charge of each molecule, in your case qmcharge_com = -4.0, qmcharge_rec = -5.0, and qmcharge_lig = 1.0. From version v1.5.5 onwards you don't have to worry about this, since gmx_MMPBSA determines these parameters automatically (but you can still change them manually).
We recommend using the latest version. Note that it requires Python 3.9. You can perform the installation using conda as described here, which does not require administrator permissions, so it works both locally and on HPC.

HTH!

Mario S.

[Mario Sergio Valdes]

This is the output of your calculation using exactly the files you sent me. Note that there are several (6) warnings, some of them explicitly related to QM in the definition of the charge of the molecules, since they are not explicitly defined in the input file. If these values are incorrect, you can define it explicitly in the input file using the variables I mentioned above.

Important links:

Installation: https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/installation/#installing-gmx_mmpbsa-v15x

Input file:

     General options: https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#basic-input-options

     QM option: https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#qm-options

[maral afshinpour]

Thanks for the reply, I add this parameter to mmpbs.in file but this time I got error:

 app.read_input_file()
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/main.py", line 755, in read_input_file
    self.INPUT = self.input_file.Parse(infile)
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/input_parser.py", line 445, in Parse
    self.namelists[declared_namelists[i]].variables[key]. \
  File "/gpfs/shared/apps_local/python/3.7-1.20/envs/ambertools/lib/python3.8/site-packages/GMXMMPBSA/input_parser.py", line 113, in SetValue
    self.value = self.datatype(value)
ValueError: invalid literal for int() with base 10: '-4.0'
Exiting. All files have been retained.

[Mario Sergio Valdes]

Oh, in that version it must be an integer, so it should be qmcharge_com = -4, qmcharge_rec = -5, and qmcharge_lig = 1

[maral afshinpour]

Can you please use this index file and cg 10 11. because in command you chose 17 and 18.

Thanks,

Maral

[Mario Sergio Valdes]

I generated a new index that did not have the file. In my case, 17 corresponds to the receptor which in your case should be 11, and 18 corresponds to the ligand which in your case would be 10. You can use any index as long as the groups you select are consistent with the trajectory.

Note that this output corresponds to version v1.5.6, which is much more informative and consistent than version v1.4.1.

[maral afshinpour]

thank you so much, it works but in mmpbsa.in the file, I specified different residues but in the Final result file it shows different residues.

|qm_residues="A/924-942 B/4"
|# qm_residues selection with "within" keyword (new in v1.5.0)
|#qm_residues="within 5"
|/--------------------------------------------------------------
|gmx_MMPBSA Version=v1.4.1 based on MMPBSA.py v.16.0
|Complex topology file:           COM.prmtop
|Receptor topology file:          REC.prmtop
|Ligand topology file:            LIG.prmtop
|Initial mdcrd(s):                COM_traj_0.xtc
|
|Receptor mask:                  ":108-114"
|Ligand mask:                    ":1-107"
|
|Calculations performed using 11 complex frames.
|
|Generalized Born ESURF calculated using 'LCPO' surface areas
|
|All units are reported in kcal/mole.
|QM/MM: Residues 30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,111 are treated with the Quantum Hamiltonian PM3

[Mario Sergio Valdes]

The way you defined it is based on the notation we created that includes the chains, while Amber works with mask starting from residue 1. This means that no matter the residue number in the PDB, amber will always start with 1 and number continuously until it includes all the residues ignoring the chains as well.
According to your output file, you are inversely defining the receptor and ligand. This is not a problem but can be confusing.

[maral afshinpour]

thanks, that makes more sense. so it doesn't need to change the residue in the mmpbsa.in?

[Mario Sergio Valdes]

No, since the selection notation is conditioned by the complex structure itself, not by the way the groups are defined, i.e., the notation does not take into account whether the complex is defined as A+B or B+A. The way you defined QM residues it is not wrong, but it is confusing since your selection of residues is defined as A+B while the complex is defined as B+A. In this case, ideally, the complex should be defined as A+B and it makes more sense because A is much larger than B.

HTH!

Mario S.