Electrostatic Free Energies for the PB Eq.

Suppose that we have an ensemble of macroions immersed in an electrolytic solution, and we solve the nonlinear PB Eq. to determine the electric potential field $\phi$ anywhere in space. Of how much value is this information? In fact, this information can be processed in a simple way to give the free energy of the macroion/electrolyte system. The free energy as a function of macroion configuration (position and orientation of the various macroions in the system) prescribes the potential of mean force [43] for the macroion assembly. This is a key quantity, since it is an effective interaction potential between the macroions. Knowledge of this interaction potential can then be used to determine the properties of the macroion assembly. For example, a minimum in the potential of mean force corresponds to stable configurations, e.g. a ligand bound to a protein or an enzyme, or a crystalline structure of a charged colloid particle assembly.

How can we process the PB field to obtain a thermodynamic free energy? It turns out that the ``action" S associated with the minimizing $\phi$ configuration is a dimensionless free energy, i.e. $S = - \beta \Phi$ where $\Phi$ is the ``Grand Potential" associated with the Grand Canonical Ensemble (constant $\mu, V, T$) [44]. Such an ensemble is physically appropriate when the system is comprised of a finite set of macroions immersed in a large reservoir of simple ions. In this case the chemical potential is equivalent to the bulk simple ion density, which is a preset input parameter. Sometimes it is more convenient to work in the Canonical Ensemble [43,44,45], in which the volume, temperature and number of ions is specified, e.g., in the case of a liquid or solid of charged colloid particles, with a prescribed number of simple ions (counterions associated with each polyball plus impurity ions that may have been dissolved into the solution). The activity coefficients $\gamma_{\pm}$ (equivalent to chemical potentials) must then be adjusted to put the desired number of particles $n_{\pm}$ in the simulation box. Hence, we ``guess" initial values for $\gamma_{\pm}$, solve the lattice PB Equation (39), then ``refresh" the $\gamma$'s according to:

$$\gamma_{\pm} = n_{\pm} / \sum_n e^{\pm \psi_n - v_n}$$ (41)

The lattice PB Eq. is then solved again until self consistency is achieved (i.e., the differential equation is solved and the output $\psi$ field corresponds to a distribution with the desired total number of simple ions). The Helmholtz free energy can be extracted from the converged $\psi$ field and the appropriate values of $\gamma$ using standard thermodynamic relations, namely [46]:

$$\beta A = n_+ ln \gamma_+ + n_- ln \gamma_- - S$$ (42)