first of all - thank you so much for the excellent tutorial, as well as the biobb package. it really helped me get set up. I am wondering if its possible to extend this to systems with multiple ligands/cofactors and metals? Do you have a tutorial about it, or how would i go about implementing that myself?
thank you so much.
Hi David,
unfortunately, the addLigand building block is only working with one ligand. If more than one ligand is present in your structure, a further modification of the topology file needs to be performed, and this is not implemented in the current version. You can try to follow a tutorial on how to manually modify the gromacs topology for that (e.g. Justin Lemkul’s tutorial).
Something similar happens with cofactors and metals, which usually require coordination with protein atoms. In this case, you can manually modify the special bonds file in GROMACS (specbond.dat, see the manual pages), adding the required bonds. In this case, you can still use the BioBB library, wrapping a local GROMACS installation, where you can manually modify the specbond.dat file. This is what we do for example when dealing with structural Zinc ions.
Hope this helps!
thank you for your response:) i am quite the novice when it comes to MD. I guess i will just do things manually in those cases then. May i ask you one more question? The add hydrogen script sometimes does not add any/all hydrogens. I wonder if that is an error or supposed to be like that? from my understanding all hydrogens have to be described. thank you
Hi David,
can I ask you which building block are you using to add hydrogens to your system? And what type of system are you dealing with? The library has different building blocks to add hydrogens, wrapping different tools. Depending on the system you are working with, the “proper” building block to add hydrogens can change.
i am following this tutorial:
https://mmb.irbbarcelona.org/biobb/availability/tutorials/protein-complex_md_setup
so for example on phenol it might only add hydrogens in the ortho/para positions. what i do now is i just add hydrogens with avogadro- that works and then i continue as in the tutorial
Hello,
could you please clarify how to use a local GROMACS installation with the BioBB library? I am working on a similar case and need to modify the specbond.dat file and work with a local GROMACS install.
Thank you.
Hi Juan,
every building block from the biobb_md module accepts a couple of properties to indicate the wrapper where to find the gromacs binary and the associated force-fields:
See for example the pdb2gmx documentation.
In principle, you should add these two properties in all the gromacs building blocks, with the corresponding paths, gmx_path to the local gromacs binary, and gmx_lib to the local gromacs force-field folder, where your modified specbond.dat file is.
Hope it helps!
Thanks for using our tools!
-Adam-
Thank you Adam!
I did this:
from biobb_md.gromacs.pdb2gmx import Pdb2gmx
output_pdb2gmx_gro = pdbCode+’_pdb2gmx.gro’
output_pdb2gmx_top_zip = pdbCode+’_pdb2gmx_top.zip’
prop = {
‘force_field’ : ‘amber99sb’,
‘water_type’: ‘spce’,
‘gmxlib’ :’/Applications/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/’,
‘gmx_path’:’/Applications/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin/gmx’
}
Pdb2gmx(input_pdb_path=fixed_pdb,
output_gro_path=output_pdb2gmx_gro,
output_top_zip_path=output_pdb2gmx_top_zip,
properties=prop).launch()
and got this error:
Inconsistency in user input:
Force field ‘amber99sb’ occurs in 2 places, but not in the current directory.
Run without the -ff switch and select the force field interactively.
I don´t understand what´s the problem as amber99sb.ff is in the gmxlib path and here is also the modified specbond.dat file.
Any clues?
Thank you!
Juan
Hi Juan,
could you please try it without the gmxlib property? Gromacs is finding more than one force field named amber99sb and doesn’t know which one should be selected. That usually happens when there’s more than one path to gromacs force fields, and gmxlib is adding these paths.
Waiting for your reply!
-Adam-
Hello Adam,
sorry for the delay in my reply.
Somehow it finally worked without errors. Anyway I wonder why would the removal of the gmxlib property be a reasonable approach. If I remove gmxlib I would not be using the modified force field including the Zn parameters that I need right?
On a different note, despite including positional restraints on my Zn ion I can see it moving outside its coordination site during the NVT equilibration, so this is quite puzzling…
Thank you
Regards
Juan
Hi Juan,
sorry I didn’t realise you needed the gmx_lib for the specbonds, I forgot about your previous e-mail, my apologies.
In this case we should probably remove all possible paths other than the one assigned in the gmx_lib variable in the set of available environment variables, so that gromacs is only finding one force field folder. This problem is easy to solve when working in an interactive way, as you can use a number as an input for the force-field selector, but is not the case for the -ff parameter that is the one used by the building block. Just out of curiosity, how did you solve it?
Regarding the Zn coordination, the behaviour of the special bonds is not always easy to follow. I would recommend to take a careful look at the log file of the pdb2gmx tool, and identify if the restraint is applied. A line like this one should appear (similar to the ones appearing for the disulfide bonds):
Linking HID-356 NE2-2891 and ZN2-597 ZN-4863...
If this is not present, then you should try to play with the distance assigned in the specbond file, which if I’m not wrong gromacs is using to try to find an atom within ±10%
of this value, otherwise the bond is not created.
Hope it helps.
Regards,
-Adam-
Hi Adam,
thank you. I did it as you said calling pdb2gmx of grommacs from the command line and I don´t get the line that you mention but the following (also same output when I do it with the biobb_md.gromacs.pdb2gmx)
"Residue ZN443 has type ‘Ion’, assuming it is not linked into a chain.
Problem with chain definition, or missing terminal residues. This chain does not appear to contain a recognized chain molecule. If this is incorrect, you can edit residuetypes.dat to modify the behavior."
The definition of the Zn ion in the pdb file is HETATM and it´s separeted from the protein entries and in residuetypes.dat is Ion. Also, in the specbond.dat I have added:
CYZ SG 1 ZN ZN 4 2.52 CYZ ZN
2.52 ±10% covers the distance range for the 4 cysteine residues coordinating Zn in my pdb file. CYZ is a new residue type that I´ve also included in the residuetypes.data file. So I don´t truly understand the meaning of this error message.
Sorry for so many questions.
Kind regards
Juan
Hi Juan,
we are dangerously entering into the gromacs field 
, but maybe you can use the pdb2gmx -merge parameter, to ensure gromacs is checking distances between any pair of atoms/residues, including those in different chains.
-Adam-
Thanks a lot for all your comments Adam! I´ll give it a try and will talk to GROMACS people if it doesn´t work.
Cheers
Juan