<jupyter_start><jupyter_text>Submission Preparation We'll start from the JSON files submitted in the `2019-07-25-phenyl-set`:<jupyter_code>!ls *.json<jupyter_output>biphenyls_set_input.json  phenyl_set_torsiondrive_inputs.json<jupyter_text>Our first task will be to take some of the machinery from the original submission script `03_create_torsiondrive_dataset.py` and use it here to create a QCSubmit dataset.**Important**: Must add 1 to every initial_molecule geometry in order to trick the deduplication on the server into letting us submit these molecules.<jupyter_code>import os
import json

import numpy as np

from qcsubmit.factories import TorsiondriveDatasetFactory
from openforcefield.utils.toolkits import RDKitToolkitWrapper, UndefinedStereochemistryError
def read_selected_torsions(input_json):
    """ Read data generated by select_torsions.py

    Returns
    -------
    selected_torsions: dict
        Dictionary for selected torsions, has this structure:
        {
            canonical_torsion_index1: {
                'initial_molecules': [ Molecule1a, Molecule1b, .. ],
                'atom_indices': [ (0,1,2,3) ],
                'attributes': {'canonical_explicit_hydrogen_smiles': .., 'canonical_isomeric_smiles': .., ..}
            },
            ..
        }
    """
    with open(input_json) as infile:
        selected_torsions = json.load(infile)
        
    return selected_torsions<jupyter_output><empty_output><jupyter_text>First, are existing indices unique across both input files?<jupyter_code>inputs = ['biphenyls_set_input.json',  'phenyl_set_torsiondrive_inputs.json']
set(read_selected_torsions(inputs[0]).keys()) & set(read_selected_torsions(inputs[1]).keys())<jupyter_output><empty_output><jupyter_text>Should be good! No overlap!<jupyter_code>def get_qcelmol(torsion_data):
    from qcelemental.models import Molecule
    qcmol = torsion_data['input_molecules']
    
    if isinstance(qcmol, list):
        qcelmols = []
        for qcm in qcmol:
            qcm['extras'] = torsion_data['cmiles_identifiers']
            qcm['geometry'] = (np.array(qcm['geometry']) + 1.0).tolist()
            qcelmols.append(Molecule(**qcm))
    else:
        qcmol['extras'] = torsion_data["cmiles_identifiers"]
        qcmol['geometry'] = (np.array(qcmol['geometry']) + 1.0).tolist()

        qcelmols = [Molecule(**qcmol)]
        
    return qcelmols
factory = TorsiondriveDatasetFactory()
factory.scf_properties
dataset = factory.create_dataset(dataset_name="OpenFF Substituted Phenyl Set 1 v2.0",
                                 tagline="Torsiondrives for selected dihedrals of various phenyl fragments",
                                 description="Torsiondrives for selected dihedrals of various phenyl fragments",
                                 molecules=[])
dataset.metadata.long_description_url = "https://github.com/openforcefield/qca-dataset-submission/tree/master/submissions/2020-10-06-OpenFF-Phenyl-Set"
dataset.metadata.submitter = 'dotsdl'<jupyter_output><empty_output><jupyter_text>Check that molecules don't already exist in QCArchive; that is, that our strategy of shifting the coordinates by 1.0 bohr worked to avoid deduplication.<jupyter_code>from qcportal import FractalClient
client = FractalClient()
def populate_dataset(dataset, input_files, client):
    matches = []
    for input_file in input_files:
        print(f"extracting from '{input_file}'")
        input_mols = read_selected_torsions(input_file)

        for key, moldata in input_mols.items():
            mol = get_qcelmol(moldata)
            
            # we have three known molecules with iodines
            # that we can't converge and want to leave out
            if key in ('CC[CH2:4][NH:3][c:2]1ccnc([cH:1]1)I',
                       'C[C:4](=O)[NH:3][c:2]1[cH:1]cnc(c1)I',
                       '[CH3:4][O:3][c:2]1[cH:1]cnc(c1)I'):
                continue
            
            for m in mol:
                matches.append(client.query_molecules(molecule_hash=m.get_hash()))

            #attributes = factory.create_cmiles_metadata(offmol)
            attributes = mol[0].extras
            print(f"...adding '{key}'")

            dihedral_indices = moldata['dihedral']
            if not isinstance(dihedral_indices[0], list):
                dihedral_indices = [dihedral_indices]

            if len(dihedral_indices) == 0:
                print(f"...NO DIHEDRAL INDICES FOR '{key}'")
                continue

            dihedral_indices = [list(indices) for indices in dihedral_indices]

            dataset.add_molecule(index=key,
                                 molecule=None,
                                 dihedrals=dihedral_indices,
                                 attributes=attributes,
                                 initial_molecules=mol)
    return dataset, matches
