Physical Properties

A core philosophy of this framework is that users should be able to seamlessly curate data sets of physical properties and then estimate that data set using computational methods without significant user intervention and using sensible, well validated workflows.

This page aims to provide an overview of which physical properties are supported by the framework and how they are computed using the different calculation layers.

In this document \(\left<X\right>\) will be used to denote the ensemble average of an observable \(X\).

Density

The density (\(\rho\)) is computed according to

\[\rho = \left<\dfrac{M}{V}\right>\]

where \(M\) and \(V\) are the total molar mass and volume the system respectively.

Direct Simulation

The density is estimated using the default simulation workflow without modification. The estimation of liquid densities is assumed.

MBAR Reweighting

The density is estimated using the default reweighting workflow without modification. The estimation of liquid densities is assumed.

Dielectric Constant

The dielectric constant (\(\varepsilon\)) is computed from the fluctuations in a systems dipole moment (see Equation 7 of [1]) according to:

\[\varepsilon = 1 + \dfrac{\left<\vec{\mu}^2\right> - \left<\vec{\mu}\right>^2} {3\varepsilon_0\left<V\right>k_bT}\]

where \(\vec{\mu}\), \(V\) are the systems dipole moment and volume respectively, \(k_b\) the Boltzmann constant, \(T\) the temperature, and \(\varepsilon_0\) the permittivity of free space.

Note

In v0.2.2 and earlier of the framework the variance was computed as \(\left<{\left(\vec{\mu} - \left<\vec{\mu}\right>\right)}^2\right>\) in order to match the mdtraj implementation which has been used in previous studies by the OpenFF Consortium (see for example [2]). The two approaches should be numerically indistinguishable however.

Direct Simulation

The dielectric is estimated using the default simulation workflow which has been modified to use the specialized AverageDielectricConstant protocol in place of the default AverageObservable protocol. The estimation of liquid dielectric constants is assumed.

MBAR Reweighting

The dielectric is estimated using the default reweighting workflow which has been modified to use the specialized ReweightDielectricConstant protocol in place of the default ReweightObservable protocol. It should be noted that the ReweightDielectricConstant protocol employs bootstrapping to compute the uncertainty in the average dielectric constant, rather than attempting to propagate uncertainties in the average dipole moments and volumes. The estimation of liquid dielectric constants is assumed.

Enthalpy of Vaporization

The enthalpy of vaporization \(\Delta H_{vap}\) (see [3]) can be computed according to

\[\Delta H_{vap} = \left<H_{gas}\right> - \left<H_{liquid}\right> = \left<E_{gas}\right> - \left<E_{liquid}\right> + p\left(\left<V_{gas}\right>-\left<V_{liquid}\right>\right)\]

where \(H\), \(E\), and \(V\) are the enthalpy, total energy and volume respectively.

Under the assumption that \(V_{gas} >> V_{liquid}\) and that the gas is ideal the above expression can be simplified to

\[\Delta H_{vap} = \left<U_{gas}\right> - \left<U_{liquid}\right> + RT\]

where \(U\) is the potential energy, \(T\) the temperature and \(R\) the universal gas constant. This simplified expression is computed by default by this framework.

Direct Simulation

  • Liquid phase: The potential energy of the liquid phase is estimated using the default simulation workflow, and divided by the number of molecules in the simulation box using the divisor input of the AverageObservable protocol.

  • Gas phase: The potential energy of the gas phase is estimated using the default simulation workflow, which has been modified so that

    • the simulation box only contains a single molecule.

    • all periodic boundary conditions have been disabled.

    • all simulations are performed in the NVT ensemble.

    • the production simulation is run for 15000000 steps at a time (rather than 1000000 steps).

    • all simulations are run using the OpenMM reference platform (CPU only) regardless of whether a GPU is available. This is fastest platform to use when simulating a single molecule in vacuum with OpenMM.

The final enthalpy is then computed by subtracting the gas potential energy from the liquid potential energy (SubtractValues) and adding the \(RT\) term (AddValues). Uncertainties are propagated through the subtraction by the normal means using the uncertainties package.

MBAR Reweighting

  • Liquid phase: The potential energy of the liquid phase is estimated using the default reweighting workflow, and divided by the number of molecules in the simulation box using an extra DivideValue protocol.

  • Gas phase: The potential energy of the gas phase is estimated using the default reweighting workflow, which has been modified so that all periodic boundary conditions have been disabled.

The final enthalpy is then computed by subtracting the gas potential energy from the liquid potential energy (SubtractValues) and adding the \(RT\) term (AddValues). Uncertainties are propagated through the subtraction by the normal means using the uncertainties package.

Enthalpy of Mixing

The enthalpy of mixing \(\Delta H_{mix}\left(x_0, \cdots, x_{M-1}\right)\) for a system of \(M\) components is computed according to

\[\Delta H_{mix}\left(x_0, \cdots, x_{M-1}\right) = \dfrac{\left<H_{mix}\right>}{N_{mix}} - \sum_i^M x_i \dfrac{\left<H_i\right>}{N_i}\]

where \(H_{mix}\) is the enthalpy of the full mixture, and \(H_i\), \(x_i\) are the enthalpy and the mole fraction of component \(i\) respectively. \(N_{mix}\) and \(N_i\) are the total number of molecules used in the full mixture simulations and the simulations of each individual component respectively.

When re-weighting cached data to compute \(H_{mix}\) we make the approximation that the kinetic energy contributions cancel out between the mixture and each of the components, and hence can be computed by only re-weighting the NPT reduced potential:

