|
|
- WARNING The treatment of normalization constants for rigid bonds and rigid angles changed 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 multi-box swap move is attempted, and also chempotperstep times at the end of every Monte Carlo cycle.
μtotal(i) = - kB T ln[ < W * V / ( [N(i)+1] * Λ3(i) ) > ]
where kB 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) = - kB T ln[ < V / N(i) > / Λ3(i)) ]
-
The isolated molecule portion of the chemical potential is the non-density 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 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 configurational-bias.
μisolation(i) = -kB T ln[ < Wisolation > ].
-
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) = - kB 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 isobaric-isothermal 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) = - kB 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 the number of molecules of a given type 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. In the canonical ensemble the residual chemical potential is
μresidual(i) = μNVT Insertion(i) - μIsolation(i)
and in the isobaric-isothermal ensemble (the typical ensemble for computing Henry's law) the residual chemical potential is
μresidual(i) = μNpT Insertion(i) - μIsolation(i).
Return to the Towhee algorithm page
|