QUBEKit fails with linear molecules

[UnWaDo]

Assume you want to parametrise Br₂ (for clathrates calculations or whatever else). QUBEKit fails at stage bonded_parameters.

Reproducing the problem

qubekit config create -p 5b
qubekit run -c test.conf -sm BrBr -m 1

Error traceback

``` Traceback (most recent call last): File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/workflow/workflow.py", line 335, in _run_stage result_mol = stage.run( File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/bonded/mod_seminario.py", line 231, in run self._modified_seminario_method(molecule=molecule, hessian=hessian) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/bonded/mod_seminario.py", line 264, in _modified_seminario_method self.calculate_angles(eigenvals, eigenvecs, molecule, bond_lens) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/bonded/mod_seminario.py", line 283, in calculate_angles for count, angle in enumerate(molecule.angles): TypeError: 'NoneType' object is not iterable The above exception was the direct cause of the following exception: Traceback (most recent call last): File "/home/unwado/.local/miniconda3/envs/qubekit/bin/qubekit", line 8, in sys.exit(cli()) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/click/core.py", line 1130, in __call__ return self.main(*args, **kwargs) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/click/core.py", line 1055, in main rv = self.invoke(ctx) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/click/core.py", line 1657, in invoke return _process_result(sub_ctx.command.invoke(sub_ctx)) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/click/core.py", line 1404, in invoke return ctx.invoke(self.callback, **ctx.params) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/click/core.py", line 760, in invoke return __callback(*args, **kwargs) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/cli/run.py", line 135, in run workflow.new_workflow(molecule=molecule, skip_stages=skip_stages, end=end) File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/workflow/workflow.py", line 257, in new_workflow return self._run_workflow( File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/workflow/workflow.py", line 295, in _run_workflow molecule = self._run_stage( File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/workflow/workflow.py", line 359, in _run_stage raise WorkFlowExecutionError( qubekit.utils.exceptions.WorkFlowExecutionError: The workflow stopped unexpectedly due to the following error at stage: bonded_parameters ```

What to do?

Main part of traceback (in my opinion)

``` File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/bonded/mod_seminario.py", line 283, in calculate_angles for count, angle in enumerate(molecule.angles): TypeError: 'NoneType' object is not iterable ```

According to spolier above, the main problem is in mod_seminario.py file, where molecule.angles returns None if there are no angles in molecule. I'm not sure if it is critical somewhere, but may be Molecule.angles should return empty list instead of None. If it is important for some check, probably the condition should be changed from molecule.angles is not None to len(molecule.angles) != 0. I can implement changes myself, but would like to hear your opinion on the problem

NB

HHal molecules (Hal = F, Cl, Br, I) fail on virtual_sites stage due to error below, but I have no idea why are they different from Hal₂

Error with HHal

``` Traceback (most recent call last): File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/workflow/workflow.py", line 335, in _run_stage result_mol = stage.run( File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/nonbonded/virtual_sites.py", line 153, in run self._save_virtual_sites() File "/home/unwado/.local/miniconda3/envs/qubekit/lib/python3.10/site-packages/qubekit/nonbonded/virtual_sites.py", line 1180, in _save_virtual_sites close_b_coords = self._coords[closest_atoms[1]] IndexError: list index out of range ```

[UnWaDo]

In case of HHal (and I₂ also) the problem is caused by the fact, that while saving virtual site, the program requires at least three atoms index, which is impossible in case of diatomic molecule.

Possible ways to solve the problem is either to find a way to introduce virtual site to OpenMM's forcefield based on only two atoms or to ignore virtual site for diatomic molecules.

I created an issue in OpenMM repository with the question about virtual sites in diatomic molecules.

[jthorton]

Hi, @UnWaDo thanks for the excellent error report! The angle part is defiantly a bug which I will fix today this was definitely an oversite on our part as we did not envision this use case. With regards to the virtual sites we have tried to use a single type known as the LocalCoordinateSite which is based on three atoms as you point out to make it easier to write out our OpenMM input files. I think a TwoParticleAverageSite in openmm should do the trick here but this will take a little longer to work out, for now, I recommend running without them by choosing a protocol optimised without sites.

[UnWaDo]

Hello, @jthorton! Thank you for your response!

I made some modifications by myself (made .angles to return empty list instead of None and _save_virtual_sites to do nothing in case of diatomic molecules) and it worked fine for me, but surprisingly enough it gave charges for iodine equal to −1.38369 and +1.38369. Other halogens are fine, so it is kind of unexpected for me

Not sure if it is a bug, but may be you should check it after implementing changes, that you noted

[jthorton]

Hi @UnWaDo, I have now made those changes and after some testing, I have found that the H-Hal molecules have an initial ESP error on the hal atom higher than our threshold which is why QUBEKit tries to fit vsites and causes an error (Ill fix this in another PR). Whereas the hal-hal molecules seem to not need the vsites.

For I-I I am not able to run this as the basis set we use in protocol 5b has no parameters for iodine, can you share the QM method/basis and QM program you used to get those charges so I can check what happened?

[UnWaDo]

Hi, @jthorton, I have a new problem: it seems to me, that new fragmentation method with the use of MMFF is not working for HHal molecule

After mmff_properties = AllChem.MMFFGetMoleculeProperties is called, the resulting mmff_properties is null. Setting verbosity level of MMFF to 2, I was managed to find that the problem is with MMFF inability to understand type of hydrogen atom (see spoiler). The problem persists with all HHal molecules

Message from AllChem.MMFFGetMoleculeProperties(mmffVerbosity=2)

``` A T O M T Y P E S A N D C H A R G E S ATOM FORMAL PARTIAL ATOM TYPE CHARGE CHARGE ----------------------------------- Br #1 13 0.000 0.000 H #2 0 0.000 0.000 Missing atom types - charges were not computed ```

Can something be done with it from QUBEKit's side? Is there a way to use another fragmentation method in the latest version?

For I−I I wasn't able to reproduce early mentioned error in modern version of QUBEKit