Multi-Molecule Torsiondrives

[jthorton]

Description

This PR adds support for multi molecule inputs to torsion drives bringing the QCFractal and QCEngine behaviour in sync.

Changelog description

The torsiondrive procedure can use multiple molecules as input to torsiondrives.

Status

• [x] Code base linted
• [x] Ready to go

[codecov[bot]]

Codecov Report

Merging #325 (96bfd3e) into master (0612e81) will not change coverage. The diff coverage is 100.00%.

[jthorton]

requires https://github.com/MolSSI/QCElemental/pull/281 for the tests to pass.

[loriab]

I'm going to merge the qcel 0.24.0 feedstock PR shortly, so the conda package ought to be available soon, so you can make this PR pass.

I'd like to get the next qcng release with this and native_files and probably nwc restart minted soon so I can do a psi4 release next week. please lmk if that timeframe doesn't work for you wrt this PR.

[xiki-tempula]

@loriab Thanks for implementing the Torsiondrives. I'm interested in running a dihedral scan of a certain dihedral of a molecule with GFN2-xTB. It seems that the documentation for the Torsiondrives is not there yet. I wonder if you could give me some help in doing this, please? Thank you.

Suppose I have a molecule with the coordinate

0 1
H           -0.900956416662    -0.508511199761    -0.767347451313
O           -0.728052921632     0.024961188945     0.023986108954
O            0.727620291227     0.033161864577    -0.026966996721
H            0.907822575145    -0.413943838488     0.814656303345

The dihedral to be scanned is 1 2 3 4 with 1-based indexing and I would like to have a step size of 15 degree. The output would be dihedral_list which is a list with a length of 24 indicating the corresponding dihedral angle. energy_list which is a list with a length of 24 indicating the corresponding energy. coord_list which is a list with a length of 24, where each element is an array with the shape of (4,3) indicating the coordinate.

Thank you.

[pavankum]

@xiki-tempula the example in there might be of use to you https://github.com/MolSSI/QCEngine/pull/305#issue-936896576

[xiki-tempula]