\[\Delta H_{mix}\left(x_0, \cdots, x_{M-1}\right) \approx \dfrac{1}{\beta} \left( \dfrac{\left<u_{mix}\right>}{N_{mix}} - \sum_i^M x_i \dfrac{\left<u_i\right>}{N_i} \right)\]

where \(u \equiv \beta \left(U + pV\right)\) is the NPT reduced potential, \(U\) the potential energy, \(p\) the pressure and \(V\) the volume.

Direct Simulation

  • Mixture: The enthalpy of the full mixture is estimated using the default simulation workflow and divided by the number of molecules in the simulation box using the divisor input of the AverageObservable protocol.

  • Components: The enthalpy of each of the components is estimated using the default simulation workflow, divided by the number of molecules in the simulation box using the divisor input of the AverageObservable protocol, and weighted by their mole fraction in the mixture simulation box using the WeightByMoleFraction protocol.

The final enthalpy is then computed by summing the component enthalpies (AddValues) and subtracting these from the mixture enthalpy (SubtractValues). Uncertainties are propagated through the summation and subtraction by the normal means using the uncertainties package.

MBAR Reweighting

  • Mixture: The reduced potential of the full mixture is estimated using the default reweighting workflow and divided by the number of molecules in the reweighting box using an extra DivideValue protocol.

  • Components: The reduced potential of each of the components is estimated using the default reweighting workflow, divided by the number of molecules in the reweighting box using an extra DivideValue protocol, and weighted by their mole fraction using the WeightByMoleFraction protocol.

The final enthalpy is then computed by summing the component enthalpies (AddValues), subtracting these from the mixture enthalpy (SubtractValues), and multiplying by \(1 / \beta\) (MultiplyValue). Uncertainties are propagated by the normal means using the uncertainties package.

Excess Molar Volume

The excess molar volume \(\Delta V_{excess}\left(x_0, \cdots, x_{M-1}\right)\) for a system of \(M\) components is computed according to

\[\Delta V_{excess}\left(x_0, \cdots, x_{M-1}\right) = N_A \left( \dfrac{\left<V_{mix}\right>}{N_{mix}} - \sum_i^M x_i \dfrac{\left<V_i\right>}{N_i} \right)\]

where \(V_{mix}\) is the volume of the full mixture, and \(V_i\), \(x_i\) are the volume and the mole fraction of component \(i\) respectively. \(N_{mix}\) and \(N_i\) are the total number of molecules used in the full mixture simulations and the simulations of each individual component respectively, and \(N_A\) is the Avogadro constant.

Direct Simulation

  • Mixture: The molar volume of the full mixture is estimated using the default simulation workflow and divided by the molar number of molecules in the simulation box using the divisor input of the AverageObservable protocol.

  • Components: The molar volume of each of the components is estimated using the default simulation workflow, divided by the molar number of molecules in the simulation box using the divisor input of the AverageObservable protocol, and weighted by their mole fraction in the mixture simulation box using the WeightByMoleFraction protocol.

The final excess molar volume is then computed by summing the component molar volumes (AddValues) and subtracting these from the mixture molar volume (SubtractValues). Uncertainties are propagated through the summation and subtraction by the normal means using the uncertainties package.

MBAR Reweighting

  • Mixture: The enthalpy of the full mixture is estimated using the default reweighting workflow and divided by the molar number of molecules in the reweighting box using an extra DivideValue protocol.

  • Components: The enthalpy of each of the components is estimated using the default reweighting workflow, divided by the molar number of molecules in the reweighting box using an extra DivideValue protocol, and weighted by their mole fraction using the WeightByMoleFraction protocol.

The final enthalpy is then computed by summing the component enthalpies (AddValues) and subtracting these from the mixture enthalpy (SubtractValues). Uncertainties are propagated through the summation and subtraction by the normal means using the uncertainties package.

Solvation Free Energies

Solvation free energies are currently computed using the Yank free energy package using direct molecular simulations. By default the calculations attempt to use 2000 solvent molecules, and the alchemical lambda spacings are selected using the built-in ‘trailblazing’ algorithm.

See the Yank documentation for more details.

Host-Guest Binding Free Energy

Warning

The computation of this property is still in beta. Users are heavily recommended to validate any calculations involving this property.

Host-guest binding free energies are currently computed using the attach-pull-release (APR) method [4] through integration with the pAPRika framework.

References

1

Alice Glättli, Xavier Daura, and Wilfred F van Gunsteren. Derivation of an improved simple point charge model for liquid water: spc/a and spc/l. The Journal of chemical physics, 116(22):9811–9828, 2002.

2

Kyle A Beauchamp, Julie M Behr, Ariën S Rustenburg, Christopher I Bayly, Kenneth Kroenlein, and John D Chodera. Toward automated benchmarking of atomistic force fields: neat liquid densities and static dielectric constants from the thermoml data archive. The Journal of Physical Chemistry B, 119(40):12912–12920, 2015.

3

Junmei Wang and Tingjun Hou. Application of molecular dynamics simulations in molecular property prediction. 1. density and heat of vaporization. Journal of chemical theory and computation, 7(7):2151–2165, 2011.

4

David R Slochower, Niel M Henriksen, Lee-Ping Wang, John D Chodera, David L Mobley, and Michael K Gilson. Binding thermodynamics of host–guest systems with smirnoff99frosst 1.0. 5 from the open force field initiative. Journal of Chemical Theory and Computation, 15(11):6225–6242, 2019.