Problem when loading TorsionDrive qcmol from QCA

[wutobias]

Hi folks,

please have a look at the following code:

import qcportal as ptl
import cmiles

print(ptl.__version__)
print(cmiles.__version__)

qc_client  = ptl.FractalClient("https://api.qcarchive.molssi.org:443/")
ds         = qc_client.get_collection('TorsionDriveDataset', "OpenFF Substituted Phenyl Set 1")
qcrecord   = ds.get_record("CCC(=O)Nc1[cH:1][c:2](ccn1)[NH:3][CH2:4]C", 
                           specification="default")
qcmols     = qcrecord.get_final_molecules()

for torsion_crd, qcmol in qcmols.items():
    print(qcmol["symbols"].shape)
    print(qcmol["geometry"].shape)
    print(qcmol["symbols"].shape[0]==qcmol["geometry"].shape[0])
    oemol = cmiles.utils.load_molecule(qcmol, toolkit="openeye")

Fails with the following error:

v0.14.0
v0.1.6
(29,)
(29, 3)
True
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-7-f79a5682ac30> in <module>
     16 qcmols     = qcrecord.get_final_molecules()
     17 for torsion_crd, qcmol in qcmols.items():
---> 18     oemol = cmiles.utils.load_molecule(qcmol, toolkit="openeye")

~/progs/anaconda3/envs/atomtyping/lib/python3.7/site-packages/cmiles/utils.py in load_molecule(inp_molecule, toolkit, **kwargs)
     65     if isinstance(inp_molecule, dict):
     66         # This is a JSON molecule.
---> 67         molecule = mol_from_json(inp_molecule, toolkit=toolkit, **kwargs)
     68 
     69     elif isinstance(inp_molecule, str):

~/progs/anaconda3/envs/atomtyping/lib/python3.7/site-packages/cmiles/utils.py in mol_from_json(inp_molecule, toolkit, **kwargs)
    126     geometry = np.asarray(inp_molecule['geometry'], dtype=float)*BOHR_2_ANGSTROM
    127     if len(symbols) != geometry.shape[0]/3:
--> 128         raise ValueError("Number of atoms in molecule does not match length of position array")
    129 
    130     if toolkit == 'openeye' and has_openeye:

ValueError: Number of atoms in molecule does not match length of position array

The test qcmol in the example should be fine and the number of symbols matches the number of coordinates. I think the offending line must be changed to something like:

if symbols.shape[0] != geometry.shape[0]:
    raise ValueError("Number of atoms in molecule does not match length of position array")

Thanks, Tobias

[mattwthompson]

I can reproduce this with a few modifications. I think your suggested fix is good (it's actually somewhat confusing how this was working before).

Could you post the output of conda list? I'm mostly curious what version of Pydantic you have installed.

[wutobias]

Thanks for looking into this. I think this only broke recently, since the geometry in the qcmols seems now to be parsed as (natoms, 3) np.ndarray. I found at least two more lines that may suffer from the same issue (current master branch):

• cmiles/_cmiles_oe.py line 42

  # Add geometry
  if molecule.NumAtoms() != geometry.shape[0]/3:
      raise ValueError("Number of atoms in molecule does not match length of position array")

• cmiles/_cmiles_rd.py line 29

  geometry = geometry.reshape(int(len(geometry)/3), 3)

  My version of Pydantic is 1.8.1, see attached file.

env.txt

[mattwthompson]

Could you look at the changes in #59 and see if those will fix your issues? I think I found the same changes you did 👍

[wutobias]

I can confirm this fixes #55. The fix also works with qcmols from OptimizationDatasets. For instance:

import qcportal as ptl
import cmiles

qc_client  = ptl.FractalClient("https://api.qcarchive.molssi.org:443/")
ds         = qc_client.get_collection('OptimizationDataset', "FDA Optimization Dataset 1")
qcrecord   = ds.get_record("C[NH2+]C[C@@H](c1ccc(c(c1)O)O)O-0", 
                           specification="default")
qcmol      = qcrecord.get_final_molecule()
oemol      = cmiles.utils.load_molecule(qcmol.dict(), toolkit="openeye")
rdmol      = cmiles.utils.load_molecule(qcmol.dict(), toolkit="rdkit")

The oemol is returned successfully, however the rdmol fails due to another reason not related to this issue I think RuntimeError: Could not sanitize molecule. I'll open another issue for this.