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\)) can be computed from 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.

The framework currently computes \(\varepsilon\) according to

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

making use of the fact that

\[\left<\vec{\mu}^2\right> - \left<\vec{\mu}\right>^2 = \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]).

Direct Simulation

The dielectric is estimated using the default simulation workflow which has been modified to use the specialized ExtractAverageDielectric protocol in place of the default ExtractAverageStatistic. 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 ReweightStatistics. 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 assumtion 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 ExtractAverageStatistic 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.

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 ExtractAverageStatistic 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 ExtractAverageStatistic 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 enthalpy 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 enthalpy 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) and subtracting these from the mixture enthalpy (SubtractValues). Uncertainties are propagated through the summation and subtraction 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 ExtractAverageStatistic 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 ExtractAverageStatistic 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.

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.