Openff-toolkit: Automating QCArchive dataset submission

[jthorton]

Is your feature request related to a problem? Please describe.

In preparation for the bespoke workflow and to also help with general dataset creation we plan to automate the submission of QCArchive datasets to the public and local instances. I think we have decided that this functionality would live in the toolkit but would only be possible to use if users have all the dependencies installed possibly acting through a toolkit wrapper.

I want to restart the conversation around the API from here so I can begin working on this, initially, we would want to start with 1D torsiondrives but the methods should be flexible enough to handle creating any required dataset type.

I like the idea of having a small library (like that in @jchodera example) that can be imported to create custom submission scripts but also including a CLI entry point which should run with all of the default OFF settings. This should allow us to meet the needs of both standard and non-standard submissions, one thing worth pointing out is that currently, we can not specify that we wish to calculate the Hessian following an optimisation at the point of submission see the issue here which will complicate the automation slightly. In the meantime, we may have to use a two-step submission to generate Hessians, where once the user has confirmed the optimisations have completed they run a new CLI command which will create the new data from the optimisation set.
qcfractal-submission --name "My optimisation set" --generator HessianDataset
this would then create a new basic dataset collection with the same name as the optimisation dataset containing the calculations.

[jthorton]

I think we also need to think about how deduplication would work during submissions as when we submit the molecules we may want to put them in canonical order which will be different depending on the backend toolkit used. This means QCFractal would miss duplicate entries and we would waste resources computing the same torsiondrive or optimisation twice.

Recently I experimented with a two-stage searching method to check if a molecule has already been submitted regardless of input ordering using the molecular formula in the first pass and then a strict isomorphism check through the new api Molecule.are_isomorphic(mol1, mol2). This could work provided we cache the archive records to make searching quick.

If we then find a record which matches we then have some logic choices:

• If the record has been seen and is in the current collection we collect and return the results instantly and place no new calculation into the queue.
• If the record has been seen and is in a different collection we use the order in the other record and starting conformation and resubmit the calculation to the new collection, whereby we guarantee QCArchive will determine the procedure to be a duplicate. QCArchive will then _compute_ the procedure where it just sets up links to the old results which will still be quick depending on a worker availability? This may result in a molecule appearing twice in a general force field fitting which uses multiple collections unless we check for duplicates on the collection of the results.
• If the record has been seen and is in a different collection we just return the results along with the name of the collection it came from.

I think the bespoke workflow would always want to return the result instantly to start fitting if possible however, the general fitting collections would probably want all of the results together (option 2) so that it can be tagged with a DOI corresponding to the release. It may be best to let the behaviour at this point be defined by some user option.

[jchodera]

@jthorton : For reference, the toy submission tool I was tinkering with following our discussion in https://github.com/openforcefield/qca-dataset-submission/issues/53 is here:
https://github.com/openforcefield/qcsubmit/blob/master/qcsubmit/factories.py

[jthorton]

Thanks for posting to this, it aligns with my thinking that we would have one base general dataset class (that might be useful to be able to submit as well for benchmark applications as it can hold a lot of single point energies on molecules) that the others would derive from. This could add some strange cases though where users accidentally fragment molecules for an optimization dataset which generally would not be wanted? So would we only add this ability to the torsiondrive datasets?

Also, the programs combined in the workflow bring with them many different running options which we need to control. It may be a good idea to separate out these controls into subcategories corresponding to fragmenter, geometric, torsiondrive and QM to make them easy for users to find and change. For example, the base dataset may have something like this which could be extended to each of the programs mentioned.

@dataclass
class QMOptions:
    theory: str = 'B3LYP-D3BJ'
    basis: str = 'dzvp'
    program: str = 'psi4'
    maxiter: int = 200
    driver: str = 'gradient'
    scf_properties: List = field(default_factory=lambda: ['dipole', 
                                                               'quadrupole',
                                                               'wiberg_lowdin_indices'])
    def request_wbo(self):
        if 'wiberg_lowdin_indices' not in self.scf_properties:
            self.scf_properties.append('wiberg_lowdin_indices')

@dataclass
class GeometricOptions:
    program: str = "geometric"
    coordsys: str = "tric"
    enforce: float = 0.1
    reset: bool = True
    qccnv: bool = True
    epsilon: float = 0.0
    # plus many more 

    def use_tric_coords(sef):
        self.coordsys = 'tric'


