Trouble Reproducing Results - Bond Energy Function definition

[manglav]

I'm running through the code and am a bit confused on some of the code choices. For example, For example, B1 and B2 solvent files have different lambda restraints.

bond_energy_function = "lambda_restraints*(K/2)*(r-r0)^2;"
bond_energy_function = "lambda_restraints"

I would expect these to be the same. I'm assuming I'm missing something - any pointers to what it is?

[djhuggins]

In the context of the solvent calculations, there are no harmonic bonds because there is no restraint. These are only used in the complex and restraint simulations. The setting of bond_energy_function should be the same across all ligands in the complex simulations and the restraint simulations. The bond_energy_function was set here due to the automated way I wrote the scripts.

[manglav]

Ah makes sense. I cleaned it up for the dead code parts and I somehow fixed my error, so thank you!

I have one more quick question - I'm trying to reproduce your paper on OpenMM8, but the correction term seems hardcoded in reach results.py. Is the source of those calculations available? I've seen a few equations floating around either relating the bounding box volume or concentrations ratios, but thought I'd ask. Thank you very much for this! This is the most comprehensive example of alchemical free energies using pure OpenMM/Tools on the entire internet.

[djhuggins]

I am using what are commonly called Boresch restraints or VBA. Eqn 14 in this paper:

https://pubs.acs.org/doi/10.1021/jp0217839

V is set as 1661 cubic Angstroms. To make the units kcal/mol, the prefactor is set at 0.593 kcal/mol.