bespokefit 0.4.3 failed with "forming a nearly linear angle" for the compound working fine with bespokefit 0.2.3

[okbckim]

Description

bespokefit job failed for the structure attached when the version 0.4.3 is used. This job worked fine with the version 0.2.3

Error message:

 [null, null, {"type": "ValueError", "message": "TorsionDrive error at
 60:\ngeomeTRIC run_json error:\nTraceback (most recent call last):\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/run_json.py\", line 269, in geometric_run_json\n
 progress = geometric.optimize.Optimize(coords, M, IC, engine, dirname, params,
 print_info = (ic==0))\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1116, in Optimize\n
 return optimizer.optimizeGeometry()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1029, in optimizeGeometry\n
 self.calcEnergyForce()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 306, in calcEnergyForce\n
 self.checkStructure()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1067, in checkStructure\n
 raise LinearTorsionError(\"A constrained torsion has three consecutive
 atoms\\n\"\ngeometric.errors.LinearTorsionError: A constrained torsion has
 three consecutive atoms\nforming a nearly linear angle, making the torsion
 angle poorly defined.\n> Atoms Angle\n> 4-5-6 179.99\n\n", "traceback":
 "Traceback (most recent call last):\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/celery/app/trace.py\", line 453, in trace_task\n    R =
 retval = fun(*args, **kwargs)\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/celery/app/trace.py\", line 736, in __protected_call__\n
 return self.run(*args, **kwargs)\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/openff/bespokefit/executor/services/qcgenerator/worker.py\
 ", line 151, in compute_torsion_drive\n    return_value =
 qcengine.compute_procedure(\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/qcengine/compute.py\", line 181, in compute_procedure\n
 return handle_output_metadata(output_data, metadata, raise_error=raise_error,
 return_dict=return_dict)\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/qcengine/util.py\", line 186, in handle_output_metadata\n
 raise ValueError(output_fusion[\"error\"][\"error_message\"])\nValueError:
 TorsionDrive error at 60:\ngeomeTRIC run_json error:\nTraceback (most recent
 call last):\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/run_json.py\", line 269, in geometric_run_json\n
 progress = geometric.optimize.Optimize(coords, M, IC, engine, dirname, params,
 print_info = (ic==0))\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1116, in Optimize\n
 return optimizer.optimizeGeometry()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1029, in optimizeGeometry\n
 self.calcEnergyForce()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 306, in calcEnergyForce\n
 self.checkStructure()\n  File
 \"/home/byungchan_kim_neurolambdatx_com/miniforge3/envs/bespokefit-env/lib/pyth
 on3.10/site-packages/geometric/optimize.py\", line 1067, in checkStructure\n
 raise LinearTorsionError(\"A constrained torsion has three consecutive
 atoms\\n\"\ngeometric.errors.LinearTorsionError: A constrained torsion has
 three consecutive atoms\nforming a nearly linear angle, making the torsion
 angle poorly defined.\n> Atoms Angle\n> 4-5-6 179.99\n\n\n"}, null, null, null,
 null, null, null, null]

Please describe what you were trying to do, what you expected to happen, and what happened instead.
structure.offxml

Reproduction

Please include a minimally reproducing example of this bug.

Output

Please include the output, including full tracebacks of any errors, resulting from the reproduction.

Software versions

openff-bespokefit         0.4.3              pyhd8ed1ab_0    conda-forge

• Which operating system and version are you using? Linux
• How did you install BespokeFit?

mamba create -n bespokefit-env -y -c conda-forge  python=3.10 openff-bespokefit "qcportal <0.50" psi4 dftd3-python ambertools cuda-version=12.4

• Are you using Apple Silicon? If so, are you running BespokeFit in Rosetta (osx-64) or natively (osx-arm64)? No
• What is the output of running conda list?

Output of conda list

Please place the output of conda list here
openff-bespokefit         0.4.3              pyhd8ed1ab_0    conda-forge
forcebalance              1.9.6           py310h256d06d_2    conda-forge
openff-forcefields        2024.09.0          pyhff2d567_0    conda-forge
xtb                       6.7.1                h8876d29_4    conda-forge
xtb-python                22.1            py310ha75aee5_2    conda-forge
openff-amber-ff-ports     0.0.4              pyhca7485f_0    conda-forge
openff-bespokefit         0.4.3              pyhd8ed1ab_0    conda-forge
openff-forcefields        2024.09.0          pyhff2d567_0    conda-forge
openff-fragmenter-base    0.2.2              pyhd8ed1ab_0    conda-forge
openff-interchange        0.3.29             pyhd8ed1ab_0    conda-forge
openff-interchange-base   0.3.29             pyhd8ed1ab_0    conda-forge
openff-models             0.1.2              pyhca7485f_0    conda-forge
openff-qcsubmit           0.50.0             pyhd8ed1ab_0    conda-forge
openff-toolkit            0.16.5             pyhd8ed1ab_0    conda-forge
openff-toolkit-base       0.16.5             pyhd8ed1ab_0    conda-forge
openff-units              0.2.2              pyhd8ed1ab_1    conda-forge
openff-utilities          0.1.15             pyhd8ed1ab_0    conda-forge

023.tgz
043.tgz

[j-wags]

Summary

Hi @okbckim, thanks for writing in. I've been unable to reproduce this problem using xTB-GFN2, a fast semiempirical QM method (this molecule runs successfully in bespokefit with xTB using both openff-bespokefit 0.2.3 and 0.4.3). My top suspicion is that there's a molecular rearrangement happening somewhere in the psi4 optimizations (that doesn't happen in xTB) that causes a torsion to become linear. If this is the case, it's unfortunately out of scope for us to fix, since

1. the QM knows nothing about bonds, and it would be really complicated to add constraints to the optimizations in the QM torsion drive (also it would no longer really be the requested QM method/it would be a big behavior change), and
2. throwing out fragments and continuing with the fit when QM rearrangements happen would also be extremely complicated (one torsiondrive may affect several torsion parameters, since some fragments have multiple torsions)

It also seems weird to me that this torsion rearrangement happens in one version of bespokefit and not another, since the geometry rearrangement would be entirely within psi4/geometric/torsiondrive. Maybe there's been a fine numerical change in one of those packages, or some of the QM geometries are near a saddle point and some small changes cause them to go the wrong way. It would be really expensive to verify this, and even if we did it wouldn't point us toward a course of action that we could afford to follow.

Details

My initial hypothesis was that something changed between openff-bespokefit 0.2.3 and 0.4.3 that cause this linear nitrile to be considered rotatable.

It would be great if this were the bug, since it'll be easy to debug+fix this problem in cheminformatics-land.

I tested this hypothesis by creating envs with openff-bespokefit 0.2.3 and 0.4.3 and xtb-python and trying to run this molecule. If the problem is that the linear nitrile is being considered rotatable, the QC method shouldn't matter.

openff-bespokefit 0.2.3

The 0.2.3 env was created by running:

micromamba env create -n temp-bf-023 -c conda-forge xtb-python openff-bespokefit=0.2.3 python=3.10
micromamba activate temp-bf-023
micromamba install -c conda-forge ambertools=23
micromamba install -c conda-forge rdkit=2023

The output of running micromamba env export in this environment is 023.txt

and bespokefit ran using (don't read too far into acetaminophen, I just copied that from the quickstart guide):

openff-bespoke executor run --file               18660-1.sdf    \
                            --workflow           "default"              \
                            --output             "acetaminophen.json"   \
                            --output-force-field "acetaminophen.offxml" \
                            --n-qc-compute-workers 2                    \
                            --qc-compute-n-cores   1                    \
                            --default-qc-spec xtb gfn2xtb none

This succeeded.

openff-bespokefit 0.4.3

The 0.4.3 env was created by running:

micromamba create -n temp-bf-043 -c conda-forge openff-bespokefit xtb-python
micromamba activate temp-bf-043

The output of running micromamba env export in this environment is 043.txt

and bespokefit run using:

BEFLOW_KEEP_TMP_FILES=True BEFLOW_OPTIMIZER_KEEP_FILES=True    \
                            openff-bespoke executor run --file               18660-1.sdf    \
                            --workflow           "default"              \
                            --output             "acetaminophen.json"   \
                            --output-force-field "acetaminophen.offxml" \
                            --n-qc-compute-workers 2                    \
                            --qc-compute-n-cores   1                    \
                            --default-qc-spec xtb gfn2xtb none --directory bespoke-tmp

This also succeeded.

Since changing the version of bespokefit didn't reproduce the bug, my next best hypothesis is that this is something deeper in QM-land, and we unfortunately can't afford to delve deeper on it. It's good that some error was emitted during the scan and bespokefit didn't try to fit a FF to the mangled molecule, however I'd interpret this as meaning that there's a stochastic error with this molecule, and that when bespokefit finishes and produces an OFFXML it's safe to use, but if you don't have time to run multiple times you should just use general Sage parameters instead.