Redundant parameters in Sage 2.1
Several concerns were raised in a private email by Paul Labute (CCG) about Sage 2.1:
t49 vs t84
t49 "*~[#7a]:[#6a:3]~*" in Sage 2.1.- seems to be handled by later t84; should t49 be deleted?
t84 is the below:
"[*:1]~[#7X2,#7X3$(*~[#8X1]):2]:[#6X3:3]~[*:4]"In SMIRNOFF, later parameters (i.e. t84) “override” earlier ones. If we can find an example molecule that t49 is applied to, we should keep it; if not, it’s taking up space.
An example of checking torsions is in the code snippet below:
from openff.toolkit import Molecule, ForceField
sage = ForceField("openff-2.1.0.offxml")
molecule = Molecule.from_smiles("c1cncnc1", allow_undefined_stereo=True)
all_labels = sage.label_molecules(molecule.to_topology())[0]
torsions = all_labels["ProperTorsions"]
for torsion in torsions.values():
print(torsion.id)From looking at the pattern, a charged aromatic nitrogen with a non-Oxygen substituent might be what t49 is applicable to.
A potential solution is iterate through the training datasets from QCArchive to iterate through all molecules to see if there are matches. Getting the coverage of the parameter would be generally interesting in seeing what it applies to vs. t84.
To download from QCArchive:
# Create a client which allows us to connect to the main QCArchive server.
qcarchive_client = FractalClient()
# Retrieve the data set containing the molecules of interest.
from openff.qcsubmit.results import TorsionDriveResultCollection
td_result_collection = TorsionDriveResultCollection.from_server(
client=qcarchive_client,
datasets=[
"OpenFF Gen 2 Torsion Set 1 Roche 2",
"OpenFF Gen 2 Torsion Set 2 Coverage 2",
"OpenFF Gen 2 Torsion Set 3 Pfizer Discrepancy 2",
"OpenFF Gen 2 Torsion Set 4 eMolecules Discrepancy 2",
"OpenFF Gen 2 Torsion Set 5 Bayer 2",
"OpenFF Gen 2 Torsion Set 6 supplemental 2",
],
spec_name="default"
)
# tqdm is for progress bars -- very useful
records_and_molecules = td_result_collection.to_records()
for _, molecule in tqdm.tqdm(records_and_molecules, desc="checking"):
all_labels = sage.label_molecules(molecule.to_topology())[0]However, given that downloading is liable to take a very long time and that Pavan has already put up the records he uses online at , we can just download that instead. (This will be larger than the collection above, as Pavan added additional torsions for training).
Solution
t49 applies to any aromatic CN bond. t84 applies to any CN bond where the N has 2 substituents and the C has 3, or if the N has 3 substituents and one is O. t49 should apply to any molecule which has an aromatic right with a charged N+ in the ring, where the N has 3 substituents but none are O.
We can verify this by inspecting the training dataset to see which molecules need t49 (results are shown in this notebook ). Visualizing the molecules that use t49 shows that it applies to aromatic rings with a charged N, for example:
Therefore, ost aromatic N atoms are treated by t84, but t49 is not redundant, as it applies to aromatic charged N.
t123
t123 "[*:1]~[#15:2]-[#6:3]-[*:4]" in Sage 2.1.0 seems entirely contained in t123a and t124 - should t123 be deleted? The V1 value is suspicious too.
Parameters:
( per is short for "periodicity", ph is short for “phase”. k values have been truncated to the 6th decimal place for conciseness. They’re in kcal/mol. The phase is in degrees. The torsional term has the functional form k*(1+cos(periodicity*theta-phase)))
ID | SMIRKS | per1 | ph1 | k1 | per2 | ph2 | k2 |
|---|---|---|---|---|---|---|---|
t123 |
| 1 | 0 | -10.84539 |
|
|
|
t123a |
| 3 | 0 | 0.112496 |
|
|
|
t124 |
| 2 | 0 | -2.188333 | 3 | 0 | 0.281732 |
A good way to tackle this would be the “coverage” approach above – seeing what kind of training data was used for this parameter might explain the relatively steep force constant.
Solution
t123 corresponds to a C that is single bonded to a P and at least one other group. t123a corresponds to a C that is single bonded to a P and one other substituent, and has 1 other substituents that can have any kind of bond. t124 corresponds to a C that is single bonded to a P and has 2 other substituents that can have any kind of bond.
It is hard to imagine a molecule that has a C with two single bonds that doesn't have either one or two other substituents (e.g. a double bond or two single bonds), though maybe it could happen in some exotic molecule.
We can verify by checking whether any molecules in the training set use this parameter (results in this notebook ):
It appears that no molecules in the training set use t123, which would make it redundant (at least based on this dataset).
It also appears that no molecules in the OpenFF Industry Benchmark Season 1 v1.1 or OpenFF-benchmark-ligand-fragments-v1.0 datasets use t123, as seen in this notebook .
t129
Covered 47 times by Sage 2.1.0/td-set-for-fitting.json and by my filtered version and by my filtered version with torsion multiplicity additions. These 47 torsions cover the 13 molecules here:
And shown with their SMILES here: .
The lingering question is whether or not the large force constant is justified by the torsion potential energy surface.
PESs
To answer this, below are the molecules for which the TorsionDriveRecord dihedrals match the SMIRKS for t129 in the training data set:
Each pair of images shows the molecule on the left, and the corresponding torsion drive scan on the right.
From these data the force constant value of -20 kcal/mol looks reasonable and could even be a bit low for these cases.