relentless.optimize.RelativeEntropy#

class relentless.optimize.RelativeEntropy(target, simulation, potentials, thermo, extensive=False, T=None)#

Relative entropy.

The relative entropy \(S_{\rm rel}\) (or Kullback-Leibler divergence) quantifies the overlap of two probability distributions. For a known target statistical mechanical ensemble having distribution \(p_0\) and a simulated model ensemble having distribution \(p\) and parametrized on a set of design variables \(\mathbf{x}\), the relative entropy from the model to the target is:

\[S_{\rm rel} = -\int d\Gamma p_0(\Gamma)\ln\left(\frac{p(\Gamma)}{p_0(\Gamma)}\right)\]

where \(\Gamma\) is an element of phase space. The relative entropy is zero when the two ensembles overlap completely, and it is positive otherwise. Hence, minimization of \(S_{\rm rel}\) can be used to find parameters of a model that reproduce a known target ensemble.

The value of the relative entropy is not readily determined in molecular simulations, so this ObjectiveFunction does not return a value. However, the gradient of the relative entropy with respect to the design variables \(\mathbf{x}\) is much easier to compute as ensemble averages.

The gradient of the relative entropy can be computed through direct ensemble averaging or through integration of the radial distribution function (rdf) of the target and model ensembles.

Direct ensemble average

When target is a str and thermo is a WriteTrajectory object, the relative entropy is computed through direct ensemble averaging. In this case, the relative entropy gradient is computed as:

\[\nabla_\mathbf{x} S_{\rm rel} = \beta \left( \left\langle \frac{\partial U}{\partial \mathbf{x}} \right\rangle_0 - \left\langle \frac{\partial U}{\partial \mathbf{x}} \right\rangle \right)\]

where \(\langle\cdot\rangle_0\) and \(\langle\cdot\rangle\) are ensemble averages in the target and model ensembles, respectively, \(\beta=1/(k_{\rm B}T)\), and \(U\) is the total potential energy.

It is recommended to use the direct ensemble average method in most cases, as it allows for design of PairPotential and bonded interactions. This method also easily supports design of pair potentials with bonded exclusions.

Radial distribution function

The RelativeEntropy objective function supports a method for designing PairPotential only interactions. These interactions are characterized by \(g_{ij}(r)\), an RDF for each pair of interacting types \((i,j)\) in each Ensemble. The gradient of \(S_{\rm rel}\) is then:

\[\nabla_\mathbf{x} S_{\rm rel} = -\frac{1}{2} \sum_{i,j} \int dr 4\pi r^2 \left[\frac{\beta N_i N_j}{V} g_{ij}(r) \right. \left. -\frac{\beta_0 N_{i,0} N_{j,0}}{V_0} g_{ij,0}(r)\right] \nabla_\mathbf{x} u_{ij}(r)\]

where \(\beta=1/(k_{\rm B}T)\), \(N_i\) is the number of particles of type \(i\), \(V\) is the extent, and \(u_{ij}(r)\) is the pair potential in the model ensemble. The corresponding properties of the target ensemble are denoted with subscript \(0\).

This method is performed when Ensemble and EnsembleAverage objects are provided for the target and thermo parameters, respectively. While it is possible to use this method with bonded exclusions by setting the target ensemble’s rdf and the model ensemble’s rdf calculation to have the same exclusions and sampling, in practice this proves to be difficult.

\(S_{\rm rel}\) is extensive as written, meaning that it depends on the size of the system. This can be undesirable for optimization because it means certain hyperparameters are system-dependent, so the default behavior is to normalize \(s_{\rm rel}=S_{\rm rel}/V_0\). To use the extensive relative entropy set extensive=True.

Parameters:

  • target (Ensemble or str) – The target Ensemble (must have specified V and N). The str should be the path to a gsd trajectory file containing the target ensemble.

  • simulation (Simulation) – The simulation engine to use, with specified simulation operations.

  • potentials (Potentials) – The pair potentials to use in the simulations.

  • thermo (EnsembleAverage or) – WriteTrajectory The thermodynamic analyzer operation for the simulation ensemble. The model ensemble will be extracted from this operation.

  • extensive (bool) – Specification of whether the relative entropy is extensive (defaults to False).

  • T (float) – Temperature of the simulation ensemble. If not specified with the direct ensemble average method, it will be computed from the simulation trajectory.

Methods

compute(variables[, directory])

Evaluate the value and gradient of the relative entropy function.

compute_gradient(ensemble, variables)

Computes the relative entropy gradient for an ensemble.

Attributes

target

The target ensemble.
compute(variables, directory=None)#

Evaluate the value and gradient of the relative entropy function.

The value of the relative entropy is not computed, but the gradient is. Calculating the gradient requires running a simulation, which may be computationally expensive.

Optionally, a directory can be specified to write the simulation output as defined in run(), namely the values of the pair potential design variables, which are written to pair_potential.i.json for the \(i\)th pair potential, and the ObjectiveFunctionResult, which is written to result.json.

Parameters:

  • variables (Variable or tuple) – Variables with respect to which to compute gradient.

  • directory (str or Directory) – The output directory

Returns:

The result, which has unknown value None and known gradient.

Return type:

ObjectiveFunctionResult

compute_gradient(ensemble, variables)#

Computes the relative entropy gradient for an ensemble.

Parameters:

  • ensemble (Ensemble) – The ensemble for which to evaluate the gradient.

  • variables (Variable or tuple) – Variables with respect to which to compute gradient.

Returns:

The gradient, keyed on the variables.

Return type:

KeyedArray

Raises:

ValueError – If the variables have non-zero derivatives with respect to bonded potentials.

property target#

The target ensemble. Must have both V and N parameters set.

Type:

Ensemble