class QCFractalDatasetFactory(object):
    """
    Base class with helper functions for creating QCFractal submission datasets.
    Objects from this class should never be created on their own.
    Use their subclasses instead.
    Attributes
    ----------
    fragment : bool, default=False
        If True, will fragment the molecule prior to submission
    enumerate_stereochemistry : bool, default=True
        If True, will enumerate unspecified stereochemistry
    enumerate_tautomers : bool, default=False
        If True, will enumerate tautomers
    max_conformers : int, default=20
        Maximum number of conformers generated per fragment
    input_filters : list, default=[]
        Filters to be applied to input molecules
    submit_filters : list, default=[]
        Filters to be applied to molecules prior to submission
    compute_wbo : bool, default=False
        If True, compute Wiberg-Lowdin bond orders
    """
    def __init__(self):
        """Base class constructor for generating QCFractal datasets.
        """
        self.fragment = False
        self.enumerate_stereochemistry = True
        self.enumerate_tautomers = False
        self.max_conformers = 20
        self.input_filters = list()
        self.submit_filters = list()
        self.compute_wbo = False
        self.qm_options = QMOptions()
        self.geometric_options = GeometricOptions()

With regards to the filters would we want to create an OFF set that could be applied evenly between the toolkit backends as RDKit and OpenEye seem to have different standard filters available.

[j-wags]

Per @jchodera @jthorton and my check-in this morning, This API discussion will be extended until the end of Friday, Feb 21.

@jthorton will make API use cases examples for the input used and the output produced by the following workflows:

• Dump in a bunch of molecules from SDF, have them deduplicated, submitted, retrieve the results
• "Just submit" a dirty SMILES molecule set for QC calculation (some of the SMILES are bad in various ways -- Missing stereochemistry, weird elements, any other common cases we can think of-- What does the workflow do in these cases?)
• Specify computation on a local instance of QCFractal vs. a public server
• Submitting a job which has already been seen before, which immediately retrieves results
• Submission with blocking or non-blocking behavior

[jthorton]

Here are my takes on some example uses of the tools:

• "Just submit" a dirty SMILES molecule set for QC calculation (some of the SMILES are bad in various ways -- Missing stereochemistry, weird elements, any other common cases we can think of-- What does the workflow do in these cases?)

I think the use of the toolkit to import the smiles would cover this for us and handle any errors from the file. The way I see submitter working is to just do the middle work for us so if a user gives molecules that they have from the toolkit using the allow_undefined_stero=True then we will submit any that pass through the chosen filters and make no assumptions about the quality of the molecules.

• Dump in a bunch of molecules from SDF, have them deduplicated, submitted, retrieve the results.

This point I think covers most use cases:

from openforcefield.toolkit import Molecule
from openforcefield.qcsubmit import OptimizationDatasetFactory, Filters

# load the molecules into the toolkit, assuming they are all sensible
molecules = Molecule.from_file('my_collection_of_mols.sdf')

# set up the factory 
opt_generator = OptimizationDatasetFactory()
# now set some of the settings
opt_generator.fragment = False
opt_generator.compute_hessian = False
opt_generator.qm_options.theory = 'wb97x-d'
opt_generator.optimisation_options.maxiter = 350
opt_generator.enumerate_tautomers = True
opt_generator.deduplicate = True
opt_generator.max_conformers = 2
opt_generator.filters = Filters.blockbuster
# add the client information
opt_generator.client = ''public'

# if the user wants to reuse these settings they could then write them to file
opt_generator.export_settings('my_settings.yaml')

# now create the dataset ready for submission
dataset = opt_generator.create_dataset(dataset_name='Optimising my collection', molecules=molecules)

# now we can check details of the generation
# get the total amount of conformers being submitted 
dataset.n_conformers

100

# get the amount of molecules being submitted
dataset.n_molecules

50

# get which molecules failed the filtering
dataset.filtered

Molecule with name 'test 1' and SMILES '[H]C([H])([H])C([H])([H])[H]'