@pavankum Thanks for the reply. The comment in https://github.com/MolSSI/QCEngine/pull/305#issue-936896576 would not run for qcengine 'v0.21.0' and qcelemental 'v0.24.0'. With error

    input_data=TorsionDriveInput(
  File "pydantic/main.py", line 331, in pydantic.main.BaseModel.__init__
pydantic.error_wrappers.ValidationError: 1 validation error for TorsionDriveInput
initial_molecule
  value is not a valid list (type=type_error.list)

[jthorton]

Hi, @xiki-tempula The torsiondrive interface has since been updated to take a list of input molecule geometries to seed the calculation which helps with the wavefront propagation try just using a list initial_molecule=[molecule, ].

[xiki-tempula]

@jthorton Thank you very much for the help. I have changed the initial molecule in the to https://github.com/MolSSI/QCEngine/pull/305#issue-936896576

molecule = Molecule(
    symbols=["C", "C", "H", "H", "H", "H", "H", "H"],
    connectivity=[
        (0, 1, 1),
        (0, 2, 1),
        (0, 3, 1),
        (0, 4, 1),
        (1, 5, 1),
        (1, 6, 1),
        (1, 7, 1),
    ],
    geometry=[
        1.5403406836914142, -1.0173082391323545, 0.9312810207342496,
        4.0719763300123235, -0.0975682592642413, -0.02203578938791481,
        0.0002563605701713713, 0.0013953403968729996, 0.0011121160323317373,
        1.309831306165048, -3.03614919350581, 0.5491856718564878,
        1.3800394103640543, -0.7181256543708329, 2.9707878359388222,
        5.6120991748009645, -1.1161249890160694, 0.907991575289461,
        4.302418801484787, 1.921022388748466, 0.36057345099334853,
        4.232223312568672, -0.3961916040297606, -2.061588178357897]
)

result = compute_procedure(
    input_data=TorsionDriveInput(
        keywords={
            "dihedrals": [(2, 0, 1, 5)],
            "grid_spacing": [60]
        },
        input_specification=QCInputSpecification(
            driver=DriverEnum.gradient,
            model=Model(method="UFF", basis=None)
        ),
        initial_molecule=[molecule,],
        optimization_spec=OptimizationSpecification(
            procedure="geomeTRIC",
            keywords={
                "coordsys": "dlc",
                "maxiter": 300,
                "program": "rdkit",
            }
        )
    ),
    procedure="torsiondrive"
)

print(result)

which still give the error

  File "/Volumes/GoogleDrive-108223189278844282544/My Drive/Simulations/sampl8/qforce_error/xtb_TD/run.py", line 59, in <module>
    input_data=TorsionDriveInput(
  File "pydantic/main.py", line 329, in pydantic.main.BaseModel.__init__
  File "pydantic/main.py", line 1022, in pydantic.main.validate_model
  File "pydantic/fields.py", line 865, in pydantic.fields.ModelField.validate
  File "pydantic/fields.py", line 898, in pydantic.fields.ModelField._validate_sequence_like
  File "pydantic/fields.py", line 1064, in pydantic.fields.ModelField._validate_singleton
  File "pydantic/fields.py", line 854, in pydantic.fields.ModelField.validate
  File "pydantic/fields.py", line 1071, in pydantic.fields.ModelField._validate_singleton
  File "pydantic/fields.py", line 1118, in pydantic.fields.ModelField._apply_validators
  File "pydantic/class_validators.py", line 313, in pydantic.class_validators._generic_validator_basic.lambda12
  File "pydantic/main.py", line 684, in pydantic.main.BaseModel.validate
  File "/Users/zhiyiwu/miniconda3_x86/envs/torsiondrive/lib/python3.10/site-packages/qcelemental/models/molecule.py", line 343, in __init__
    from_schema(kwargs, nonphysical=nonphysical), dtype=kwargs["schema_version"], copy=False, np_out=True
  File "/Users/zhiyiwu/miniconda3_x86/envs/torsiondrive/lib/python3.10/site-packages/qcelemental/molparse/from_schema.py", line 46, in from_schema
    frag_pattern = [np.arange(len(ms["symbols"]))]
KeyError: 'symbols'

[jthorton]

@xiki-tempula can you try using this code which works for me I think you might have a geometric molecule as the input instead of a qcelemental molecule.

from qcelemental.models import Molecule, DriverEnum
from qcelemental.models.common_models import Model
from qcelemental.models.procedures import OptimizationSpecification, QCInputSpecification, TDKeywords, TorsionDriveInput
from qcengine import compute_procedure

molecule = Molecule(
    symbols=["C", "C", "H", "H", "H", "H", "H", "H"],
    connectivity=[
        (0, 1, 1),
        (0, 2, 1),
        (0, 3, 1),
        (0, 4, 1),
        (1, 5, 1),
        (1, 6, 1),
        (1, 7, 1),
    ],
    geometry=[
        1.5403406836914142, -1.0173082391323545, 0.9312810207342496,
        4.0719763300123235, -0.0975682592642413, -0.02203578938791481,
        0.0002563605701713713, 0.0013953403968729996, 0.0011121160323317373,
        1.309831306165048, -3.03614919350581, 0.5491856718564878,
        1.3800394103640543, -0.7181256543708329, 2.9707878359388222,
        5.6120991748009645, -1.1161249890160694, 0.907991575289461,
        4.302418801484787, 1.921022388748466, 0.36057345099334853,
        4.232223312568672, -0.3961916040297606, -2.061588178357897]
)

result = compute_procedure(
    input_data=TorsionDriveInput(
        keywords={
            "dihedrals": [(2, 0, 1, 5)],
            "grid_spacing": [60]
        },
        input_specification=QCInputSpecification(
            driver=DriverEnum.gradient,
            model=Model(method="UFF", basis=None)
        ),
        initial_molecule=[molecule,],
        optimization_spec=OptimizationSpecification(
            procedure="geomeTRIC",
            keywords={
                "coordsys": "dlc",
                "maxiter": 300,
                "program": "rdkit",
            }
        )
    ),
    procedure="torsiondrive"
)

print(result)

[xiki-tempula]

@jthorton You are right. Sorry, Torsiondrive use geometric molecule as input so I think this Torsiondrive uses geometric molecule as input as well.