I am planning on docking a ligand in a cavity with a heme and a constrained oxygen. Whereas I can make params of ligands and non-canonical residues, I cannot for molecular oxygen (OXY). It is not in the database either. Molfile_to_params.py gives the following
$ python molfile_to_params.py oxy.sdf -n OXY Centering ligands at ( 0.000, 0.000, 0.000) Atom names contain duplications -- renaming all atoms. WARNING: structure contains double bonds but no aromatic bonds Aromatic bonds must be identified explicitly -- alternating single/double bonds (Kekule structure) won't cut it. This warning does not apply to you if your molecule really isn't aromatic. Total naive charge -1.520, desired charge 0.000, offsetting all atoms by 0.760 Fragment 1: [' O1 ', ' O2 '] Traceback (most recent call last): File "molfile_to_params.py", line 1323, in <module> sys.exit(main(sys.argv[1:])) File "molfile_to_params.py", line 1295, in main num_frags = fragment_ligand(m) File "molfile_to_params.py", line 535, in fragment_ligand raise ValueError("Fragment %i has %i atoms; merge with another fragment or add virtual atoms to make 3 total" % (frag_id, num_atoms)) ValueError: Fragment 1 has 2 atoms; merge with another fragment or add virtual atoms to make 3 total
The molecular oxygen is from https://pubchem.ncbi.nlm.nih.gov/compound/977#section=3D-Conformer although there is a structure of oxy-myoglobin in the PDB, but the latter has only a single CONECT bonding*.
What is a `virtual atom` that the py-script is requesting?
* I know that molecular oxygen actually resonates between O=O and O.-O. making it blue and paramagnetic. But for now, I am happy with a sloppy double bonded oxygen.
EDIT. I forgot to mention that running `relax` crashes due to the `OXY` so its not handled implicitly. Running the same on a structure with `OXY` removed but with other ligands (`-extra_res_fa HEM.params -extra_res_fa LG.params`) works fine.
You have two atoms. You need 3 or more. Each residue needs to be able to define a local coordinate frame, which requires 3 atoms for 3 DOFs for 3 dimensions. It's a consequence of Rosetta using internal coordinates (length-angle-torsion) for most of its work; you need at least 3 atoms per residue to make that work.
I don't know if there are any diatomic types already in the database, but I know there are monatomic metals. If you look at them you will find they fill out their parameters with "virtual" atoms - they're ghosts w/r/t phyiscs, no effects, but they provide the extra handles the math requires. I don't know if there is a way to do this automatically, but that's where I'd start. (Since your system is two atoms, only the both length matters...you can probably just hack up the ZN params into O=O params if you need to? That's where I'd start.)
molfile_to_params.py recognizes "V" atoms as virtual ones. You should be able to add a spurious V atom to your input sdf/mol2 file, connect it to the oxygens with a pseudo-bond (a single bond will work), and molfile_to_params.py should be able to go from there.
Otherwise, hacking one of the monoatomic metal types probably would work too.
Thanks! I did as R Moretti recomended —adding a vanadium bounded to one of the oxygen— and doublechecked against zinc as suggested. The params file works and does behave weirdly.
Hopefully, I will not have to work with an organovanadate and stumble across the inverse issue!
Here is the OXY.params (O1=O2 molecule) for any future reader.