# lets write a receipt for the process
dataset.create_receipt(include_json=True, make_pdf=True, write_settings=True, coverage_report=['openff-1.0.0.offxml', 'openff-1.1.0.offxml''])

# now lets submit the molecules and wait for the results, if the molecules have been seen before 
# the results should be returned very quickly
results = dataset.submit(await_result=True)
# Here the blocking/non-blocking mode can easily be set with this variable on submission. 

# lets look at the optimised energies
results.final_energies

dict{dataset_index: -62.255142.....}

In the case of optimisations, we may have many starting conformers that result in different final geometries and energies, so here we would have a result for each conformer separately according to the dataset index.

The instance of fractal that we submit to could also be assigned two ways I think:

import qcportal as ptl
client = ptl.FractalClient('localhost:7777', 'user', 'password')
opt_generator.client = client

# or we could have from file
opt_generator.load_client('client.yaml')

# where the default would be 
opt_generator.client = 'public'

It might be nice to also have the ability to create a dataset from another users json receipt

from openforcefield.qcsubmit import Dataset

# this will contain all of the molecule inputs and any qm settings ready for submission
dataset = Dataset.from_file('my_collection.json')

A bonus example maybe when we just want to run the same computation again at another level of theory, for example during MM benchmarking with a new forcefield.

generator = OptimizationDatasetFactory()

# lets load in some new settings
generator.load_settings('parsley-2.yaml')

# set the client to public
generator.client = 'public'

# ask the generator to bypass the creation of the dataset 
# just add a new optimisation for each molecules in the set at the new settings. 
results = generator.add_compute(dataset_name='Optimising my collection', await_result=False)

[jthorton]

I guess with this we are also creating a new object that can hold all of the results of a collection as well so it might be nice to be able to make an instance of this for any past dataset with a convenient method rather than having to remake the full dataset and submit it again with the same settings and wait for the results to be collected and returned. I think that would be a case where a slightly miss-typed basis set or theory could result in running a lot of calculations again by mistake, so to avoid this we could have

from openforcefield.qcsubmit import Results
import qcportal as ptl
# load the public client
client = ptl.FractalClient()

# find the dataset you want
ds = client.get_collection('OptimizationDataSet', 'My collection')
# identify the spec name you want 
spec_name = 'wb97x-d'
results = Results.from_collection(collection=ds, spec=spec_name)

# get a json of all of the results together in the whole collection
results.to_json('my_collection_results.json')

This could also help a lot with the benchmarking dashboard as the data will be collected in a consistent way between the datasets and could be loaded into a convenient dictionary, we could also have some quality assurance checks built into the class such as checking the WBO on the optimisation results or warning about internal hydrogen bonds.

[jchodera]

This is pretty great!

My initial thoughts:

I had forgotten that we may need to enumerate protonation states too at some point soon. Extending the API to do this would require adding another opt_generator.expand_protonation_states=True argument, which is totally possible, but leads to the question of whether we want to think about whether we want some way to build composable workflow components rather than just monolithic factories. We can still have monolithic factories that wrap workflows, but there is value in being able to use the components independently via APIs or the CLI.

For example, since much of the pipeline is SMILES-based up until the last step, one could even imagine doing mixing up the CLI with

cat ligands.smi | expand_protonation_states | expand_tautomers | fragment | filter | unique | submit -c torsion_scan -n my_torsion_scan

Having composable workflow components that could be used via their API or CLI (via the entry_points mechanism) could be neat if we can get it essentially for free.

The other thought was that OptimizationDataset and TorsionScan will include many common workflow steps, so we'd want to make sure to avoid having to duplicate code for both.

Otherwise, this is looking great!

[jthorton]

Thanks, glad you like the layout I have in mind.

Ahh yes there will be many other options and operations included and this list could expand a lot over time. I like the sound of your idea my initial thought was that the yaml settings file would essentially encode the workflow for users but this other CLI option would be great. Would this require separate entry points for each of the possible operations to work? instead of say one entry point with multiple arguments like openff-submit -i ligands.smi -settings my_settings.yaml -operations fragment filter expand_tautomers? My only worry with having them all separate is the sheer amount of variable options this would introduce that could be typed for each entry point making something like this might be needed instead of just hiding it all in a settings file unless we made each entry point accept the same settings file?

cat ligands.smi | fragment -wbo_threshold=0.1 -keep_non_rotor_ring_substituents  | filter -blockbuster ... | unique | submit -client client_file.yaml -qm_theroy=wb97x-d -qm_basis=6-311G -c torsion_scan -n my_torsion_scan

One last thing would be where we want this code to live, as you have already started another repo with a good name qcsubmit and we plan to move to namespace style add ons this is probably best suited outside of the toolkit but depending directly on it. So I think the idea would be to continue from what you have already started and move future API discussions to that repo so we can separate them out. If we are all okay with this ill start working on implementing this next week.

[jchodera]

I think we should definitely provide a monolithic workflow driver API and CLI---I was just pointing out the opportunity for also exposing modular workflow components, which would also make it easier to add more components to the workflow later. We can create CLI entry points for individual components later, perhaps even metaprogramming their entry points in setup.py by just discovering all defined subclasses of the common component and exposing their arguments as options since they will support a common API.

A yaml to specify options is very useful, but we want to be careful to not fall into the trap of having to write and extend a large parser that dispatches yaml options to the arguments of different methods using a big translation table. If we kept the workflow components modular, the yaml could be structured like

metadata:
    key: value

workflow_step1_module:
    kwargs1: value
    kwargs2: value

workflow_step2_module:
    kwargs1: value

...

This makes dispatch of the workflow components and options simple and also makes it easy to add new workflow components later since no old code needs updating. Workflow components could every simply be identified at runtime by finding all component subclasses that have been imported.

Some of the things we want to do, such as expanding protonation and tautomeric states, should clearly be implemented in the Molecule class to make our lives easier, and then we can call them in simple composable workflow components.

Not everything needs to be implemented at once. Once we have a modular, extensible design, we can implement the bare minimum needed now and expand as needed later.

I think we are all happy with letting you get started next week!

[j-wags]

The way I see submitter working is to just do the middle work for us so if a user gives molecules that they have from the toolkit using the allow_undefined_stero=True then we will submit any that pass through the chosen filters and make no assumptions about the quality of the molecules.

I think we need to delve into this a bit more. If they set allow_undefined_stereo=True, what stereochemistry do we write for the CMILES?

[j-wags]

My above comment is not blocking for beginning to work on this (sorry, I set today as the deadline for API discussion, and then I waited until the end of the UK workday to post). But this will be a blocking question/potential issue on a final PR.

[jchodera]

I think we need to delve into this a bit more. If they set allow_undefined_stereo=True, what stereochemistry do we write for the CMILES?

We would never generate a CMILES block for a molecule with undefined stereochemistry.

The way this pipeline should work is that _most_ of the methods simply ingest and emit SMILES. At some point, we need to expand the stereochemistry for any SMILES with undefined stereochemistry in order to generate input to one of the QCFractal calculation types.

For molecules with a single stereocenter, simply enumerating one enantiomer should be sufficient. For molecules with two stereocenters, we'll need to enumerate unique combinations.

[j-wags]

We would never generate a CMILES block for a molecule with undefined stereochemistry.

I think this is a good chance to make a OFFMol.is_valid method, where we could query molecules and see if they're missing stereochemistry or other features (like, whether their aromaticity model is actually OEAroModel_MDL). The return type could be an enum of everything wrong with the molecule.

[davidlmobley]

@j-wags note that the OFFMol.is_valid thing you just proposed would also be a good way to resolve some of the "sanity checks" issues I've mentioned another times in previous comments (to a large extent, that users -- academic ones more than pharma -- may feed in molecules that have all kinds of problems without being aware of it).

[jaimergp]

I think this is a good chance to make a OFFMol.is_valid method.

How would you use this in a boolean context? .is_valid() screams for code like:

if molecule.is_valid():
    # proceed
else:
    # report problems

But by returning an Enum the boolean behavior would be inverted. I think it's just a matter of naming... Some ideas:

problems = mol.sanity_checks()
if problems:
    # ooops
else:
    # yay!

---
# using :=

if problems := mol.sanity_checks():
     # oops
else:
    # yay!

Native speakers would probably be able to come up with a better name, but you get the idea.

[jchodera]

I think we want two separate methods here:

• A method that raises an exception of the appropriate type to indicate the kind of problem the molecule has.
• A Boolean property that just returned True (if the molecule is valid) and False (if not)

The boolean property can call the first method and catch the exceptions to return False if any occur.

The first exception-throwing method can be called at the beginning of any operation that requires a valid molecule, and will ensure we are able to throw a consistent set of exceptions that inform the caller of exactly what they did wrong.

[jaimergp]

{Can we / Is it Pythonic to} add extra information to custom exceptions? (aligned with comments in https://stackoverflow.com/questions/6062576/adding-information-to-an-exception).

Something like this:

class CustomError(Exception):
    def __init__(self, details: Dict):
        self.details = details

# we raise something...
raise CustomError({'data': 5})

# and catch it later
except CustomError as e:
    # Do whatever you want with the exception instance
    print(e.details)

[jthorton]

I think it might be best to have different custom enums/exceptions corresponding to each error type so that people can use them in scripts to attempt to automatically clean the molecule like this example using enums:

from openforcefield.descriptors import Warnings 

# hold the problem molecules
problem_molecules = {}

for molecule in molecules:
    # check for problems
    problems = molecule.sanity_checks()

    if problems is not None:
        if Warnings.Stereo in problems:
            # do something to fix it 
            # remove the problem 
            problems.remove(Warnings.Stereo) 

        if Warnings.Valence in problems:
            # collect the molecules for later 
            problem_molecules[molecule] = problems

    else:
        # molecule is fine

[jaimergp]

I like how requests approaches this. If you want to raise an exception they have response.raise_for_errors(). Otherwise, you just have a response.errors list.

That allows for patterns like:

r = requests.get("https://google.com")
r.raise_for_errors()

instead of the more verbose:

try:
    r = requests.get("https://google.com")
except RequestException:
    do something

In our case it'll be:

m = Molecule("some.sdf")
m.raise_for_errors()
# or 
errors = m.sanity_checks(raise_for_errors=True) # default value True or False?
if Warnings.Stereo in errors:
    # fix it

[jaimergp]

Btw, this might be worth a separate issue at this point.