
 WARNING The treatment of normalization constants for rigid bonds and rigid angles changed yet again in version 6.2.10. Recomputing the
chemical potential for such molecules is required when comparing between versions.
 WARNING Prior to the 6.0.4 release the treatment of some normalization constants for rigid bond and angles was different. Chemical
potentials for molecules with rigid bonds or angles computed in earlier versions are going to give a different result in newer versions. This is
usually just a constant shift so a quick recalculation of the chemical potential is advised before using old values in a new Grand Canonical
simulation.
 WARNING Prior to the 4.5.2 release all chemical potentials (except the Gibbs total chemical potential) were defined differently than
documented here.
Chemical Potential Definitions
There are several different types of chemical potential output by the Towhee program.
The various chemical potentials reported by the code are listed and defined in this essay.
Consider the following equation as a starting point.
μ_{total}(i) = μ_{residual}(i) + μ_{density}(i)
+ μ_{isolation}(i)
μ_{total}(i) is the total chemical potential of molecule i and is the appropriate value to
consider when looking for total free energies of a system. This total chemical potential is then commonly broken
down into several subcomponents either for computational convenience, or because the subcomponents relate to
certain thermodynamic quantities. Towhee defines, and calculates, these subcomponents as follows.

The total Gibbs chemical potential is computed using the insertion Rosenbluth weight (containing all interactions  intermolecular and intramolecular)
and is computed each time a multibox swap move is attempted, and also chempotperstep times at the end of every Monte Carlo cycle.
μ_{total}(i) =  k_{B} T ln[ < W * V / ( [N(i)+1] * Λ^{3}(i) ) > ]
where k_{B} is Boltzmann's constant, T is the temperature, V is the volume of the simulation box,
N(i) is the number of molecules of type i in the simulation box, and Λ(i) is the
thermal de Broglie wavelength for molecule type i. The < > brackets indicate an ensemble average of the quantity inside the brackets.

The ideal density dependent portion of the chemical potential is
μ_{density}(i) =  k_{B} T ln[ < V / N(i) > / Λ^{3}(i)) ]

The isolated molecule portion of the chemical potential is the nondensity dependent chemical potential of an isolated molecule in an ideal gas
(which means it does not interact with any other molecule, but does interact with itself). This quantity is not normally computed by Towhee.
However, this is computed in the special case where the total number of molecules (nmolectyp) of one of the species is set
to zero. In that case, the isolated molecule chemical potential is computed chempotperstep times per cycle. It uses the full Rosenbluth
weight (W) computed for growing the molecule in an extremely large box that has the Ewald sum and the van der Waals tail corrections disabled.
This uses the standard formula for computing the chemical potential via configurationalbias.
μ_{isolation}(i) = k_{B} T ln[ < W_{isolation} > ].

The NVT insertion chemical potential is computed in a very similar manner as the Gibbs total chemical potential. The only difference is the lack of
any number density terms in the average. This chemical potential calculation is only formally correct for a single box canonical ensemble. However,
experience has shown that this chemical potential is often better when the average number of molecules of the type of interest drops below 1. It turns
out that systematic errors in the number density term from using the Gibbs total chemical potential (essentially the use of N+1) are significant in
this case. This is discussed further in Martin and Siepmann 1998 (TCA).
μ_{NVT Insertion}(i) =  k_{B} T ln[ < W > ]

The NpT insertion chemical potential is computed in a very similar manner as the NVT insertion chemical potential. The only difference is the
inclusion of a volume term in the average, that is then removed later by dividing out the average volume. This chemical potential calculation is
only formally correct for a single box isobaricisothermal ensemble. However, experience has shown that this chemical potential is often better when
the average number of molecules of the type of interest drops below 1. It turns out that systematic errors in the number density term from using the
Gibbs total chemical potential (essentially the use of N+1) are significant in this case. This is discussed further in
Martin and Siepmann 1998 (TCA).
μ_{NpT Insertion}(i) =  k_{B} T ln[ < V W > / < V > ]

The residual chemical potential is not computed directly in Towhee. Instead it is inferred by computing the full insertion chemical potential and
then substracting the isolated molecule chemical potential. This chemical potential is only computed in the special case when we are interested in
Henry's Law coefficients and nmolectyp is set to 0 for at least one of the components. The residual chemical potential is
computed in two slightly different ways depening on the ensemble. If we are using the canonical ensemble then the residual chemical potential is
μ_{residual}(i) = μ_{NVT Insertion}(i)  μ_{Isolation}(i)
and if we are using the isobaricisothermal ensemble (the typical ensemble for computing Henry's law) then the residual chemical potential is
μ_{residual}(i) = μ_{NpT Insertion}(i)  μ_{Isolation}(i).
Return to the Towhee algorithm page
