Reproducibility

As part of the Torsion multiplicity work, I attempted to reproduce the original Sage 2.1.0 fit and benchmark to make sure my scripts were functioning correctly. As shown in the figure below, my attempt to do so (sage-sage) performed considerably worse than the original Sage 2.1.0 fit. As one check, I cloned the Sage 2.1.0 repo and ran ForceBalance with my conda environment from the torsion multiplicity work, giving the my-sage-2.1.0 line, which again seems worse than Sage 2.1.0. I also had to modify the torsion multiplicity environment with an older version of OpenEye to get ForceBalance to run successfully (see this issue for more details). Other investigation into the torsion multiplicity results did reveal a data set filtering issue in my scripts, but correcting this and rerunning the fit didn’t lead to much improvement (sage-sage-new).

Open

After these observations, I abandoned my scripts and tried to reproduce the original Sage 2.1.0 fit exactly using the input files from the Sage 2.1.0 repo and Pavan’s original conda environment on HPC3 (/dfs4/dmobley-lab/pbehara/conda-env/fb_193). The results are shown in the figure below. My first attempt (Pavan env) performed a bit worse than Sage again, but at Lily’s suggestion I repeated the run with the same inputs and environment to measure the variability across runs. This second attempt (Pavan repeat) performed even better than the original Sage force field, so there is a substantial amount of variation between ForceBalance runs with the same inputs.

Open

I also noticed that some of the ForceBalance runs were not converging fully. The table below shows a collection of force fields, their full and un-penalized final objective function values, and their final convergence statuses.

The poor performance of my attempts to reproduce Sage seems to be caused by a lack of convergence, but this lack of convergence might also be associated with the new environment. The pavan-repeat force field had the best performance, despite requiring the fewest number of optimization iterations. On the other hand, the three force fields optimized in the newer environment took as many, or more, steps as any other force field but failed to converge and had final objective function values nearly double those of the successful runs.

The table below contains a few measures for quantifying the difference between the various Sage 2.1.0 repeats in terms of the DDEs (in kcal/mol) since those were the most obvious differences above.

All of the data reported here is available in my benchmarking repo. In particular, the output/industry subdirectory contains subdirectories for each of the force fields described herein with the correspondences between the names used here and there given in the table below.

Benchmark Repeat

Because these benchmarks are computed with Matt’s new ib benchmarking package, we also wanted to check that the variation was not coming from that. As shown in the figure below, repeating just the benchmark on the pavan-2.1.0 force field produced visually identical results, suggesting that the variation is due to ForceBalance itself, not an artifact of the benchmarking process.

Despite the visually-identical performance here, there are some small differences. Below is a histogram of all of the DDEs differing by more than 1 between the two runs:

These account for 72 records out of 68837 in the full benchmark, and plotting histograms of the results with differences greater than 0.1, for example, is pointless because the bar at 0 completely distorts the rest of the graph. The one record way over near 20 is 36967197 with a difference of 19.79927.

The results are similar for the RMSD, this time plotting the records greater than 0.01:

where the most extreme outlier is again 36967197, and the TFDs also greater than 0.01:

For this, there are actually two records in the bin near 0.4: 36998992 and 36998994.

These three records, 36967197, 36998992, and 36998994 have the following SMILES strings, respectively (in a file because CONFLUENCE wouldn’t stop converting them to links):

Open

Other differences

Similar trends are observed for my run of the Sage 2.1.0 refitting with Pavan’s original input files, but my updated environment with ForceBalance 1.9.5 instead of 1.9.3 used in Pavan’s environment. The figure below shows the DDE differences greater than 8 kcal/mol, which account for 176/68837 records.

On the other hand, there are 8222 entries with differences greater than 1 kcal/mol or about 12%. I think this means that many molecules have small differences rather than a few molecules exhibiting huge differences. The differences are fairly evenly split between being better in the new and old versions, however. The absolute value of the new DDE is lower in 27111 cases and lower in the old data in 35071. The two are equal in the remaining 6655 cases. Restricting the plot above to the cases where the old data is better produces a very similar distribution, but obviously reduces the counts. Notably, this also removes the greatest outlier with a deviation greater than 100 kcal/mol.

The records for the cluster around 80 kcal/mol are 36975868, 36983564, and 36997441. These are clearly much worse in the new force field because the original Sage values are --7.8, --6.0, and --11.6, compared to the new values of --96.6, --94.7, and --101.1, respectively.

I also plotted a CDF for these, but there wasn’t much to gain from it. The old force field increases slightly faster than the new, as expected.

Charge Models

Tom Potter reported an issue when running the Sage 2.1.0 refit with AmberTools/RDKit instead of OpenEye where the initial objective function value was “1.67619e+04, compared to 1.09618e+04 for the original Sage 2.1 run.” Initially he was using versions of the toolkit affected by the charge caching bug, but even after updating to the new versions, he “found it gave essentially the same results.” This was somewhat in line with what I had observed with my own Sage 2.1.0 reproduction runs before the charge caching fix, with some example objective function values from a single run below:

These are somewhat similar to what Tom observed.

Still, we wanted to see what difference the charge model itself would cause (OpenEye AM1-BCC ELF10 vs AmberTools/RDKit AM1-BCC), so I started another Sage 2.1.0 refit starting from the targets and input file from the Sage 2.1.0 repo but with my OpenEye license unexported. I only gave it 6 days of walltime, so it died after iteration 4, but I think that gives us most of the data we wanted to see:

Overall these seem a bit closer to the original Sage 2.1.0 values, both initially and as they converge, than what Tom was observing or what I saw with the charge caching issue.

In case we want to follow up on this in the future, here are my input (and main output) files for the run: and the path on HPC3 is /dfs9/dmobley-lab/bwestbr1/refits/ambertools/new. The targets therein are the same as those used by Tom, which in turn are the same as those used for Sage 2.1.0 minus these targets:

which cause conformer generation errors in RDKit when running ForceBalance. I also found that I had to remove opt-geo-batch-51/18438154-19 for the same reason, and torsion-18537145 also failed 6 times in the initial iteration before finally succeeding, but it didn’t cause any further problems.