ds, matches = populate_dataset(dataset, inputs, client)
sum([len(m) for m in matches])<jupyter_output><empty_output><jupyter_text>We're good to go! No matching molecules in the database.<jupyter_code># upping SCF convergence in case we hit issues on the iodine cases
dataset.maxiter = 800
dataset.scf_properties
dataset.export_dataset('dataset.json.bz2')
dataset.visualize('molecules.pdf')
dataset.molecules_to_file('molecules.smi', 'smi')
confs = np.array([len(mol.conformers) for mol in dataset.molecules])
print("Number of unique molecules       ", dataset.n_molecules)
print("Number of filtered molecules     ", dataset.n_filtered)
print("Number of conformers             ", dataset.n_records)
print("Number of conformers min mean max", 
      confs.min(), "{:6.2f}".format(confs.mean()), confs.max())
dataset.metadata.elements
dataset.qc_specifications<jupyter_output><empty_output><jupyter_text>[scratchspace] Testing convergence of iodine cases Create a dataset with just the problem cases, submit to local QCArchive for test. Instead of filtering here, just use the same functions from above, but only take problem molecules with iodine.<jupyter_code>factory = TorsiondriveDatasetFactory()
factory.scf_properties
dataset = factory.create_dataset(dataset_name="OpenFF Substituted Phenyl Set 1 v2.3",
                                 tagline="Torsiondrives for selected dihedrals of various phenyl fragments",
                                 description="Torsiondrives for selected dihedrals of various phenyl fragments",
                                 molecules=[])
dataset.metadata.long_description_url = "https://github.com/openforcefield/qca-dataset-submission/tree/master/2020-10-06-OpenFF-Phenyl-Set"
dataset.metadata.submitter = 'dotsdl'
def populate_dataset(dataset, input_files):
    for input_file in input_files:
        print(f"extracting from '{input_file}'")
        input_mols = read_selected_torsions(input_file)

        for key, moldata in input_mols.items():
            mol = get_qcelmol(moldata)
            if not 53 in mol[0].atomic_numbers:
                continue
            
            # we have three known failures we want to focus on
            if key not in ('CC[CH2:4][NH:3][c:2]1ccnc([cH:1]1)I',
                           'C[C:4](=O)[NH:3][c:2]1[cH:1]cnc(c1)I',
                           '[CH3:4][O:3][c:2]1[cH:1]cnc(c1)I'):
                continue
                    
            #attributes = factory.create_cmiles_metadata(offmol)
            attributes = mol[0].extras
            print(f"...adding '{key}'")

            dihedral_indices = moldata['dihedral']
            if not isinstance(dihedral_indices[0], list):
                dihedral_indices = [dihedral_indices]

            if len(dihedral_indices) == 0:
                print(f"...NO DIHEDRAL INDICES FOR '{key}'")
                continue

            dihedral_indices = [list(indices) for indices in dihedral_indices]

            dataset.add_molecule(index=key,
                                 molecule=None,
                                 dihedrals=dihedral_indices,
                                 attributes=attributes,
                                 initial_molecules=mol)
    return dataset
