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 66 forks source link

AMBER Topology from GROMACS simulation with modified forcefiel #287

Closed wulfbear closed 2 years ago

wulfbear commented 2 years ago

My Question is...

Dear all, for my gromacs simulation of a polymer and a protein I have added some modified forcefield parameters for the polymer atoms into the standard AMBER99SBildn forcefield. So far, everything has worked out well in the GROMACS universe.

I now came across your promising tool and would like to use my simulational data with the non-standard atom and residue definition of my "GROMACS" ff with gmx_MMPBSA. I have already tried to include all atomtype, bond, dihedral and angle definitions in the respective frcmod.ff99SBildn and parm99.dat files (the mmpbsa.in file I used defines these two as forcefields). However, once I run the following command: gmx_MMPBSA -O -i mmpbsa_Protein_Ligand_ILDN_GAFF.in -cs complex_bound.pdb -cp topol.top -ci complex_index_mmpbsa.ndx -cg 1 12 -ct complex_bound_traj.pdb -rs rec_protein.pdb -ri rec_protein_index.ndx -rg 10 -rt rec_protein_traj.pdb -ls lig_polymer.pdb -li lig_polymer_index.ndx -lg 5 -lt lig_polymer_traj.pdb -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv

gmx_MMPBSA throws an error upon building the normal complex AMBER Topology:

gmx_MMPBSA -O -i mmpbsa_Protein_Ligand_ILDN_GAFF.in -cs complex_bound.pdb -cp topol.top -ci complex_index_mmpbsa.ndx -cg 1 12 -ct complex_bound_traj.pdb -rs rec_protein.pdb -ri rec_protein_index.ndx -rg 10 -rt rec_protein_traj.pdb -ls lig_polymer.pdb -li lig_polymer_index.ndx -lg 5 -lt lig_polymer_traj.pdb -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv [ INFO ] Starting gmx_MMPBSA v1.5.6
[INFO ] Command-line gmx_MMPBSA -O -i mmpbsa_Protein_Ligand_ILDN_GAFF.in -cs complex_bound.pdb -cp topol.top -ci complex_index_mmpbsa.ndx -cg 1 12 -ct complex_bound_traj.pdb -rs rec_protein.pdb -ri rec_protein_index.ndx -rg 10 -rt rec_protein_traj.pdb -ls lig_polymer.pdb -li lig_polymer_index.ndx -lg 5 -lt lig_polymer_traj.pdb -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv [INFO ] Checking mmpbsa_Protein_Ligand_ILDN_GAFF.in input file...
[WARNING] The startframe variable must be >= 1. Changing startframe from -2 to 1
[INFO ] Checking mmpbsa_Protein_Ligand_ILDN_GAFF.in input file...Done. [INFO ] Checking external programs...
[INFO ] cpptraj found! Using /home/kai/anaconda3/envs/gmxMMPBSA/bin/cpptraj
[INFO ] tleap found! Using /home/kai/anaconda3/envs/gmxMMPBSA/bin/tleap
[INFO ] parmchk2 found! Using /home/kai/anaconda3/envs/gmxMMPBSA/bin/parmchk2
[INFO ] sander found! Using /home/kai/anaconda3/envs/gmxMMPBSA/bin/sander
[INFO ] Using GROMACS version > 5.x.x!
[INFO ] gmx_mpi found! Using /usr/local/bin/gmx_mpi
[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...
[INFO ] Normal Complex: Saving group 1_12 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb [INFO ] A receptor structure file was defined. Using MT approach...
[INFO ] Normal Receptor: Saving group 10 in rec_protein_index.ndx file as _GMXMMPBSA_REC.pdb
[INFO ] A ligand structure file was defined. Using MT approach...
[INFO ] Normal Ligand: Saving group 5 in lig_polymer_index.ndx file as _GMXMMPBSA_LIG.pdb
[INFO ] Checking the structures consistency...
[WARNING] No reference structure was found and the complex structure not contain any chain ID. Assigning chains ID automatically...
[INFO ]
[INFO ] Using topology conversion. Setting radiopt = 0...
[INFO ] Building Normal Complex Amber topology... File "/home/kai/anaconda3/envs/gmxMMPBSA/bin/gmx_MMPBSA", line 8, in sys.exit(gmxmmpbsa())
File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 101, in gmxmmpbsa app.make_prmtops()
File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 537, in make_prmtops self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology() File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 118, in buildTopology tops = self.gmxtop2prmtop() File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 461, in gmxtop2prmtop com_top = self.cleantop(self.FILES.complex_top, self.indexes['COM']['COM'])
File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 748, in cleantop rtemp_top = parmed.gromacs.GromacsTopologyFile(ttp_file.as_posix())
File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py", line 251, in init self.read(fname, defines, parametrize) File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py", line 461, in read self.parametrize() File "/home/kai/anaconda3/envs/gmxMMPBSA/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py", line 966, in parametrize atom.atom_type = params.atom_types[atom.type] KeyError: 'c3' Exiting. All files have been retained.

It seems like the modified atomtype of my polymer residues, e.g. c3, is causing issues during the topology build. Did I miss anything when modifying the ff used by gmx_MMPBSA or is this not done as easily in AMBER as I am used to in GROMACS?

Thank you all in advance for any advice or help that you can offer.

Warmest regards,

Kai.

marioernestovaldes commented 2 years ago

Could you send us the files to check what's going wrong!?

wulfbear commented 2 years ago

I have attached the files to this post via Link:

https://1drv.ms/u/s!AlxGtPgNe0LzhLIahdlTCq8dyjaE2w?e=WabD4C

Thank you!

marioernestovaldes commented 2 years ago

I see you are using a modified version of the force filed... could you please send us the folder with the force field files?

marioernestovaldes commented 2 years ago

please, send also the "CgDef_K10L.itp" and "PLA.itp" files

wulfbear commented 2 years ago

ff and itp.zip Hi, included both itp files and the $AMBERHOME/dat/leap/cmd and /dat/leap/parm folders Please let me know if you need anything else.

Thanks for your time!

marioernestovaldes commented 2 years ago

the [ atomtypes ] directive is missing in the PLA.itp file... those atoms are not in the force field, hence an atomtype definition is needed... how did you generate this .itp file for PLA? in this case, the amber files are not needed since gmx_MMPBSA takes care of that for you

wulfbear commented 2 years ago

The topology was created using gromacs pdb2gmx functionality. This gave the PLA.top file which was manually edited into PLA.itp (only removing the #include statements for forcefield, water / ion topology and position restraints, [system] and [molecules] definition). In the original PLA.top I can also not find the [ atomtypes ] directive, only the column "type" in the [atoms] directive:

; This is a standalone topology file ; ; Created by: ; :-) GROMACS - gmx pdb2gmx, 2020.1-Ubuntu-2020.1-1 (-: ;
; Executable: /usr/bin/gmx ; Data prefix: /usr ; Working dir: /mnt/f/01_PhD Bioinformatics Backup/02_Polymer/05_Gusava_paper/03_helix_ordered ; Command line: ; gmx pdb2gmx -f helix_layer.gro -o helix_layer.gro -i ; Force field data was read from: ; /mnt/c/Users/Kai/Documents/03_PhDBioinformatics/02_Polymer/05_Gusava_paper/top ; ; Note: ; This might be a non-standard force field location. When you use this topology, the ; force field must either be present in the current directory, or the location ; specified in the GMXLIB path variable or with the 'include' mdp file option. ; ; Include forcefield parameters

include "amber99sb-ildn.ff/forcefield.itp"

[ moleculetype ] ; Name nrexcl Polymer 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB ; residue 1 PLAO rtp PLAO q -0.2 1 c3 1 PLAO C3 1 0.110767 12.01 2 o 1 PLAO O4 2 -0.535335 16 3 h1 1 PLAO H5 3 0.111367 1.008 4 oh 1 PLAO O10 4 -0.578135 16 5 ho 1 PLAO H11 5 0.419667 1.008 6 c3 1 PLAO C6 6 -0.131433 12.01 7 hc 1 PLAO H7 7 0.0617 1.008 8 hc 1 PLAO H8 8 0.0617 1.008 9 hc 1 PLAO H9 9 0.0617 1.008 10 os 1 PLAO O1 10 -0.418234 16

Where does gmx_MMPBSA expect the [ atomtypes ] directive?

marioernestovaldes commented 2 years ago

I see... so, you did use the standard pdb2gmx in GROMACS to generate the topology file for this molecule? did you modify something in the force field folder "mber99sb-ildn.ff"

these atoms are not present by defult in the force field:

c3 o h1 oh ho c3 hc hc hc os c

that's why the [ atomtypes ] directive is needed in the topology, unless you modified the force field and then, I will need that folder with the modified files

wulfbear commented 2 years ago

Hi, yes I modified the amber99sb-ildn.ff and included all non-present atom types. Heres the modified ff: modified amber99sb-ildn.zip

Sorry, I though you meant the modified ff in the Amber directories previously.

marioernestovaldes commented 2 years ago

place the "amber99sb-ildn.ff" folder with the modified files in the same folder you are running gmx_MMPBSA and everything should run smoothly... I did notice your trajectory is not fit properly... please, remove PBC conditions and fit the trajectory, otherwise you will encounter an error like this one:

'ValueError: could not convert string to float: '*****''

you can check this page (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/Q%26A/calculations/#possible-solutions) for more help...

please, let me know if you have any other questions