How to generate GROMACS format files

[howlfwq]

Hi, I'm trying to use qubekit to get bespoke force filed parameters. I tested it with ethanol. It worked without error. However, I checked the folder of "final parameter". The default format seems to be OpenMM. I created a config file and found the following command: "parametrisation": { "type": "OpenFF", "force_field": "openff_unconstrained-1.3.0.offxml" Is this the correct spot to change the format of derived parameters? If it is, how to change it for GROMACS? I have read the paper and it said that qubekit supports BOSS/MCPRO, OpenMM and GROMACS. But I could not find the instruction. Thanks a lot.

[jthorton]

Hi @howlfwq thanks for the question this is something we defiantly need to document as it will help a lot of users! The best way to convert to GROMACS is via ParmEd using the pdb and XML files in the final_parameters folder, I would create an openmm system and pass this to parmed and then use it to write out the GROMACS files, note that a force field with virtual sites probably will not pass through this conversion correctly. Attached is a code example which should allow you to convert a non-vsite force field to GROAMCS format.

from openmm import app
from parmed.openmm import load_topology

pdbfile = app.PDBFile('mol01.pdb')
force_field = app.ForceField('mol01.xml')

system = force_field.createSystem(pdbfile.topology)
parm_top = load_topology(topology=pdbfile.topology, system=system, xyz=pdbfile.positions)

parm_top.save('molecule.top')
parm_top.save('molecule.gro')

The parametrisation flag you found in the config refers to the initial parameterisation engine used to get reference parameters for the molecule. If we take benzene, for example, we borrow ridged torsions and improper terms from the openff-1.3.0 as QUBEKit does not fit or derive these terms by default. I hope this helps!

[howlfwq]

Hi @jthorton thank you for your help. I have tried using the ParmEd to convert the '.xml' and 'pdb' to '.top' and '.gro' file. But I got this error: Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead. Traceback (most recent call last): File "/public/home/user/qubekit/ETOH/QUBEKit_ETOH_2022_05_07/final_parameters/parm.py", line 11, in parm_top.save('ETOH_convert.gro') File "/public/home/user/anaconda3/envs/qubekit/lib/python3.9/site-packages/parmed/structure.py", line 1486, in save gromacs.GromacsGroFile.write(self, fname, **kwargs) File "/public/home/user/anaconda3/envs/qubekit/lib/python3.9/site-packages/parmed/gromacs/gromacsgro.py", line 303, in write raise RuntimeError("Could not find %s" % atom) RuntimeError: Could not find <Atom C1 [0]; In UNL 0>

I still got the top and gro file. but these two files both have one hydrogen atom missing. I have attached the files and changed the format to txt. ETOH_convert_top.txt ETOH_convert_gro.txt

[jthorton]

Very strange can you share the follwing information which should help me work out whats going on:

1. the qubekit command used to run ethanol
2. the output pdb and xml file
3. the output from conda list Thanks!

[howlfwq]

1. qubekit run -i ETOH.pdb -s non_bonded (If I don't skip non-bonded, qubekit will report and error.
2. I have attached all the output files from qubekit. QUBEKit_ETOH.zip
3. (qubekit) user@***08:~/qubekit/ETOH$ qubekit run -i ETOH.pdb -s non_bonded Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. How ever, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eye sopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free To olkit licenses for academics: https://www.eyesopen.com/academic-licensing Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead. If QUBEKit ever breaks or you would like to view timings and loads of other info , view the log file. Our documentation (README.md) also contains help on handling the various command s for QUBEKit.

Parametrising molecule with openff_unconstrained-1.3.0.offxml. [11:00:09] WARNING: not removing hydrogen atom without neighbors Molecule parameterised and values stored. Optimising the molecule using pre_optimisation method gfn2xtb. [11:00:11] Molecule does not have explicit Hs. Consider calling AddHs() [11:00:11] WARNING: not removing hydrogen atom without neighbors Optimising conformers with gfn2xtb: 100%|█████████| 1/1 [00:02<00:00, 2.38s/it] Partial optimisation complete Optimising molecule using gaussian with wb97x-d/6-311++G(d,p) Optimising conformer: 0%| | 0/1 [02:20<?, ?it/s] Molecule optimisation complete. Calculating Hessian matrix. Symmetry check successful. The matrix is symmetric within an error of 1e-05. Hessian matrix calculated and confirmed to be symmetric. Calculating charges using chargemol and ddec6. Charges calculated and AIM reference data stored. Calculating new bond and angle parameters with the modified Seminario method. Bond and angle parameters calculated. Performing QM-constrained optimisation with Torsiondrive and gaussian No rotatable bonds found to scan! Torsiondrive finished and QM results saved. Performing torsion optimisations using ForceBalance. Torsion optimisation complete.

[jthorton]

Thanks for sharing the QUBEKit folder @howlfwq that helped identify the issue, it seems that the input pdb file is missing a bond between the oxygen and hydrogen which is causing the warnings that are printed to screen about one hydrogen having no neighbours etc. Please can you try running it again starting from smiles or from a pdb/sdf with full explicit connectivity and let me know how you get on? I have opened an issue for this in #315 to try and catch this early with an informative error which should help other users in future!

[howlfwq]

I have tried ethanol with connectivity. I got correct parameter files and Parmed converted them successfully. Thank you again!