ds = populate_dataset(dataset, inputs)<jupyter_output>extracting from 'biphenyls_set_input.json'
extracting from 'phenyl_set_torsiondrive_inputs.json'
...adding 'CC[CH2:4][NH:3][c:2]1ccnc([cH:1]1)I'
...adding 'C[C:4](=O)[NH:3][c:2]1[cH:1]cnc(c1)I'
...adding '[CH3:4][O:3][c:2]1[cH:1]cnc(c1)I'<jupyter_text>This dataset also had iodine cases that won't converge with the current SCF iteration number.We could try upping the SCF iteration number.<jupyter_code># want to test this out on the iodine cases to ensure convergence
dataset.maxiter = 3200<jupyter_output><empty_output><jupyter_text>Try doubling, run on an iodine case, see if we get convergence.<jupyter_code>dataset.scf_properties
from qcportal import FractalClient
client = FractalClient("https://localhost:7777", verify=False)
client
dataset.submit(client, ignore_errors=True)
client.list_collections()
ds = client.get_collection("TorsionDriveDataset", "OpenFF Substituted Phenyl Set 1 v2.3")
ds.status()
ds.df<jupyter_output><empty_output><jupyter_text>1, 4, and 6 failed. We'll focus on re-running these.<jupyter_code>import sys
sys.path.append('../../management/')
from management import get_unfinished_torsiondrive_optimizations, merge, restart_optimizations, restart_torsiondrives, get_unique_optimization_error_messages
opts = get_unfinished_torsiondrive_optimizations(ds, 'default', client)
erred = [opt for opt in merge(opts) if opt.status == 'ERROR']
erred
print(list(get_unique_optimization_error_messages(erred, full=True))[0])
restart_optimizations(erred, client)
erred_tdrs = [tdr for tdr in ds.df.default.values if tdr.status == 'ERROR']
restart_torsiondrives(erred_tdrs, client)
ds = client.get_collection("TorsionDriveDataset", "OpenFF Substituted Phenyl Set 1 v2.3")
ds.status()
opts = get_unfinished_torsiondrive_optimizations(ds, 'default', client)
erred = [opt for opt in merge(opts) if opt.status == 'ERROR']
erred
print(list(get_unique_optimization_error_messages(erred, full=True))[0])<jupyter_output>geomeTRIC run_json error:
Traceback (most recent call last):
  File "/home/david/.conda/envs/qcarchive-worker-openff-psi4/lib/python3.7/site-packages/geometric/run_json.py", line 225, in geometric_run_json
    geometric.optimize.Optimize(coords, M, IC, engine, None, params)
  File "/home/david/.conda/envs/qcarchive-worker-openff-psi4/lib/python3.7/site-packages/geometric/optimize.py", line 1331, in Optimize
    return optimizer.optimizeGeometry()
  File "/home/david/.conda/envs/qcarchive-worker-openff-psi4/lib/python3.7/site-packages/geometric/optimize.py", line 1298, in optimizeGeometry
    self.calcEnergyForce()
  File "/home/david/.conda/envs/qcarchive-worker-openff-psi4/lib/python3.7/site-packages/geometric/optimize.py", line 1002, in calcEnergyForce
    spcalc = self.engine.calc(self.X, self.dirname)
  File "/home/david/.conda/envs/qcarchive-worker-openff-psi4/lib/python3.7/site-packages/geometric/engine.py", line 873, in calc
    return self.calc_new(coords, dirname)
  File "/hom[...]
{'fraction_non_alphanumeric': 0.108002896976281, 'fraction_numerical': 0.015842839036755388, 'mean_word_length': 3.2817829457364343, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '274317', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3621, 'n_tokens_neox': 3376, 'n_words': 726, 'path': 'submissions/2020-10-06-OpenFF-Phenyl-Set/Submission Preparation.ipynb', 'repo_name': 'openforcefield/qca-dataset-submission', 'size': '11046'}