oz_doc.tex

% This file is part of SunlightDPD - a home for open source software
% related to the dissipative particle dynamics (DPD) simulation
% method.

% Original copyright (c) 2009-2018 Unilever UK Central Resources Ltd
% (Registered in England & Wales, Company No 29140; Registered Office:
% Unilever House, Blackfriars, London, EC4P 4BQ, UK).  Later
% modifications copyright (c) 2020-2025 Patrick B Warren
% <patrick.warren@stfc.ac.uk> and STFC.

% SunlightDPD is free software: you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation, either version 3 of the
% License, or (at your option) any later version.

% SunlightDPD is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.

% You should have received a copy of the GNU General Public License
% along with SunlightDPD.  If not, see <http://www.gnu.org/licenses/>.

\documentclass[12pt,a4paper]{article}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
%\usepackage{color}
\usepackage{upquote}

\DeclareMathOperator{\erf}{erf}
\DeclareMathOperator{\erfc}{erfc}
\DeclareMathOperator{\imaginary}{Im}
\DeclareMathOperator{\Tr}{Tr}

\newcommand{\latin}[1]{\emph{#1}}
\newcommand{\german}[1]{$\frak{#1}$}
\newcommand{\etal}{\latin{et al.}}
\newcommand{\etc}{\latin{et\,c.}}
\newcommand{\eg}{\latin{e.\,g.}}
\newcommand{\cf}{\latin{c.\,f.}}
\newcommand{\ie}{\latin{i.\,e.}}
\newcommand{\via}{\latin{via}}
\newcommand{\aposteriori}{\latin{a posteriori}}
\newcommand{\adhoc}{\latin{ad hoc}}
\newcommand{\viz}{\latin{viz.}}
\newcommand{\viceversa}{\latin{vice versa}}
\newcommand{\mutmut}{\latin{mutatis mutandis}}
\newcommand{\ansatz}{{ansatz}}
\newcommand{\perse}{\latin{per se}}

\newcommand{\french}[1]{\emph{#1}}
\newcommand{\ala}{\french{\`a la}}

\newcommand{\half}{{\textstyle\frac{1}{2}}}

\newcommand{\alt}{\lesssim}
\newcommand{\agt}{\gtrsim}

\newcommand{\kB}{k_{\mathrm{B}}}
\newcommand{\kT}{\kB T}
\newcommand{\myex}{^{\mathrm{ex}}}
\newcommand{\pex}{p\myex}
\newcommand{\Uex}{U\myex}
\newcommand{\uex}{u\myex}
\newcommand{\Aex}{A\myex}
\newcommand{\aex}{a\myex}
\newcommand{\muex}{\mu\myex}
\newcommand{\myvec}[1]{{\mathbf #1}}
\newcommand{\rvec}{\myvec{r}}
\newcommand{\kvec}{\myvec{k}}
\newcommand{\lr}{^{\mathrm{L}}}
\newcommand{\vlr}{v\lr}
\newcommand{\tildevlr}{{\tilde v}\lr}
\newcommand{\myprime}{^{{}\,\prime}}
%\newcommand{\sigmap}{{\overline\sigma}}
%\newcommand{\sigmap}{{\Sigma}}
\newcommand{\sigmap}{{\sigma'}{}}
\newcommand{\gc}{g^{c}}
\newcommand{\href}{h^{(0)}}

\newcommand{\varbeta}{\text{\ss}}

% The following are used to typeset a box of width character '0'
\newlength{\zerolen}
\settowidth{\zerolen}{0}
\newcommand{\zeroset}[1]{\makebox[\zerolen][c]{#1}}

\newcommand{\Eqref}[1]{Eq.~\eqref{#1}}
\newcommand{\Eqsref}[1]{Eqs.~\eqref{#1}}

\newcommand{\Refcite}[1]{Ref.~\cite{#1}}

\newcommand{\myav}[1]{\langle #1\rangle}

\newcommand{\FLAG}[1]{{\color{red} #1}}

\newcommand{\FORTRAN}{{\sc fortran}}
\newcommand{\python}{{\tt python}}
\newcommand{\LAPACK}{{\sc lapack}}
\newcommand{\FFTW}{{\sc fftw}}

\begin{document}

\noindent{\bf\large FORTRAN 90 Ornstein-Zernike solver :: version 1.13}\\[6pt]
\noindent{\it Patrick B Warren, Hartree Centre, STFC (July 2021)}\\[6pt]

\noindent The code and this accompanying documentation is released
under the GPL.  You are welcome to use and modify it for your own
projects. If you end up publishing results based on this, please cite\\
%
\begin{center}
  \fbox{\parbox{0.9\textwidth}{P.~B.~Warren, A.~Vlasov, L.~Anton and
    A.~J.~Masters ``Screening properties of Gaussian electrolyte
    models, with application to dissipative particle dynamics'',
    J. Chem. Phys. {\bf138}, 204907 (2013).}}
\end{center}

\vspace{0.25in}

\noindent Features of version 1.13:
\begin{itemize}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
\item modern \FORTRAN\ 90 based, with example \python\ driver scripts;
\item HNC, MSA, RPA, and EXP closures;
\item fast Ng solver, and Ng decomposition for electrostatics;
\item multicomponent (arbitrary number of components);
\item hard core (RPM-like) and soft core (DPD-like) potentials;
\item thermodynamics computed (full HNC free energy added in v1.12);
\item fully open source.
\end{itemize}

\noindent
{\small Versions 1.13 and 1.12 were feature enhancements. Versions
  1.10, 1.11 fixed bugs: the contact virial pressure contribution was
  missing from the HNC free energy, and the compressibility and HNC
  chemical potentials were missing a mean-field term for the off-SYM
  RPM potential.  Added in version 1.7--1.9: MSA closure and EXP;
  exact and approximate pair potentials for Slater smearing; EXP for
  hard core potentials fixed, new examples added; source code brought
  up to a modern \FORTRAN\ 90 standard.}

\newpage
{\small\tableofcontents}
\newpage

\section{Introduction}
%
The \FORTRAN\ 90 module \verb+oz_mod.f90+ implements an
Ornstein-Zernike solver for multicomponent mixtures of possibly
charged particles, returning both structural and thermodynamic
information \cite{WVA+13}.  The closures currently implemented are the HNC, MSA,
RPA, and EXP.  Though we say it ourselves, the code is \emph{fast}
since it is written in native \FORTRAN\ and implements the Ng
acceleration schemes \cite{Ng74}.  The interaction potentials are
specified in their own routines and the code could be used for other
potentials with minor modifications.  Further technical background to
liquid state integral equation methods can be found in Hansen and
McDonald \cite{HM06}, Kelley and Montgomery Pettitt \cite{KMP04}, and
Vrbka \etal\ \cite{Vrbka09}.  We note that Vrbka provides an
independent package called `pyOZ', implemented entirely in
\python.  There are some differences with the present
approach: the Ornstein-Zernike relation is normalised slightly
differently, a conjugate gradient method is used to accelerate
convergence, and pyOZ provides other closures in addition to the ones
provided here.  The present code is based on an HNC code developed by
Lucian Anton, modified for dissipative particle dynamics (DPD) by
Andrey Vlasov with help from Lucian and Andrew Masters.  The code was
later converted to use the open source libraries \FFTW\ and
\LAPACK\ by Patrick Warren, and rewritten to be compatible with
\python\ using the \verb+f2py+ interface generator.

\section{Theoretical background}
%
\subsection{Single component HNC}
%
We suppose that the interaction between a pair of particles is given
by the potential $v(r)$.  We introduce the following quantities:
%
\begin{equation}
  \begin{array}{lll}
    g(r)        &\quad& \text{pair distribution function},   \\
    h(r)=g(r)-1      && \text{total correlation function},   \\[6pt]
    c(r)             && \text{direct correlation function},  \\
    e(r) = h(r)-c(r) && \text{indirect correlation function},
  \end{array}
\end{equation}
%
where $c(r)$ is defined by the Ornstein-Zernike relation below.  Note
that the indirect correlation function is called $b(r)$ by Hansen and
McDonald \cite{HM06} and $\gamma(r)$ by Vrbka \etal\ \cite{Vrbka09}.
We also introduce the Fourier transforms, \eg\
%
\begin{equation}
\tilde h(k) = \int\!d^3\rvec\, e^{-i\kvec\cdot\rvec} \,h(r)\,,
\end{equation}
%
which simplifies to the Fourier-Bessel transform
%
\begin{equation}
\tilde h(k) = \frac{4\pi}{k} \int_0^\infty \!\! dr\, \sin(kr)\, r\, h(r)\,,
\label{eq:fFB}
\end{equation}
%
and
%
\begin{equation}
h(r) = \int\! \frac{d^3\kvec}{(2\pi)^3} \,e^{i\kvec\cdot\rvec} \,\tilde h(k)\,,\label{eq:ftinv}
\end{equation}
%
which simplifies to 
%
\begin{equation}
h(r) = \frac{1}{2 \pi^2 r} \int_0^\infty \!\!dk\, \sin(kr)\, k\, \tilde h(k)\,.
\label{eq:bFB}
\end{equation}
%
In these terms the Ornstein-Zernike (OZ) relation is~\cite{HM06}
%
\begin{equation}
\tilde h(k) = \tilde c(k) + \rho\,
\tilde h(k)\, \tilde c(k)
\label{eq:oz1a}
\end{equation}
%
where $\rho$ is the density.  This confirms the relevance of the
standard choice of normalisation of the 3d Fourier transform pair
which puts the factor $1/(2\pi)^3$ into the back transform.  The OZ
relation can be rearranged to 
%
\begin{equation}
\tilde h(k) = \frac{\tilde c(k)}{1-\rho \tilde c(k)}
\label{eq:oz1b}
\end{equation}
%
(\cf\ Eq.~(3.5.13) in \Refcite{HM06}).
The hyper-netted chain (HNC) closure is defined in real space and is
%
\begin{equation}
g(r)=\exp[-\beta v(r)+h(r)-c(r)]
\label{eq:hnc1a}
\end{equation}
%
(\cf\ Eq.~(4.3.19) in \Refcite{HM06}) where $\beta=1/\kT$.  It
amounts the neglect of the bridge diagrams.  For hard spheres HNC is
known to be inferior to Percus-Yevick, but for soft potentials like
DPD it generally gives excellent results.  In the presence of hard
cores, we note that $\exp[-\beta v]$ is still well defined (and zero,
in fact) where $\beta v$ is infinite.  We can accommodate hard cores
therefore by making sure that we always work with $\exp[-\beta v]$
instead of $\beta v$ in the numerics.  The function $\exp[-\beta v]$
can be pre-computed for a given potential.

To obtain the algorithm implemented in the code we rewrite
\Eqsref{eq:oz1b} and \eqref{eq:hnc1a} in terms of $c(r)$ and
$e(r)$.  The OZ relation becomes
%
\begin{equation}
\tilde e = \frac{\tilde c}{1-\rho\tilde c}-\tilde c
\label{eq:oz1c}
\end{equation}
%
and the HNC closure becomes
%
\begin{equation}
c=\exp[{-\beta v}+e]-e-1\,.
\label{eq:hnc1b}
\end{equation}
%
The algorithm is as follows.  We start with some guess for the direct
correlation function $c(r)$.  We Fourier transform this to $\tilde
c(k)$ and solve the OZ relation in \Eqref{eq:oz1c} to get the indirect
correlation function (in reciprocal space) $\tilde e(k)$.  We Fourier
back-transform to get $e(r)$ and plug this into the HNC closure in the
form of \Eqref{eq:hnc1b} to get a new direct correlation function
$c(r)$.  This cycle is iterated until $c(r)$ and $e(r)$ converge to a
self-consistent solution pair.  A direct approach like this is usually
numerically unstable so in the Picard scheme we mix a little of the
new $c(r)$ into the old $c(r)$, and iterate until we converge to a
solution.  In the Ng acceleration scheme we keep track of the
convergence trajectory and use this to accelerate convergence to the
solution.  The details of the Ng scheme will not be described here,
rather we point the interested reader to \Refcite{Ng74} and the code
itself.  The initial trajectory for the Ng scheme is generated by a
few Picard iterations.  Some possible initial choices for the direct
correlation function $c(r)$ are zero, $-\beta v(r)$ and $\exp[-\beta
  v(r)]-1$.  For soft potentials any of these will do in principle but
the initial rate of convergence may differ.  For hard core potentials,
the middle option cannot be used.
 
Now we describe Ng's splitting method for handling long range forces.
We partition the interaction potential into a short range part and a
long range part,
%
\begin{equation}
v(r)=v'(r)+v\lr(r)\,.
\end{equation}
%
It is generally accepted that $c(r)\sim-\beta v\lr(r)$ for
$r\to\infty$ so we make this explicit by partitioning both the direct
and indirect correlation functions,
%
\begin{equation}
c(r)=c'(r)-\beta v\lr(r)\,,\quad
e(r)=e'(r)+\beta v\lr(r)\,.
\end{equation}
%
These are taken to provide the definitions of the offset functions
$c'(r)$ and $e'(r)$.  We rewrite the OZ relation and the HNC closure
in terms of these,
%
\begin{equation}
\tilde e\myprime=\frac{\tilde c\myprime-\beta\tildevlr}{1-\rho(\tilde
  c\myprime-\beta\tildevlr)}-\tilde c\myprime
\label{eq:oz1d}
\end{equation}
and
\begin{equation}
c'=\exp[-\beta v'+e']-e'-1\,.
\label{eq:hnc1d}
\end{equation}
%
The first of these requires that we know the Fourier transform of the
long range part of the potential.  The second requires that we know
the exponentiated short range potential.  Any hard core interaction
can be parked safely in the latter but the consequential freedom to
specify $v\lr(r)$ inside the hard core should be used wisely!

The main advantage of the Ng split is that $c'(r)$ and $e'(r)$ are now
genuinely short-ranged so the numerical behaviour is much better,
certain thermodynamic integrals are guaranteed convergence, and the
expected asymptotic behaviour of $c(r)$ is explicitly incorporated.
The numerical solution scheme based on \Eqsref{eq:oz1d} and
\eqref{eq:hnc1d} goes through as before.  The initial choices for
$c'(r)$ are replicated from above with $v$ replaced by $v'$.

\subsection{Multicomponent HNC}
%
Now we turn to the multicomponent problem.  The pair functions
become for example $g_{\mu\nu}(r)$, \etc, and the problem is specified
by the interaction potential $v_{\mu\nu}(r)$ and the species densities
$\rho_\mu$.  The total density is $\rho=\sum_\mu \rho_\mu$ and the
species mole fractions are $x_\mu=\rho_\mu/\rho$.  We again split the
potential into short and long range contributions, 
%
\begin{equation}
v_{\mu\nu}(r)=v'_{\mu\nu}(r)+ v\lr_{\mu\nu}(r)
\end{equation}
%
where typically $v\lr_{\mu\nu}(r)$ incorporates the long range part of
the charge interactions.  Often the long range part satisfies a
symmetry condition (SYM) where
%
\begin{equation}
v\lr_{\mu\nu}(r)=z_\mu z_\nu v\lr(r)\quad\text{(SYM)}
\label{eq:sp}
\end{equation}
in which $z_\mu$ is the valency and $v\lr(r)$ is the interaction
potential between unit charges of the same sign.  However we will not
assume that $v\lr_{\mu\nu}(r)$ necessarily meets the SYM condition.  As
above we define $c'_{\mu\nu} = c_{\mu\nu} + \beta v\lr_{\mu\nu}$ and
$e'_{\mu\nu} = e_{\mu\nu} - \beta v\lr_{\mu\nu}$.  The multicomponent
HNC closure becomes
%
\begin{equation}
c_{\mu\nu}'=\exp[-\beta v_{\mu\nu}'+e_{\mu\nu}']-e_{\mu\nu}'-1\,.
\label{eq:hncnd}
\end{equation}
%
To derive the OZ relation in a usable form we start from 
the multicomponent analogue of \Eqref{eq:oz1a},
%
\begin{equation}
{\tilde h}_{\mu\nu} = {\tilde c}_{\mu\nu}+\rho\,{\textstyle\sum_\lambda}
\,x_\lambda\,{\tilde c}_{\mu\lambda}\,{\tilde h}_{\lambda\nu}\,.
\end{equation}
%
We write this as a matrix relation, introducing $R$ as a diagonal
matrix with entries $\rho_\mu=\rho x_\mu$,
%
\begin{equation}
{\tilde H}={\tilde C}+{\tilde C}\cdot R\cdot{\tilde H}\,.
\label{eq:oz2}
\end{equation}
%
Note that all the matrices in this are symmetric, so we can equally
well write this as
%
\begin{equation}
{\tilde H}={\tilde C}+{\tilde H}\cdot R\cdot{\tilde C}\,.
\label{eq:oz2a}
\end{equation}
%
Introducing next ${\tilde C}' = {\tilde C} + \beta{\tilde U}\lr$ and
${\tilde E}'={\tilde E}-\beta{\tilde U}\lr$ in \Eqref{eq:oz2}, and
rearranging, gives eventually
%
\begin{equation}
{\tilde E}'=[I-({\tilde C}'-\beta{\tilde U}\lr)\cdot R]^{-1}\cdot
[({\tilde C}'-\beta{\tilde U}\lr)\cdot R\cdot {\tilde C}'-\beta{\tilde
    U}\lr]\,.
\label{eq:oznd}
\end{equation}
%
The solution scheme is essentially as before.  Given an initial guess
for $c_{\mu\nu}'(r)$ we Fourier transform, evaluate the matrix
expression in \Eqref{eq:oznd}, and Fourier back-transform to obtain
$e_{\mu\nu}'(r)$.  We substitute this into the HNC closure in
\Eqref{eq:hncnd} to obtain a new guess for $c_{\mu\nu}'(r)$.  The
cycle is iterated to convergence and in practice a few Picard
iterations are used to pump-prime a multicomponent version of the Ng
acceleration scheme.  As before, any hard core part of the potential
should be parked in $\exp[-\beta v_{\mu\nu}']$.

\subsection{Multicomponent MSA}
%
The closure in the mean spherical approximation (MSA) is
$g_{\mu\nu}(r)=0$ for $r\le d_{\mu\nu}$, and $c_{\mu\nu}(r)= -\beta
v_{\mu\nu}(r)$ for $r> d_{\mu\nu}$.  Here, $d_{\mu\nu}$ is the range
of the hard core repulsion between species $\mu$ and $\nu$
(\ie\ $v'_{\mu\nu}=\infty$ for $r\le d_{\mu\nu}$).  Rewriting as usual
in terms of our offset functions the MSA closure becomes
%
\begin{equation}
  c'_{\mu\nu}=\Bigl\{\begin{array}{ll}
  -e'_{\mu\nu}-1 & (r<d_{\mu\nu})\,,\\[3pt]
  -\beta v'_{\mu\nu} & (r>d_{\mu\nu})\,.
  \end{array}
\end{equation}
%
This is almost a drop-in replacement for the HNC closure in
\Eqref{eq:hncnd}.  For a cold start we can initialise $c'_{\mu\nu}=-1$
for $r\le d_{\mu\nu}$, and we leave $c'_{\mu\nu}=-\beta v'_{\mu\nu}$
untouched for $r>d_{\mu\nu}$ during the iterative solution.  We make
sure that whenever the potential is changed $c'_{\mu\nu}$ for
$r>d_{\mu\nu}$ is correctly set.

\subsection{Multicomponent RPA}\label{sec:RPA}
%
The RPA closure is a special case of the MSA when there are no hard
cores, and is represented by $c_{\mu\nu}(r)=-\beta v_{\mu\nu}(r)$.
This is in fact one of the choices for initialising the HNC solver.
Given the HNC machinery, the implementation of the RPA is trivial and
comprises a single round trip through the OZ solver starting from
$c'_{\mu\nu}(r)=-\beta v'_{\mu\nu}(r)$.  No iteration is required.

\subsection{Multicomponent EXP}\label{sec:EXP}
%
The EXP approximation is a  refinement of the MSA/RPA
in which the  total correlation functions are replaced by
%
\begin{equation}
  h_{\mu\nu}\to(1+\href_{\mu\nu})\exp(h_{\mu\nu}-\href_{\mu\nu})-1\,,
  \label{eq:exp}
\end{equation}
%
in which $\href_{\mu\nu}$ will be a hard sphere reference fluid in
the case of MSA, and $\href_{\mu\nu}=0$ in the case of RPA (no hard
cores).  The EXP breaks unwarranted symmetries such as
$h_{++}=-h_{+-}$ and ensures that $h_{\mu\nu}(r)>-1$ which is not
always satisfied by the MSA/RPA.  In practice EXP can be nearly as
good as HNC and extends to state points where HNC doesn't have a
solution.

Since the EXP approximation is defined as an action on the correlation
functions in real space, a round trip through the OZ relation is
required to compute the corresponding offset direct and indirect
correlation functions.  To do this we rewrite the OZ relation in
\Eqref{eq:oz2a} as
%
\begin{equation}
{\tilde C} = (I+{\tilde H}\cdot R)^{-1}\cdot{\tilde H}\,.
\label{eq:ozndb}
\end{equation}
%
With this implemented, $c'_{\mu\nu}=c_{\mu\nu}+\beta
v\lr_{\mu\nu}$ and $e'_{\mu\nu}= h_{\mu\nu} - c'_{\mu\nu}$ follow.

\subsection{Structural properties}
%
\subsubsection{Pair functions}
%
Given a converged solution pair $(c_{\mu\nu}', e_{\mu\nu}')$, the
total correlation functions can be evaluated from
%
\begin{equation}
h_{\mu\nu}(r) = c'_{\mu\nu}(r) +  e'_{\mu\nu}(r)\label{eq:hrs}
\end{equation}
%
(note that $\beta v_{\mu\nu}\lr$ cancels).  The pair functions are
given by $g_{\mu\nu}=\delta_{\mu\nu}+h_{\mu\nu}$.

\subsubsection{Structure factors}
%
For the partial structure factors, there is a choice of normalisation.
The present code defines
%
\begin{equation}
    S_{\mu\nu}(k)
    = \rho_\mu\delta_{\mu\nu} + \rho_\mu\rho_\nu{\tilde h}_{\mu\nu}
    \qquad \text{(present definition)}\,,
  \label{eq:skdef}
\end{equation}
%
where ${\tilde h}_{\mu\nu}={\tilde c}_{\mu\nu}\myprime + {\tilde
  e}_{\mu\nu}\myprime$ (again, $\beta v_{\mu\nu}\lr$ cancels) and
(recall) $\rho_\mu=\rho x_\mu$\,; the Fourier transforms ${\tilde
  c}_{\mu\nu}\myprime$ and ${\tilde e}_{\mu\nu}\myprime$ are available
as a byproduct of solving the OZ relation.  With this normalisation
the structure factors obey $S_{\mu\nu}(k)\to\rho_\mu\delta_{\mu\nu}$
as $k\to\infty$.  Two alternative normalisations are due to Aschroft
and Langreth \cite{AL67} and to Hansen and McDonald
\cite{HM06} respectively:
%
\begin{equation}
  \begin{array}{llll}
    S_{\mu\nu}(k)/\rho^{1/2}_\mu \rho^{1/2}_\nu
    &= \delta_{\mu\nu} + \rho x^{1/2}_\mu x^{1/2}_\nu{\tilde h}_{\mu\nu}
    &\quad& \text{(Ashcroft-Langreth)}\,,\\[6pt]
    S_{\mu\nu}(k)/\rho
    &= x_\mu\delta_{\mu\nu} + \rho x_\mu x_\nu{\tilde h}_{\mu\nu}
    && \text{(Hansen-McDonald)}\,.
  \end{array}
\end{equation}
%
The present choice of normalisation is made to facilitate
the calculation of the density-density and charge-charge structure
factors \via\
%
\begin{equation}
  S_{NN}(k) = \frac{\sum_{\mu\nu}S_{\mu\nu}(k)}{\sum_\mu\rho_\mu}\,,\quad
  S_{ZZ}(k) = \frac{\sum_{\mu\nu}z_\mu z_\nu S_{\mu\nu}(k)}{\sum_\mu z_\mu^2\rho_\mu}\,.
  \label{eq:snnszz}
\end{equation}
%
This ensures $S_{NN}$, $S_{ZZ}\to1$ as $k\to\infty$.

\subsection{Thermodynamics}
\label{sec:thermo}
%
The above sections specify the various closures for the
integral equation problem for multicomponent systems.  We provide here
the suite of thermodynamic quantities which can be evaluated from the
converged solution pair $(c_{\mu\nu}', e_{\mu\nu}')$.  Some care is
needed in the case of hard core potentials, where $v_{\mu\nu}=\infty$
for $r<d_{\mu\nu}$.

\subsubsection{Virial pressure}
%
The so-called virial route compressibility factor is
%
\begin{equation}
\frac{\beta p}{\rho}=1-\frac{2\pi}{3}\,{\textstyle \sum_{\mu\nu}}\,
\rho x_\mu x_\nu \int_0^\infty\!\!dr\,r^3\,\frac{\partial(\beta
  v_{\mu\nu})}{\partial r} g_{\mu\nu}(r)
\end{equation}
%
(\cf\ Eq.~(1.2) in \Refcite{Vrbka09}).  To handle the hard core
contribution to this we follow the line of argument presented in
Hansen and McDonald \cite{HM06} and introduce the function
$y_{\mu\nu}=\exp[\beta v_{\mu\nu}]\,g_{\mu\nu}$, which is continuous
through into the hard core region.  The above integral can be written
%
\begin{equation}
\begin{split}
I_{\mu\nu}&=\int_0^\infty\!\!dr\,r^3\,\frac{\partial(\beta
  v_{\mu\nu})}{\partial r} \exp[-\beta v_{\mu\nu}]\,y_{\mu\nu}(r)\\[6pt]
&{}\hspace{8em}{}=-\int_0^\infty\!\!dr\,r^3\,\frac{\partial(\exp[-\beta
    v_{\mu\nu}])}{\partial r} \,y_{\mu\nu}(r)\,.
\end{split}
\end{equation}
%
Let us write
%
\begin{equation}
\exp[-\beta v_{\mu\nu}]=\Theta(r-d_{\mu\nu})\,\exp[-\beta \overline v_{\mu\nu}]
\end{equation}
%
where $\overline v_{\mu\nu}=v_{\mu\nu}$ for $r\ge d_{\mu\nu}$, and we
don't care what it does for $r<d_{\mu\nu}$ except that it should be
\emph{continuous} through $r=d_{\mu\nu}$.  The Heaviside step function
$\Theta(x)=0$ for $x<0$ and $\Theta(x)=1$ for $x\ge0$, with derivative
$\Theta'(x)=\delta(x)$. Then
%
\begin{equation}
\begin{split}
I_{\mu\nu}&=-\int_0^\infty\!\!dr\,r^3\,\Bigl(\delta(r-d_{\mu\nu})
\exp[-\beta \overline v_{\mu\nu}])\\[3pt]
&{}\hspace{8em}{}+\Theta(r-d_{\mu\nu})
\frac{\partial(\exp[-\beta
    \overline v_{\mu\nu}])}{\partial r} \Bigr)\,y_{\mu\nu}(r)\\[9pt]
&=-g_{\mu\nu}(d_{\mu\nu})\,d_{\mu\nu}^3+
\int_{d_{\mu\nu}}^\infty\!\!dr\,r^3\,\frac{\partial(\beta
  v_{\mu\nu})}{\partial r} g_{\mu\nu}(r)
\end{split}
\end{equation}
%
In the second line we have restored $v_{\mu\nu}$ and $g_{\mu\nu}$. 
Thus, finally,
%
\begin{equation}
\frac{\beta p}{\rho}=1
+\frac{2\pi}{3}\,{\textstyle \sum_{\mu\nu}}\,\rho x_\mu x_\nu 
\Bigl[g_{\mu\nu}(d_{\mu\nu})\,d_{\mu\nu}^3 
- \int_{d_{\mu\nu}}^\infty\!\!dr\,r^3\,\frac{\partial(\beta
  v_{\mu\nu})}{\partial r} g_{\mu\nu}(r)\Bigr]
\end{equation}
%
This now contains explicitly the famous contact contribution to the
pressure, and the integral is now restricted to the region external to
the hard core.  In the case of a soft potential (DPD, URPM), the
contact contribution vanishes and the integral extends to $r=0$.

In practice we write the integral contribution to the virial pressure as
%
\begin{equation}
-\frac{2\pi}{3}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^3\,
\frac{\partial({\textstyle \sum_{\mu\nu}}\,
\rho x_\mu x_\nu\,\beta
  v_{\mu\nu})}{\partial r}
+{\textstyle \sum_{\mu\nu}}\,
\rho x_\mu x_\nu t^{\mathrm{V}}_{\mu\nu}
\label{eq:vrpb}
\end{equation}
%
where 
%
\begin{equation}
t^{\mathrm{V}}_{\mu\nu}=-\frac{2\pi}{3}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^3\,
\frac{\partial(\beta v'_{\mu\nu}+\beta v_{\mu\nu}\lr)}{\partial r} 
\,(c'_{\mu\nu}+e'_{\mu\nu})
\end{equation}
%
The first term in \Eqref{eq:vrpb} is the result for $g_{\mu\nu}=1$ and
is the mean-field contribution.  For many applications it is important
to evaluate this analytically, as the individual contributions may
diverge.  Under the SYM condition the contribution from
$v_{\mu\nu}\lr$ to the mean field term vanishes since charge
neutrality implies $\sum_{\mu\nu}\rho x_\mu x_\nu v_{\mu\nu}\lr
\propto v\lr(\sum_\mu \rho_\mu z_\mu)^2 = 0$.  The second term in
\Eqref{eq:vrpb} is the correlation correction (note
$c'_{\mu\nu}+e'_{\mu\nu}=g_{\mu\nu}-1$).

\subsubsection{Compressibility}
%
The compressibility (not to be confused with the above
\emph{compressibility factor}) is
%
\begin{equation}
\frac{\partial(\beta p)}{\partial\rho}=
1-{\textstyle \sum_{\mu\nu}}\,\rho x_\mu x_\nu t^{\mathrm{C}}_{\mu\nu}
\end{equation}
%
where
%
\begin{equation}
t^{\mathrm{C}}_{\mu\nu}={4\pi}\int_0^\infty\!\!dr\,r^2\,(c_{\mu\nu}'(r)
-\beta v_{\mu\nu}\lr)
\end{equation}
%
(\cf\ Eq. (18) in \Refcite{Vrbka09}).  In terms of this, the
isothermal compressibility is $\chi_T=[\rho(\partial
  p/\partial\rho)]^{-1}$.  The contribution from $v_{\mu\nu}\lr$ can
also often be evaluated analytically and vanishes under the SYM
condition by the same argument as above.  Obviously the
compressibility can be integrated along an isotherm to obtain an
alternative expression for the pressure.  This is the so-called
compressibility route to the equation of state.  This expression is
insensitive to the presence of hard cores since they are supposed to
be contained in $v_{\mu\nu}'$ and enter \via\ the direct correlation
function.

\subsubsection{Energy}
%
The excess internal energy density is
%
\begin{equation}
  \beta\uex\equiv
  \frac{\beta \Uex}{V}=2\pi\rho\,{\textstyle \sum_{\mu\nu}}\,\rho x_\mu x_\nu 
\int_{d_{\mu\nu}}^\infty\!\!dr\, r^2\, \beta v_{\mu\nu}(r)\,g_{\mu\nu}(r)
\end{equation}
%
(\cf\ Eq. (2.5.20) in \Refcite{HM06}).  Again this can be split into
a mean field term and a correlation correction,
%
\begin{equation}
{\beta\uex}={2\pi\rho}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^2\,
({\textstyle \sum_{\mu\nu}}\,\rho x_\mu x_\nu\,\beta
  v_{\mu\nu})
+\rho\,{\textstyle \sum_{\mu\nu}}\,
\rho x_\mu x_\nu t^{\mathrm{E}}_{\mu\nu}
\end{equation}
%
where 
%
\begin{equation}
t^{\mathrm{E}}_{\mu\nu}={2\pi}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^2\,
(\beta v'_{\mu\nu}+\beta v_{\mu\nu}\lr)\,(c'_{\mu\nu}+e'_{\mu\nu})\,.
\end{equation}
%
An integration by parts shows that the mean field term here is related
to that for the compressibility factor,
%
\begin{equation}
{2\pi}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^2\,
\beta v_{\mu\nu}
=
-\frac{2\pi\beta v_{\mu\nu}(d_{\mu\nu})\,d_{\mu\nu}^3}{3}
-\frac{2\pi}{3}\int_{d_{\mu\nu}}^\infty\!\!dr\,r^3\,
\frac{\partial(\beta
  v_{\mu\nu})}{\partial r}\,.\label{eq:same}
\end{equation}
%
The equation of state can be derived from the internal energy and this
is the so-called energy route.

\subsubsection{Free energy and chemical potentials}\label{subsec:aex}
%
Last, \emph{only valid for the HNC closure}, are expressions for the
free energy and chemical potentials.  The excess free energy density
can be written as
%
\begin{equation}
  \begin{split}
    \beta\aex\equiv\frac{\beta\Aex}{V}&=\frac{1}{2}\int\!d^3\rvec\,
                        {\textstyle \sum_{\mu\nu}}\,\rho_\mu \rho_\nu
                        (\half h_{\mu\nu}^2-c_{\mu\nu})\\[6pt]
         &\qquad{}+ \frac{1}{2}\int\!\frac{d^3\kvec}{(2\pi)^3}\,
         [\ln\det(I-R\cdot {\tilde C})+ \Tr R\cdot {\tilde C} ]\,.
  \end{split}
\end{equation}
%
The first (real space) term has been simplified from the literature
expressions (\cf\ Eq.~(5.1) in \Refcite{Hir60} or Eq.~(30) in
\Refcite{LB94}), by inserting the HNC closure \Eqref{eq:hncnd}.  The
second (reciprocal space) term is couched in terms of the $R$ and
$\tilde C$ matrices introduced in \Eqref{eq:oz2}: $R$ as a diagonal
matrix with entries $\rho_\mu=\rho_\mu=\rho x_\mu$ and ${\tilde C}$ is
the matrix with entries ${\tilde c}_{\mu\lambda}$.

For the real space part,
$\beta\aex{}^{,1} = (\rho/2)\, {\textstyle \sum_{\mu\nu}}\,\rho x_\mu
x_\nu t^{\mathrm{A}}_{\mu\nu}$ where
%
\begin{equation}
  t^{\mathrm{A}}_{\mu\nu}={4\pi}\int_0^\infty\!\!dr\,r^2\,
  [\half (c_{\mu\nu}'(r)+e_{\mu\nu}'(r)]^2-
  c_{\mu\nu}'(r)+\beta v_{\mu\nu}\lr]\,.
\end{equation}
%
The last term in this is same as that appearing in the compressibility
(to within an overall sign) and often can be evaluated analytically.

We split the reciprocal space part into the separate contributions
from the determinant and the trace.  For the former we use $\ln\det
M=\half\,\sum_i\ln |\lambda_i|^2$ where the sum is over the
eigenvalues $\lambda_i$ of $M=I-R\cdot{\tilde C}$.  Note that this matrix is not
symmetric and the eigenvalues may be complex.  For the latter we write
the trace as $\sum_\mu \rho_\mu ({\tilde c}_{\mu\mu}'-\beta {\tilde
  v}_{\mu\mu}\lr)$ and we note from the definition of the inverse
Fourier transform in \Eqref{eq:ftinv} that $(2\pi)^{-3}\!\int\!
d^3\kvec\, {\tilde v}_{\mu\mu}\lr=v_{\mu\mu}\lr(0)$, \ie\ the long
range potential in real space evaluated at $r=0$.  Thus 
the reciprocal space contribution is computed from
%
\begin{equation}
  \beta\aex{}^{,2}=\frac{1}{4\pi^2}\int_0^\infty\!\!dk\,k^2\,
  \bigl[\half{\textstyle\sum_i}\ln|\lambda_i(k)|^2
    + {\textstyle\sum_\mu} \rho_\mu {\tilde c}_{\mu\mu}'(k)\bigr]
  -\half{\textstyle\sum_\mu} \rho_\mu v_{\mu\mu}\lr(0)\,.
\end{equation}
%
The final excess free energy density is $\aex = \aex{}^{,1} +
\aex{}^{,2}$.  The above expressions are insensitive to the presence
of hard cores.

The corresponding chemical potentials are
%
\begin{equation}
\beta\muex_\mu = 4\pi\,{\textstyle \sum_{\nu}}\,\rho x_\nu 
\int_0^\infty\!\!dr\, r^2\, 
    [\half\,h_{\mu\nu}(r)\,e_{\mu\nu}(r)-c_{\mu\nu}(r)]\quad
    \Bigl(\>{}\equiv\ln\gamma_\mu\>\Bigr)
\end{equation}
%
(\cf\ section II\,C in \Refcite{Vrbka09}).  This 
also defines the activity $\gamma_\mu$.  
We
write this result as $\beta\muex_\mu = {\textstyle \sum_{\nu}}\,\rho x_\nu
t^{\mathrm{M}}_{\mu\nu}$ where
%
\begin{equation}
t^{\mathrm{M}}_{\mu\nu}=
4\pi \int_0^\infty\!\!dr\, r^2\, 
[\half\,(c'_{\mu\nu}+e'_{\mu\nu})(e'_{\mu\nu}+\beta v_{\mu\nu}\lr)
-c'_{\mu\nu}+\beta v_{\mu\nu}\lr]\,.
\end{equation}
%
The contribution from $v_{\mu\nu}\lr$ (the last term) also appears in
the free energy and is computed separately (analytically).  These
expressions are also insensitive to the presence of hard cores.

\section{Potentials}
%
\subsection{DPD potential}
\label{sec:dpd}
%
The code contains routines which specify a multicomponent
potential for dissipative particle dynamics (DPD) with possible
charges.  The potential in this case consists of short range and long
range contributions which can be mapped exactly to the Ng split.  The
short range part is
%
\newcommand{\rc}{r_c}
\newcommand{\lB}{l_{\mathrm{B}}}
\begin{equation}
\beta v'_{\mu\nu}(r)=\left\{\begin{array}{ll}
\frac{1}{2}A_{\mu\nu}(1-r/\rc)^2 & r<\rc\\[3pt]
0 & r \ge \rc
\end{array}\right.\,.
\end{equation}
%
This is a conventional choice, depending on a dimensionless symmetric
repulsion amplitude matrix $A_{\mu\nu}$ and cut off beyond a
distance $\rc$.

\subsubsection{Gaussian charges}
%
Unlike the short range part there is no consensus on the best choice
for the long range interaction although the differences can always be
adsorbed into a redefinition of the short range potential.\footnote{By
  writing $\beta {\tilde v}\lr(k)=4\pi\lB
  k^{-2}(1+k^2\sigma^2/n)^{-n}$, there is a natural hierarchy going
  from Bessel ($n=1$), approximate ($n=2$) and exact ($n=4$) Slater,
  Gaussian ($n\to\infty$).}

The \emph{default} choice of Gaussian charges can be supported for a
range of practical reasons \cite{WVA+13, WV14}.  The Gaussian charges
have size $\sim\sigma$ and a density distribution
$(2\pi\sigma^2)^{-3/2} \exp({-r^2/2\sigma^2})$.  The corresponding
interaction satisfies the SYM condition and is given by
%
\begin{equation}
v_{\mu\nu}\lr(r)=z_\mu z_\nu v\lr(r)\,,\quad
\beta v\lr(r)=\frac{\lB}{r}
\erf\Bigl(\frac{r}{2\sigma}\Bigr)\,,\quad
\beta v\lr(0)=\frac{\lB}{\sigma\sqrt{\pi}}\>,
\label{eq:dpdvlr}
\end{equation}
%
where $\lB$ is the coupling strength (the Bjerrum length).  The
precise definition of $\sigma$ is made to simplify the Fourier
transform (see below).

For use with the above expressions, we need the derivative of both
parts of the potential, and the Fourier transform of the long range
part.  These are
%
\begin{equation}
\frac{\partial(\beta v_{\mu\nu}')}{\partial r}=
\left\{\begin{array}{ll}
-(A_{\mu\nu}/\rc)(1-r/\rc) & r<\rc\\[3pt]
0 & r \ge \rc
\end{array}\right.
\end{equation}
%
\begin{equation}
\frac{\partial(\beta v\lr)}{\partial r}=
-\frac{\lB}{r^2}
\erf\Bigl(\frac{r}{2\sigma}\Bigr)
+\frac{\lB}{r\sigma\sqrt\pi}\,
e^{-r^2/4\sigma^2}\,,
\end{equation}
%
and
%
\begin{equation}
\beta {\tilde v}\lr(k)=\frac{4\pi\lB}{k^2}\times e^{-k^2\sigma^2}\,.
\end{equation}
%
The mean field contributions to the virial route pressure and energy
are given by (see also \Eqref{eq:same})
%
\begin{equation}
{2\pi}\int_0^\infty\!\!dr\,r^2\,
({\textstyle \sum_{\mu\nu}}\,\rho x_\mu x_\nu\,\beta
  v_{\mu\nu})=\frac{\pi\rc^3}{30}\,{\textstyle
  \sum_{\mu\nu}}\,
\rho x_\mu x_\nu A_{\mu\nu}\,.
\end{equation}
%
Because of the SYM condition, the contribution from $v\lr$ vanishes
from this, and also vanishes from the mean field contributions to the
compressibility and chemical potentials (this is true for all the
cases in this section).

This DPD model depends on the dimensionless repulsion amplitude matrix
$A_{\mu\nu}$, the ratios $\lB/\sigma$ and $\sigma/\rc$, and the
dimensionless density $\rho r_c^3$ \cite{WVA+13}.  In these terms
$\rc=1$ is often used as the fundamental length scale.  The ratio
$\lB/\sigma$ plays the role of an inverse effective temperature for
the Coloumbic part of the potential.

\subsubsection{Bessel charges}
%
If Bessel charges are selected the charge distribution becomes
$K_1(r/\sigma) / 2\pi^2\sigma^2 r$ and the expressions change to
\cite{WV14}
%
\begin{equation}
  \beta v\lr(r)=\frac{\lB}{r}(1-e^{-r/\sigma})\,,\quad
  \beta v\lr(0)=\frac{\lB}{\sigma}\,,
\end{equation}
%
\begin{equation}
\frac{\partial(\beta v\lr)}{\partial r}=
-\frac{\lB}{r^2}\Bigl[1-e^{-r/\sigma}\Bigl(1+\frac{r}{\sigma}\Bigr)\Bigr]\,,
\end{equation}
%
\begin{equation}
\beta {\tilde v}\lr(k)=\frac{4\pi\lB}{k^2}\times\frac{1}{1+k^2\sigma^2}\,.
\end{equation}
%
Note that $\sigma$ here is chosen to match the second moment of the
charge distribution for the Gaussian case.

\subsubsection{Linear charge smearing}
%
Another alternative is linear charge smearing proposed by Groot
\cite{Groot03}.  In this case the charge distribution is $(3/\pi
R^3)(1-r/R)$ for $r<R$, vanishing for $r>R$. This gives rise to an
interaction potential in reciprocal space \cite{WV14}
%
\begin{equation}
\beta {\tilde v}\lr(k)=\frac{4\pi\lB}{k^2}
\Bigl(\frac{24-24\cos kR-12kR\sin kR}{k^4R^4}\Bigr)^2\,.
\end{equation}
%
A closed form expression in real space is not available, hence is not
implemented in the code.  This means that the thermodynamics is
\emph{not available} in this case, though the HNC solution goes
through since this only requires the reciprocal space potential.  The
term in large brackets is the charge density in reciprocal space.
Matching the second moment of the charge distribution \cite{WV14}
shows that the equivalent Gaussian smearing length is
$\sigma=R\sqrt{2/15}$.

\subsubsection{Approximate Slater smearing}
%
Another popular alternative is Slater smearing proposed by
Gonz\'alez-Melchor \etal\ \cite{GM+06}.  They use an approximate
expression for the interaction potential, for which we have
%
\begin{equation}
  \beta v\lr(r)=\frac{\lB}{r}\Bigl[1-e^{-2\varbeta r}(1+\varbeta r)\Bigr]\,,
  \quad \beta v\lr(0)=\varbeta\lB\,,
  \label{eq:slaterapprox}
\end{equation}
%
\begin{equation}
  \frac{\partial(\beta v\lr)}{\partial r}=
  -\frac{\lB}{r^2}\Bigl[1-e^{-2\varbeta r}(1+2\varbeta r(1+\varbeta r))\Bigr]\,,
\end{equation}
%
\begin{equation}
  \beta {\tilde v}\lr(k)=\frac{4\pi\lB}{k^2}
  \times\frac{1}{(1+k^2/4\varbeta^2)^2}\,.
\end{equation}
%
These actually corresponds to a charge density $\varbeta^2
e^{-2\varbeta r}/\pi r$ where $\varbeta$ is a smearing parameter, not
to be confused with $\beta=1/\kT$. The equivalent Gaussian smearing
length is here $\sigma=1/\varbeta\sqrt{2}$.

\subsubsection{Exact Slater smearing}
%
The intended charge density in Slater smearing is
$(1/\pi\lambda^3)e^{-2r/\lambda}$.  It turns out that there are exact
expressions for the interaction between such charges \cite{WV14},
%
\begin{equation}
  \beta v\lr(r)=\frac{\lB}{r}\Bigl[1-e^{-2r/\lambda}\Bigl(
    1+\frac{11r}{8\lambda}+\frac{3r^2}{4\lambda^2}
    +\frac{r^3}{6\lambda^3}\Bigr)\Bigr]\,,
  \quad \beta v\lr(0)=\frac{5\lB}{8\lambda}\>,
  \label{eq:slaterexact}
\end{equation}
%
\begin{equation}
  \frac{\partial(\beta v\lr)}{\partial r}=
  -\frac{\lB}{r^2}\Bigl[1-e^{-2r/\lambda}\Bigl(
    1+\frac{2r}{\lambda}+\frac{2r^2}{\lambda^2}
    +\frac{7r^3}{6\lambda^3}+\frac{r^4}{3\lambda^4}\Bigr)\Bigr]\,,
\end{equation}
%
\begin{equation}
  \beta {\tilde v}\lr(k)=\frac{4\pi\lB}{k^2}
  \times\frac{1}{(1+k^2\lambda^2/4)^4}\,.
\end{equation}
%
The equivalent Gaussian smearing length for these exact expressions is
$\sigma=\lambda$.  Originally, Gon\-z\'a\-lez-Melchor
\etal\ \cite{GM+06} set $\varbeta=1/\lambda$ in
\Eqref{eq:slaterapprox}, but setting
$\varbeta=5/(8\lambda)$ matches the values of $\beta v\lr(0)$.

\subsection{URPM and softened URPM potential}
\label{sec:softURPM}
%
Another potential provided by the code at present is for the URPM in
which the unlike pair potential is artificially softened.  The
standard URPM (an equimolar mixture of Gaussian charges) is given by
the DPD potential above with $z_\mu=\pm1$ and $A_{\mu\nu}=0$.  There
are two ways the softened version can be implemented.  The first way
is to work purely in reciprocal space.  It is important to note that
this takes the model away from the SYM condition expressed in
\Eqref{eq:sp}.  In this approach the short range part $v'_{\mu\nu}=0$,
and the long range part is
%
\begin{equation}
\beta v_{11}\lr=v_{22}\lr=\frac{\lB}{r}
\erf\Bigl(\frac{r}{2\sigma}\Bigr)\quad
\beta v_{12}\lr=-\frac{\lB}{r}
\erf\Bigl(\frac{r}{2\sigmap}\Bigr)
\end{equation}
%
where typically $\sigmap>\sigma$ (in the case $\sigmap=\sigma$ we have
the original URPM which is also contained in the DPD potential
described above).  The value at $r=0$, the potential in reciprocal
space, and the derivatives follow \mutmut\ from the DPD case.

In the alternative
approach we retain the SYM condition so that the long range part is
given by \Eqref{eq:dpdvlr} and the short range part is
%
\begin{equation}
\beta v'_{12}\equiv-\beta \Delta v_{12}=-\frac{\lB}{r}
\erf\Bigl(\frac{r}{2\sigmap}\Bigr)
+\frac{\lB\kT}{r}
\erf\Bigl(\frac{r}{2\sigma}\Bigr)
\label{eq:dv12}
\end{equation}
%
(obviously, $v'_{11}=v'_{22}=0$).  We define $\Delta v_{12}$ as the
potential that should be \emph{added} to the softened URPM, to recover
the standard URPM.  This is so that we can use the softened version as
a reference fluid.

For both routes the mean field contributions to the virial route
pressure and energy are non-zero only for the unlike pairs.  We can
evaluate them with $v'_{12}$ given by \Eqref{eq:dv12}
%
\begin{equation}
-\frac{2\pi}{3}\int_0^\infty\!\!dr\,r^3\,
\frac{\partial(\beta
  v'_{12})}{\partial r}={2\pi}\int_0^\infty\!\!dr\,r^2\,
\beta  v'_{12}=2\pi\lB(\sigmap^2-\sigma^2)\,.
\end{equation}
%
If we take the first, purely reciprocal space route, the failure to
satisfy the SYM condition means the compressibility and chemical
potentials need to take account of $v_{\mu\nu}\lr$.  This leads to the
long range contributions 
%
\begin{equation}
{4\pi}\int_0^\infty\!\!dr\,r^2\,\beta
v'_{12}=4\pi\lB(\sigmap^2-\sigma^2)\,.
\end{equation}
%
These should not be included if we follow the second route.

\subsection{RPM and softened RPM potential}
\label{sec:softRPM}
%
Another potential which can be solved in the code is for charged hard
spheres: both the `vanilla' RPM, and a version in which the unlike
pair potential is artificially softened, are provided.  For charged
hard spheres the short range interaction potential can be taken to be
defined by $e^{-\beta v'}=0$ for $r\le\sigma$ and $e^{-\beta v'}=1$
for $r>\sigma$, where $\sigma$ is the hard core diameter (or $\sigma_{\mu\nu}$ in general).

Given this, we might na\"\i{}vely think of using $\lB/r$ as the long
range part of the potential.  Although the divergence at $r\to0$ is
hidden inside the hard core, this choice forces the offset functions
$c'_{\mu\nu}$ and $e'_{\mu\nu}$ to try to compensate, giving rise to
numerical discretisation artefacts (note that the actual direct and
indirect correlation functions, $c_{\mu\nu}$ and $e_{\mu\nu}$, are
non-vanishing and finite inside the hard cores).  Numerically much
better is to chose the long range potential such that it is bounded
inside an inner hard core region.  By far the simplest choice is a
truncated Coulomb law,
%
\begin{equation}
  \beta v\lr(r)=\Bigl\{\begin{array}{ll}
  \lB/\sigma & r\le\sigma\,,\\
  \lB/r & r > \sigma\,,
  \end{array}\quad \beta v\lr(0)=\frac{\lB}{\sigma}\,.\label{eq:rpmhc}
\end{equation}
%
This does least damage to the potential outside the hard cores so we
can continue to use our simple definition of $e^{-\beta v'}$ above.
With this choice, the Fourier transform is easily shown to be
%
\begin{equation}
  \beta {\tilde v}\lr=\frac{4\pi\lB}{k^2}\times\frac{\sin k\sigma}{k\sigma}\,.
\end{equation}
%
The RPM potential in the code is implemented using these
choices. However this choice has a knock-on effect on the
thermodynamics discussed below.

As in the softened URPM, there are two ways the softened RPM can be
implemented.  The first way is to work purely in reciprocal space
which takes the model away from the SYM condition expressed in
\Eqref{eq:sp}.  In this approach we change the long range part between
unlike charges to
%
\begin{equation}
\beta v_{12}\lr=-\frac{\lB}{r} \erf(\kappa r)\,.
\end{equation}
%
The potential in reciprocal space and the derivatives follow
\mutmut\ from the softened URPM case (\ie\ URPM $\sigma\to0$ and
$\sigma'\to1/(2\kappa)$).  Unlike the na\"\i{}ve $\lB/r$ choice for
the vanilla RPM, this softened potential is well-behaved inside the
hard core and can be used without further modification.
  
In the alternative approach we retain the SYM condition so that the
long range part is given as in the vanilla RPM above, and the short range
part is
%
\begin{equation}
  \beta v'_{12}\equiv-\beta\Delta v_{12}=\frac{\lB}{r}\erfc(\kappa r)
  \quad(r>\sigma)\,,
\label{eq:dv12a}
\end{equation}
%
together with $\exp[-\beta v'_{12}]=0$ for $r\le\sigma$.
As in the URPM, we define $\Delta v_{12}$ as the potential that should
be \emph{added} to the softened RPM, to recover the standard RPM.

For both routes the mean field contributions to the virial route
pressure and energy are non-zero only for the unlike pairs.  We can
evaluate them with $v'_{12}$ given by \Eqref{eq:dv12}
%
\begin{eqnarray}
\displaystyle
-\frac{2\pi}{3}\int_\sigma^\infty\!\!dr\,r^3\,
\frac{\partial(\beta
  v'_{12})}{\partial r}=
\pi\lB\Bigl(\frac{\sigma\exp[-\kappa^2\sigma^2]}{\kappa\sqrt{\pi}}
+\Bigl(\frac{1}{2\kappa^2}-\frac{\sigma^2}{3}\Bigr) \erfc(\kappa
\sigma)\Bigr)\,.\\[9pt]
\displaystyle
2\pi \int_\sigma^\infty\!\!dr\,r^2\,\beta v'_{12}=
\pi\lB\Bigl(\frac{\sigma\exp[-\kappa^2\sigma^2]}{\kappa\sqrt{\pi}}
+\Bigl(\frac{1}{2\kappa^2}-\sigma^2\Bigr) \erfc(\kappa
\sigma)\Bigr)\,.
\end{eqnarray}
%
In the limit $\sigma\to0$ these both become $\pi\lB/2\kappa^2$
(\cf\ softened URPM).  In the limit $\kappa\to\infty$ (RPM case) both
vanish.

If we take the first, purely reciprocal space route, the failure to
satisfy the SYM condition means the compressibility and chemical
potentials need to take account of $v_{\mu\nu}\lr$.  This leads to a
contribution for unlike pairs
%
\begin{equation}
  {4\pi}\int_0^\infty\!\!dr\,r^2\,
  \Bigl[\frac{\lB}{r}\erfc(\kappa r)
    +\Bigl(\frac{\lB}{\sigma}
    -\frac{\lB}{r}\Bigr)\Theta(\sigma-r)\Bigr]
  =\frac{\pi\lB}{\kappa^2}-\frac{2\pi\lB\sigma^2}{3}\,.
\end{equation}
%
The second term accounts for the choice of $v\lr$ in the core.

\section{Implementation notes}
%
Most of the code is self-expanatory, given the above mathematical
background.  Correlation functions in the real space domain are
discretised in an array of size $i = 1 \dots n-1$ with a spacing
$\delta_r$ so that $r_i = \delta_r \times i$.  Usually one choses
$\delta_r$ to be a small fraction of the relevant length scale
(\eg\ $\delta_r = 0.01\,\rc$) and $n\delta_r$ to be some large multiple
of the same (\eg\ $n\delta_r\approx 40\,\rc$).  The actual value of $n$
can be chosen to optimise the fast Fourier transform algorithm, for
instance $n=4096$.  The concomitant functions in the reciprocal space
domain are discretised in an array of size $j=1\dots n-1$ with a
spacing $\delta_k \equiv \pi / (n\delta_r)$ so that $k_j = \delta_k
\times j$.  The array size $n-1$ and the value $\delta_k$ are fixed by
the demands of the fast Fourier transform algorithm.

The implementation of the Fourier transforms is as follows.  A
discrete version of \Eqref{eq:fFB} is
%
\begin{equation}
{\tilde h}_j = ({2\pi\delta_r}/{k_j}) \times 
2\, {\textstyle\sum_{i=1}^{n-1}} r_i h_i \sin(\pi i j / n)
\end{equation}
%
where $j = 1\dots n-1$.  The quantity $2\,\sum_{i=1}^{n-1} r_i h_i
\sin(\pi i j / n)$ is computed by calling the \FFTW\ routine
\verb+RODFT00+ on $r_i h_i$.  From the \FFTW\ documenation, this
specific routine works best when the array length is of the form
$2^a3^b5^c7^d11^e13^f-1$ where $e+f$ is either 0 or 1 and the other
exponents are arbitrary (this means $n$ should be of the form of the
first term since the array length is $n-1$).  For the back-transform
the corresponding discrete version of \Eqref{eq:bFB} is
%
\begin{equation}
h_i = ({\delta_k}/{(2\pi)^2 r_i})  \times 
2\,{\textstyle \sum_{j=1}^{n-1}} k_j {\tilde h}_j \sin(\pi i j / n)
\end{equation}
%
The quantity $2\,\sum_{j=1}^{n-1} k_j h_j \sin(\pi i j / n)$ is
computed by calling \verb+RODFT00+ on $k_j h_j$.

Versions $\ge1.7$ of the code can handle an arbitrary number of components.
A species pair array index $(i, j)$ is mapped to a `function' index
$s$ according to the following rule
%
\begin{equation}
  s=\biggl\{\begin{array}{ll}
  i+j(j-1)/2 & \hbox{if $i\le j$\,,}\\[3pt]
  j+i(i-1)/2 & \hbox{if $i>j$\,.}
  \end{array}
  \label{eq:ij}
\end{equation}
%
This is implemented in various forms throughout the code.  The mapping
essentially labels the entries in the upper triangular form of a
symmetric matrix as shown here (where $i$ labels rows, and $j$ labels
columns)
%
\begin{equation}
  \begin{array}{c|ccccc}
    & 1 & 2 & 3 & 4 & \cdots\\
    \hline\\[-12pt]
    1 & 1 & 2 & 4 & 7 & \cdots\\
    2 &   & 3 & 5 & 8 & \cdots\\
    3 &   &   & 6 & 9 & \cdots\\
    4 &   &   &   & \zeroset{10} & \cdots\\
    \cdots & & & & & \cdots
  \end{array}\label{eq:mat}
\end{equation}
%
If necessary, one can recover the array indices $(i, j)$ from the the
index $s$ by the following rule, where $\lfloor\cdots\rfloor$ is the
`floor' function,
%
\begin{equation}
  j = \Bigl\lfloor{\frac{1+\sqrt{8s-7}}{2}}\Bigr\rfloor\,,\quad
  i = n - {j(j-1)}/{2}\,.
\end{equation}
%
This is not used in the code though.  The proof is left as an exercise.

The matrix multiplications in the Ornstein-Zernike relation and elsewhere
are implemented using the standard \verb+MATMUL+ function in
\FORTRAN\ 90.  The matrix inversions are not implemented \perse, but
are rather done by solving the matrix equation $A\cdot X=B$ using a
bespoke subroutine \verb+axeqb_reduce+ in the code which implements a
simplified version of Gauss-Jordan elimination with pivoting (but not
$LDU$ factorisation) \cite{NR92}.  The eigenvalues of the
(non-symmetric) matrix required for the HNC free energy are found by
calling \verb+DGEEV+ in the \LAPACK\ library.

The Ng acceleration scheme is implemented by keeping track of a number
of previous iterations (typically 6), recycling the storage.  Thus the
main functions $c_{\mu\nu}'$ and $e_{\mu\nu}'$ are multidimensional:
the first axis is the spatial discretisation, the second axis is the
species pair index, and the third axis is the Ng trajectory index
(cyclical history array).  The Ng scheme involves also involves
solving $A\cdot X=B$ but this is done with a call to \verb+DSYSV+ in
the \LAPACK\ library (which does implement $LDU$ factorisation).  Now
you might think you could use \verb+DSYSV+ for the matrix operations
in the OZ inversions, or \verb+axeqb_reduce+ in the Ng accelerator,
but for some reason connected to the $LDU$ factorisation this leads to
numerical inaccuracies.  If anyone can shed any light on this please
let me know!

The thermodynamic quantities in section \ref{sec:thermo} are all
evaluated by standard trapezium rules, applied to the integrals given
therein.  

If scripting with \python\ note that the \python\ array index is one
less than the \FORTRAN\ array index.  This conversion is seamlessly
and invisibly implemented in \verb+f2py+.  Also note that the optional
parameters in the \FORTRAN\ function calls can be omitted also in \python,
and the code will gracefully fall back to the documented defaults.

\section{Usage notes}
%
The code is presented as a suite of functions in a \FORTRAN\ 90
module.  Thus it typically needs a driver code, which can be written
in \FORTRAN\ 90 or in \python\ with the use of \verb+f2py+ (see
examples below).  Thermodynamic integration schemes, for example the
compressibility route to the equation of state, are not implemented.
Rather the philosophy has been to focus on robustly solving the basic
integral equation problem.

\subsection{Routines}
%
\begin{itemize}
%
\item\verb+initialise+ -- Allocate all arrays given the number of
  components \verb+ncomp+, the number of real space points \verb+ng+,
  and the length of the Ng trajectory \verb+nps+.  Initialise the
  \FFTW\ plan for fast Fourier transforms.  This \verb+initialise+
  routine is typically called once, at the start, after the values of
  \verb+ncomp+, \verb+ng+, \verb+nps+, and \verb+deltar+ have been
  fixed.
%
\item\verb+dpd_potential(charge_type)+ -- Sets up the DPD potential
  described in section \ref{sec:dpd}.  This routine can be called
  multiple times.  The optional integer parameter indicates whether
  the potential should be initialised with Gaussian charges
  (\verb+1+), Bessel charges (\verb+2+), linear (Groot) charges
  (\verb+3+), Slater charges with approximate interaction (\verb+4+),
  Slater charges with exact interaction (\verb+5+).  The linear case
  (\verb+3+) is currently incomplete in the sense that the
  thermodynamics is not calculated.  If omitted, the
  \verb+charge_type+ defaults to Gaussian.
%
\item\verb+urpm_potential(use_ushort)+ -- Sets up the ultrasoft
  restricted primitive model (URPM) potential with optional softening
  as described in section \ref{sec:softURPM}.  This routine can be
  called multiple times.  It checks that \verb+ncomp = 2+ and forces
  $z_1=1$ and $z_2=-1$.  The optional logical parameter
  \verb+use_ushort+ (default \verb+.false.+) indicates whether the
  potential should be formulated to contain a short range real space
  correction.
%
\item\verb+rpm_potential(use_ushort)+ -- Sets up the restricted
  primitive model (RPM) with optional softening as described in
  section \ref{sec:softRPM}.  This routine can be called multiple
  times.  It checks that \verb+ncomp = 2+ or \verb+ncomp = 3+ and
  forces $z_1=1$, $z_2=-1$, and $z_3=0$ as required (so the third
  component is neutral hard spheres).  The optional logical parameter
  works the same way as in the URPM potential.  The hard core
  diameters are either taken from the array \verb+diam(:, :)+, or if
  these are all zero they are set to a common value \verb+sigma+.  The
  parameter \verb+sigma+ itself is used for the \emph{inner} core
  diameter which features in the definition of $v\lr$ in
  \Eqref{eq:rpmhc}. Thus, by setting \verb+sigma+ and separately
  \verb+diam(:, :)+, one can implement models with dissimilar hard
  sphere radii, or even non-additive hard sphere models.  For
  convenience if \verb+sigma+ was zero (\ie\ the initial default) then it
  is re-set to the minimum of \verb+diam(:, :)+.
%
\item\verb+hs_potential+ -- Sets up hard sphere potentials for an
  arbitrary number of components, with common diameter \verb+sigma+.
  This routine can be called multiple times.
%
\item\verb+hnc_solve+ -- Calculate the HNC solution.  Specifically,
  initialise $c'_{\mu\nu}$ if the logical flag \verb+cold_start+ is
  set, take enough Picard steps to pump-prime the Ng acceleration
  method, then attempt to converge to a solution.  This is the basic
  HNC solver routine and contains appropriate calls to
  \verb+oz_solve+, \verb+hnc_picard+, \verb+hnc_ng+ and
  \verb+conv_test+ below.  A warning is issued if the convergence
  criterion is not met.  If the logical flag \verb+auto_fns+ is also
  set (the default) the three \verb+make_*+ functions are also called,
  resulting in a complete solution to the problem.  The function
  \verb+hnc_solve+ should be called only after \verb+initialise+ and
  the potential has been initialised.  It may be called multiple times
  with different densities for a given potential, or the potential can
  be reset between calls.  In later calls, the existing $c'_{\mu\nu}$
  is used as a starting point unless the logical flag \verb+cold_start+ is
  reset.
%
\item\verb+msa_solve+ -- Calculate the MSA solution, using a similar
  iterative algorithm to the HNC closure. This routine is the basic
  MSA solver and contains appropriate calls to \verb+oz_solve+,
  \verb+msa_picard+, \verb+msa_ng+ and \verb+conv_test+.  A warning is
  issued if the convergence criterion is not met.  If the logical flag
  \verb+auto_fns+ is set (the default) the three \verb+make_*+
  functions are also called, resulting in a complete solution to the
  problem.  Like the HNC, the function \verb+msa_solve+ should be
  called only after \verb+initialise+ and the potential has been
  initialised.  It may be called multiple times with different
  densities for a given potential, or the potential can be reset
  between calls.  In these subsequent calls, the existing
  $c'_{\mu\nu}$ within the hard core regions is used as a starting
  point unless the logical flag \verb+cold_start+ is reset. In all
  cases though, $c'_{\mu\nu}$ outwith the hard core regions is set
  equal to $-\beta v'_{\mu\nu}$ for the current potential. Note that
  the MSA is the same as the RPA if there are no hard cores
  (therefore, use \verb+rpa_solve+ in that case!).
%
\item\verb+rpa_solve+ -- Calculate the RPA solution.  This involves
  only a single round trip through \verb+oz_solve+ from $c'_{\mu\nu} =
  -\beta v'_{\mu\nu}$.  Provided the logical flag \verb+auto_fns+ is
  set (the default), the three \verb+make_*+ functions are also
  called. The function \verb+rpa_solve+ should be called only after
  \verb+initialise+ and the potential has been initialised.  It may be
  called multiple times with different densities for a given
  potential, or the potential can be reset between calls.  The RPA
  closure cannot be used with hard cores (see section \ref{sec:RPA}).
%
\item\verb+save_reference+ -- In the presence of hard cores, EXP
  requires a reference state $\href_{\mu\nu}$, and this routine saves
  the current $h_{\mu\nu}$ solution.  Usually, this reference state
  will be the MSA solution to hard sphere problem without the tail
  potential, however this is not enforced.  See next entry for a note
  on the usage.
%
\item\verb+exp_refine+ -- Implement the EXP refinement.  Use
  \Eqref{eq:exp} to compute a refined $h_{\mu\nu}$ from the current
  solution and the reference $\href_{\mu\nu}$, then call
  \verb+oz_solve2+ to rebuild the complete solution according to
  \Eqref{eq:ozndb}.  Provided the logical flag \verb+auto_fns+ is set
  (the default), the three \verb+make_*+ functions are also called.  A
  call to \verb+initialise+ sets $\href_{\mu\nu}=0$, thus
  \verb+exp_refine+ is suitable for immediate use with soft potentials
  without hard cores.  To use \verb+exp_refine+ for potentials with
  hard cores, calculate and save the reference state, for example
  in \FORTRAN\ 90:
%
\begin{verbatim}
  [ .. remove tail, eg lb = 0.0 or call hs_potential ]
  call msa_solve
  call save_reference
  [ .. restore tail, eg lb > 0.0 ]
  call msa_solve
  call exp_refine
\end{verbatim}
%
Note that \verb+exp_refine+ restores the reference state after
completing so it can be called multiple times with the same reference
state.
%
\item\verb+oz_solve+ -- Solve the OZ relation in the form of
  \Eqref{eq:oznd} to find $e'_{\mu\nu}$ from $c'_{\mu\nu}$.  As a
  byproduct, this generates ${\tilde e}_{\mu\nu}\myprime$ and ${\tilde
    c}_{\mu\nu}\myprime$.
%
\item\verb+oz_solve2+ -- Solve the OZ relation in the form of
  \Eqref{eq:ozndb} to find $e'_{\mu\nu}$ and $c'_{\mu\nu}$ from the
  \emph{reference} state $\href_{\mu\nu}$.  As a byproduct generate
  ${\tilde e}_{\mu\nu}\myprime$ and ${\tilde c}_{\mu\nu}\myprime$.
  Store the result in the first position in the cyclical history
  array.
%
\item\verb+hnc_picard+ -- Solve the HNC closure to find $c'_{\mu\nu}$
  from $e'_{\mu\nu}$, and take a Picard step in the sense of mixing a
  fraction \verb+alpha+ of the new $c'_{\mu\nu}$ with the old
  $c'_{\mu\nu}$.
%
\item\verb+hnc_ng+ -- Solve the HNC closure to find a new
$c'_{\mu\nu}$ given $e'_{\mu\nu}$, but using the Ng acceleration scheme.
%
\item\verb+msa_picard+ -- Solve the MSA closure to find $c'_{\mu\nu}$
  from $e'_{\mu\nu}$, and take a Picard step in the sense of mixing a
  fraction \verb+alpha+ of the new $c'_{\mu\nu}$ with the old
  $c'_{\mu\nu}$.
%
\item\verb+msa_ng+ -- Solve the MSA closure to find a new
$c'_{\mu\nu}$ given $e'_{\mu\nu}$, but using the Ng acceleration scheme.
%
\item\verb+conv_test+ -- Check convergence in iterative schemes by
  calculating the difference (saved in \verb+error+) between the
  current and immediately preceding $c'_{\mu\nu}$.
%
\item\verb+make_pair_functions+ -- from \Eqref{eq:hrs}.  The
  choice to present $h_{\mu\nu}$ rather $g_{\mu\nu}$ is made because
  there may be interest in the asymptotic behaviour of
  $h_{\mu\nu}\to0$ as $r\to\infty$, and it is assumed the user can add
  `1' to get the pair distribution functions.
%
\item\verb+make_structure_factors+ -- from \Eqref{eq:skdef}.
%
\item\verb+make_thermodynamics+ -- everything in section
  \ref{sec:thermo}.
%
\item\verb+hnc_kspace_integrals+ -- specialised calculation of
  reciprocal space integrals in section \ref{subsec:aex} (called
  by \verb+make_thermodynamics+).
%
\item\verb+write_params+ -- Write out basic information.  Often useful
  as a check that the right potential is being solved.
%
\item\verb+write_thermodynamics+ -- Write out all thermodynamic
  quantities, assuming \verb+make_thermodynamics+ has been called.
%
\end{itemize}
%
Note that the three \verb+make_*+ routines are independent and can be
called in any order.

\subsection{User parameters}
%
\begin{itemize}
%
\item\verb+ncomp+ -- The number of component species, $\text{arbitrary}\ge
  1$ (default 1).
%
\item\verb+ng+ -- The real space grid size, $n$, usually a power of 2
  (default 4096).
%
\item\verb+nps+ -- The length of the Ng trajectory, equivalently the
  size of the cyclical history array (usually safe to leave at the
  default 6).
%
\item\verb+npic+ -- The number of Picard steps to pump-prime the Ng
  acceleration scheme (usually safe to leave at the default 6, but
  should be $\ge$ \verb+nps+).
%
\item\verb+maxsteps+ -- maximum number of steps to take before giving up on
convergence (usually safe to leave at default 100).
%
\item\verb+verbose+ -- Logical flag to print solver diagnostics (default
  \verb+.false.+).
%
\item\verb+silent+ -- Suppress printing of warning/error messages
  (default \verb+.true.+).
%
\item\verb+cold_start+ -- Logical flag to indicate whether the main
  solvers should execute a cold start by initialising $c'_{\mu\nu}$.
  This flag is initially set \verb+.true.+, then set \verb+.false.+
  after a main solver call since it is often the case that a
  subsequent problem can re-use the current solution as a good initial
  guess.  The flag can of course be reset \verb+.true.+ by the user at
  any point.
%
\item\verb+start_type+ -- In a HNC cold start initialise
  $c'_{\mu\nu}=0$ (start type 1), $c'_{\mu\nu}=-\beta v'_{\mu\nu}$
  (start type 2), or $c'_{\mu\nu}=\exp(-\beta v'_{\mu\nu})-1$ (start
  type 3, default).  Usually safe to leave at default.
%
\item\verb+auto_fns+ -- Logical flag (default \verb+.true.+) whether
  to call the \verb+make_*+ routines on exit of the main solver
  routines.  Since these routines are relatively inexpensive one would
  usually leave this flag set.
%
\item\verb+deltar+ -- The real space grid spacing $\delta_r$ (default 0.01).
%
\item\verb+alpha+ -- Fraction of new solution in Picard method
  (usually safe to leave at default 0.2).
%
\item\verb+tol+ -- Error tolerance for claiming convergence (usually
  safe to leave at default $10^{-12}$).
%
\item\verb+rc+ -- DPD repulsion range $r_c$ (default 1.0).
%
\item\verb+lb+ -- Coulomb coupling strength $\lB$ (default 0.0)
%
\item\verb+sigma+ -- smearing length $\sigma$ (default 1.0) in the
  case of Gaussian or Bessel charges in the DPD potential and URPM
  potentials; default or inner core diameter $\sigma$ in the case of
  the RPM, or hard sphere potentials.
%
\item\verb+kappa+ -- parameter $\kappa$ (default $-1$) in the
  softened RPM potential.  A negative value is interpreted as
  $\kappa\to\infty$ and generates the vanilla RPM.
%
\item\verb+sigmap+ -- charge size $\sigma'$ (default 1.0), used for
  the softened URPM potential.
%
\item\verb+rgroot+ -- linear (Groot) charge smearing range $R$
  (default 1.0) in the DPD potential (\verb+charge_type = 3+).
%
\item\verb+beta+ -- Slater charge smearing parameter $\varbeta$
  (default 1.0) in the DPD potential with approximate interaction
  (\verb+charge_type = 4+).  One would usually calculate this as
  $\varbeta=1/\lambda$, or (better) $\varbeta=5/(8\lambda)$.
%
\item\verb+lbda+ -- Slater charge smearing length $\lambda$ (default
  1.0) in the DPD potential with exact interaction
  (\verb+charge_type = 5+).  The name is truncated because
  \verb+lambda+ is a reserved word in python.
%
\item\verb+rho(:)+ -- An array of size \verb+ncomp+ where
  \verb+rho(mu)+ is the density $\rho_\mu$.  All entries default to
  zero after a call to \verb+initialise+.
%
\item\verb+z(:)+ -- An array of size \verb+ncomp+ where
  \verb+z(mu)+ is the valence $z_\mu$.  All entries default to
  zero after a call to \verb+initialise+.
%
\item\verb+diam(:, :)+ -- An array of size
  \verb+ncomp+$\times$\verb+ncomp+ where \verb+diam(mu, nu)+ is the
  hard core core diameter $\sigma_{\mu\nu}$ (currently used for RPM
  and HS potentials).  Only elements above the diagonal are used,
  \eg\ \verb+diam(1,2)+, \verb+diam(1,3)+, and \verb+diam(2,3)+.
  Elements below the diagonal are ignored.  All entries default to
  zero after a call to \verb+initialise+.
%
\item\verb+arep(:, :)+ -- An array of size
  \verb+ncomp+$\times$\verb+ncomp+ where \verb+arep(mu, nu)+ is the
  pairwise repulsion amplitude $A_{\mu\nu}$ used in the DPD potential.
  Only elements above the diagonal are used, \eg\ \verb+arep(1,2)+,
  \verb+arep(1,3)+, and \verb+arep(2,3)+.  Elements below the diagonal
  are ignored.  All entries default to zero after a call to
  \verb+initialise+.
%
\end{itemize}

\subsection{Outputs}
%
\begin{itemize}
%
\item\verb+deltak+ -- The reciprocal space grid spacing
  $\delta_k=\pi/(n\delta_k)$, fixed by the needs of the fast Fourier
  transform algorithm.  Available after a call to \verb+initialise+.
%
\item\verb+r(:)+ -- The array \verb+r(i)+ contains
  $r_i=i\times\delta_r$ for $i=1\dots n-1$.  Available after a call to
  \verb+initialise+.
%
\item\verb+k(:)+ -- The array \verb+k(j)+ contains
  $k_j=j\times\delta_k$ for $j=1\dots n-1$.  Available after a call to
  \verb+initialise+.
%
\item\verb+hr(:, :, :)+ -- An array containing the total correlation
  functions $h_{\mu\nu}(r)$ as in \Eqref{eq:hrs}.  Specifically
  \verb+hr(i, mu, nu)+ contains $h_{\mu\nu}(r_i)$ where
  $r_i=i\times\delta_r$ for $i=1\dots n-1$.  The values at $r=0$ can
  be obtained by linear extrapolation \cite{KMP04}.  These arrays are
  available after a call to \verb+make_pair_functions+ which is done
  automatically by the solver routines unless \verb+auto_fns+ is
  turned off.
%
\item\verb+sk(:, :, :)+ -- An array containing structure factors
  $S_{\mu\nu}(k)$ as in \Eqref{eq:skdef}.  Specifically
  \verb+sk(j, mu, nu)+ contains $S_{\mu\nu}(k_j)$ where
  $k_j=j\times\delta_k$ for $j=1\dots n-1$. These arrays are available
  after a call to \verb+make_structure_functions+ which is done
  automatically by the solver routines unless \verb+auto_fns+ is
  turned off.
%
\item\verb+error+ -- The final convergence error estimate for the
  iterative solvers, should be less than \verb+tol+ if converged.
%
\item\verb+return_code+ -- Integer indicates success if
  \verb+return_code = 0+, otherwise an error occurred (see next).
%
\item\verb+error_msg+ -- String (character array) describing the
  problem in the case that \verb+return_code+ is non-zero.  The
  message is printed out anyway unless the \verb+silent+ flag is set.
  Note that in \python\ a returned string appears as fixed
  length binary array and needs massaging to get into a usable
  form, thus for example:
%
\begin{verbatim}
#!/usr/bin/env python3
from oz import wizard as w
...
if w.return_code: 
    msg = str(w.error_msg, 'utf-8').strip()
    raise Exception(msg)
\end{verbatim}
%
\item\verb+closure_name+\,, \verb+model_name+\,, \verb+version+ --
  Strings (character arrays) giving respectively the last-used closure
  (a three-letter acronym), the potential model last used to set the
  potential arrays, and the SunlightHNC version number (added in
  1.10).
%
\end{itemize}
%
The following quantities are available after a call to
\verb+make_thermodynamics+ (which is done automatically unless
\verb+auto_fns+ is turned off).  They are also reported by
\verb+write_thermodynamics+.
%
\begin{itemize}
%
\item\verb+press+\,, \verb+pex+ -- The virial route pressure $\beta p$,
  and excess pressure $\beta \pex=\beta p-\beta\rho$
%
\item\verb+cf_mf+\,, \verb+cf_gc+\,, \verb+cf_xc+ -- The mean-field,
  contact, and correlation contributions respectively, to the virial
  route compressibility factor $\beta p/\rho$ (thus $\beta\pex$ is the sum
  of these multiplied by $\rho$).
%
\item\verb+comp+\,, \verb+comp_xc+ -- The compressibility
  $\partial(\beta p)/\partial\rho$ (not to be confused with the
  above), and the correlation contribution to the same.
%
\item\verb+uex+\,, \verb+uex_mf+\,, \verb+uex_corrn+ -- Excess internal energy
  density $\beta\uex$ (\ie\ excluding kinetic contribution); mean
  field and correlation terms.
%
\end{itemize}
%
The following are calculated only for the HNC closure.
\begin{itemize}
%
\item\verb+aex_rs+\,, \verb+aex_rl+ -- The short and long range
  ($v\lr$ term) pieces of the real space contribution to the excess
  free energy density.
%
\item\verb+aex_kd+\,, \verb+aex_ks+\,, \verb+aex_kl+ -- The
  determinant, and trace broken into short and long range ($v\lr$
  term) pieces, for the reciprocal space contribution to the excess
  free energy density.
%
\item\verb+aex+ -- The total excess free energy density
  $\beta\aex$ (sum of all the above).
%
\item\verb+muex(:)+ -- Chemical potentials, \verb+muex(mu)+ contains
  $\beta\muex_\mu$
%
\item\verb+deficit+ -- The quantity $\aex - \sum_{\mu}\rho_\mu
  \muex_\mu + \pex$ (which should vanish).  This can be viewed as a
  certificate of the numerical discretisation error.
%
\end{itemize}
%
The following internal variables may be useful.
%
\begin{itemize}
%
\item\verb+c(:, :, 1)+\,, \verb+e(:, :, 1)+ -- After a successful
  call to one of the main solver routines, these contain the most
  converged solution, such that for example \verb+c(i, s, 1)+ contains
  $c'_{\mu\nu}(r_i)$ ($=c_{\mu\nu}+\beta v_{\mu\nu}\lr$), where
  \verb+s+ is the species pair index as in \Eqref{eq:ij}.
%
\item\verb+ck(:, :)+\,, \verb+ek(:, :)+ -- Arrays containing the
  Fourier-transformed offset direct and indirect correlation
  functions, such that for example \verb+ck(j, s)+ contains ${\tilde
    c}_{\mu\nu}\myprime(k_j)$ where \verb+s+ is the species pair
  index.
%
\item\verb+h0(:, :)+ -- An array containing the reference state total
  correlation function, such that \verb+h0(i, s)+ contains
  $\href_{\mu\nu}(r_i)$ where \verb+s+ is the species pair index. Used
  by \verb+exp_refine+ and \verb+oz_solve2+.
%
\item\verb+ushort(:, :)+\,, \verb+ulong(:, :)+ --
  Arrays containing the potential, such that \verb+ushort(i, s)+
  contains $\beta v'_{\mu\nu}(r_i)$ and \verb+ulong(i, s)+ contains
  $\beta v\lr_{\mu\nu}(r_i)$, where \verb+s+ is the species index.
%
\item\verb+expnegus(:, :)+ -- A similar array containing $\exp(-\beta
  v'_{\mu\nu})$.
%
\item\verb+dushort(:, :)+\,, \verb+dulong(:, :)+ -- Similar arrays
  containing the derivatives of the potential, $\partial(\beta
  v'_{\mu\nu})/\partial r$ and $\partial(\beta v\lr_{\mu\nu})/\partial
  r$.
%
\item\verb+ulongk(:, :)+ -- Array containing the Fourier transform of
  the long range part of the potential, such that \verb+ulongk(j, s)+
  contains $\beta{\tilde v}\lr_{\mu\nu}(k_j)$, where \verb+s+ is the
  species index.
%
\item\verb+u0(:)+ -- Long range potential at overlap, 
  \verb+u0(mu)+ contains $\beta v_{\mu\mu}\lr(0)$.
%
\item \verb+pi+\,, \verb+twopi+\,, \verb+fourpi+ --- Available for
  convenience ($=\pi, 2\pi, 4\pi$).
%
\end{itemize}

Of these, the potential parts may of course be initialised in the
driver routine, as may the reference state total correlation function
in which latter case \verb+oz_solve2+ acts as a general OZ solver.
Note that the HNC, MSA and RPA solvers require \emph{only} the long
range potential in reciprocal space, $\beta{\tilde v}\lr_{\mu\nu}$
(\verb+ulongk+), and the exponentiated short range potential,
$\exp(-\beta v'_{\mu\nu})$ (\verb+expnegus+), thus a very minimalist
approach can be taken to specifying the potential.  The EXP solver
uses in addition the long range potential in real space, $\beta
v\lr_{\mu\nu}$ (\verb+ulong+).  The full suite of thermodynamics
calculations do require all the remaining potential functions though.

\section{Examples}
%
\subsection{Virial route pressure}
\label{sec:vrp}
%
As a basic example the following minimalist code solves the
HNC closure problem for standard one-component DPD at $\rho=3$ and
$A=25$ using the default grid size and spacing, and prints the virial route pressure,
%
\begin{verbatim}
#!/usr/bin/env python3
from oz import wizard as w
w.initialise()
w.arep[0,0] = A = 25.0
w.dpd_potential()
w.rho[0] = rho = 3.0
w.hnc_solve()
print(f'rho = {rho}, A = {A}')
print('pressure = %0.5f' % w.press)
print('energy density = %0.5f' % w.uex)
print('mu_excess = %0.5f' % w.muex[0])
\end{verbatim}
%
with the result
%
\begin{verbatim}
rho = 3.0  A = 25.0
pressure = 23.56415
energy density = 13.76195
mu_excess =  12.17065
\end{verbatim}
%
The first line imports the \verb+wizard+ of \verb+oz+ and renames it
\verb+w+.  The \verb+oz+ package and the \verb+wizard+ module are
automatically generated by \verb+f2py+.  The actual results from
accurate Monte-Carlo simulations are $p=23.65\pm0.02$,
$e=13.63\pm0.02$, and $\mu^{\mathrm{ex}}=12.25\pm0.05$
so the HNC virial route pressure is out by only
$\approx0.5$\%.

\begin{figure}
\begin{center}
\includegraphics[width=2.5in]{gofr.png}~~~%
\includegraphics[width=2.5in]{sofk.png}
\end{center}
\caption{HNC pair distribution function and structure factor for
  standard DPD at $\rho=3$ and $A=25$.\label{fig:gs}}
\end{figure}

A virtue of \python\ is the large package library, for instance for
graphical output.  The following code uses the \verb+matplotlib+
package to plot the pair distribution function and structure factor
(with a standard normalisation),
%
\begin{verbatim}
#!/usr/bin/env python3
import matplotlib.pyplot as plt
from oz import wizard as w
w.initialise()
w.arep[0,0] = A = 25.0
w.dpd_potential()
w.rho[0] = rho = 3.0
w.hnc_solve()
plt.figure(1)  # This will be g(r)
imax = int(3.0 / w.deltar)
plt.plot(w.r[0:imax], 1.0 + w.hr[0:imax,0,0])
plt.xlabel('$r$')
plt.ylabel('$g(r)$')
plt.figure(2)  # This will be S(k)
jmax = int(25.0 / w.deltak)
plt.plot(w.k[0:jmax], w.sk[0:jmax,0,0]/rho)
plt.xlabel('$k$')
plt.ylabel('$S(k)$')
plt.show()
\end{verbatim}
%
The result is shown in Fig.~\ref{fig:gs} (imagine trying to do this in
pure \FORTRAN!).

\subsection{Compressibility route pressure}
%
The following \python\ routine integrates the compressibility along an
isotherm to calculate the compressibility route pressure,
%
\begin{verbatim}
def cr_press(drho):
    p_xc = prev = 0.0
    n = int(w.rho[0]/drho + 0.5)
    w.cold_start = 1
    for i in range(n):
        w.rho[0] = drho * (i + 1.0)
        w.hnc_solve()
        p_xc = p_xc + 0.5 * drho * (prev + w.comp_xc)
        prev = w.comp_xc
    return w.rho[0] + p_xc
\end{verbatim}
%
Points to note are that a trapezium rule is used for the numerical
integration, and the routine actually integrates the excess
compressibility since the ideal contribution to the pressure is
trivially $\rho$.\\

When the \verb+cr_press+ routine finishes, the
system is in the same state as it started, with the thermodynamics
completely solved.  So for example
%
\begin{verbatim}
w.initialise()
w.arep[0,0] = 25.0
w.dpd_potential()
w.rho[0] = 3.0
print('CR pressure = %0.5f' % cr_press(0.05))
print('VR pressure = %0.5f' % w.press)
print('GP pressure = %0.5f' % (rho + rho*w.muex - w.aex))
\end{verbatim}
%
results in
%
\begin{verbatim}
CR pressure = 22.66623
VR pressure = 23.56415
GP pressure = 23.56441
\end{verbatim}
%
The last here is the (negative) grand potential (GP) density, which in a homogeneous system should equal the pressure; this (numerically) verifies consistency with the virial route.  The pressure calculated by the compressibility route differs by about 4\%
which appears to be typical for HNC for soft potentials.  Also,
further investigation reveals the dependence on $\delta\rho$ is linear
and extrapolates to $p=22.31$ at $\delta\rho\to0$, for $\rho=3$ and
$A=25$.  So the error from discretising the isotherm is here $\alt2$\%.

\subsection{Free energy}
\label{sec:feng}
%
One can calculate the pressure \via\ the energy route, but this really
involves the calculation of the free energy.  One standard way this
can be done is by coupling constant integration.  The basic result can
be derived by differentiating the fundamental expression
$e^{-\beta\Aex} = \int\!d\Omega\,e^{-\beta \phi}$ with respect to a
parameter in the potential $\phi$.  For instance for standard
single-component DPD, ${\partial\aex}/{\partial A}={\uex}/{A}$ where $U$
is the internal energy.  Thus
%
\begin{equation}
{\beta \aex}=\int_0^A\frac{dA'}{A'}\,\beta\uex_{A'}
\label{eq:fint}
\end{equation}
%
can be evaluated along an isochore (constant density line).
The following \python\ routine uses this to compute the excess free
energy
%
\begin{verbatim}
def energy_aex(dA):
    aex_xc = prev = 0.0
    n = int(w.arep[0,0]/dA + 0.5)
    w.cold_start = 1
    for i in range(n):
        w.arep[0,0] = dA * (i + 1.0)
        w.dpd_potential()
        w.hnc_solve()
        curr = w.uex_xc / w.arep[0,0]
        aex_xc = aex_xc + 0.5*dA*(prev + curr)
        prev = curr
    return w.uex_mf + aex_xc
\end{verbatim}
%
Points to note are again that this implements a trapezium rule, and
that the actual integration is carried out for the correlation
contribution for which $\uex/A$ vanishes as $A\to0$.  The mean field
contribution to the internal energy is strictly proportional to $A$,
so it passes unscathed through \Eqref{eq:fint} to contribute additively
to the free energy.

For the HNC specifically, we could also calculate the free energy by
integrating $\muex$ along an isotherm, giving a second routine
%
\begin{verbatim}
def mu_aex(drho):
    aex = prev = 0.0
    n = int(w.rho[0]/drho + 0.5)
    w.cold_start = 1
    for i in range(n):
        w.rho[0] = drho * (i + 1.0)
        w.hnc_solve()
        aex = aex + 0.5*drho*(prev + w.muex[0])
        prev = w.muex[0]
    return aex
\end{verbatim}
%
Both routines leave the system in the same state as it started.  As an
example, to test these,
%
\begin{verbatim}
w.initialise()
w.rho[0] = 3.0
w.arep[0,0] = 25.0
w.dpd_potential()
print('energy aex = %0.5f' % energy_aex(0.1))
print('mu aex     = %0.5f' % mu_aex(0.05))
print('direct aex = %0.5f' % w.aex)
\end{verbatim}
%
(the last of these uses the built-in free energy evaluation for HNC), giving
%
\begin{verbatim}
energy aex = 15.94765
mu aex     = 15.94818
direct aex = 15.94754
\end{verbatim}
%
We see that all three routes agree to the third decimal place, which
is a good test of numerical discretisation error: the absolute HNC
deficit is $2.6\times10^{-4}$ ($1.6\times10^{-5}$ relative).  The above
examples in sections \ref{sec:vrp}--\ref{sec:feng} have been coded up
as \verb+examples.py+.

%\newpage

\subsection{Further examples}

See comments in files for more details.

\begin{itemize}
%
\item\verb+gw_p_compare.py+ -- Plot $(p-\rho) /A\rho^2$ versus $\rho$,
  and compare with data from Fig.~4 of \Refcite{GW97}.  The agreement
  is very good.
%
\item\verb+wsg_fig5_compare.py+ -- Compute the free energy difference
  between a system of $N-1$ A particles and a B particle, and a pure
  system of $N$ A particles, as a function of $\Delta a$, where
  $A_{00}=A_{11}=a$, $A_{01}=a+\Delta a$, and $a=25$.  This concept
  was introduced in \Refcite{WSG01} and is representative of an
  effective $\chi$ value.  Here, the free energy change is calculated
  as $\muex_1-\muex_0$ in a two component problem with the mole
  fraction of the second species set to zero ($\rho_0=3$ and
  $\rho_1=0$).  The results are compared with data from Fig.~5 of
  \Refcite{WSG01}, and with independently done Monte-Carlo simulations
  using trial particle identity swaps.
%
\item\verb+vln_fig2_compare.py+ -- Compute $\log_{10}\gamma_\infty =
  (\muex_1 - \muex_0) / \ln 10$ in the infinite dilution limit for a
  binary DPD fluid, and compare to the empirical fit $0.144 \times
  \Delta a$ as in Fig.~2 in \Refcite{VLN13}.  The fit is derived by
  these authors from Monte-Carlo trial particle (Widom) insertions.
  This is essentially the same calculation as the preceding example
  except that $a=106.5$.  For such a large base repulsion amplitude,
  HNC doesn't converge from a cold start, and this example illustrates
  the trick of warming up the solver by ramping $a$ up in stages.
%
\item\verb+x_dmu_compare.py+ -- In an unpublished extension to the
  above two examples $\muex_1-\muex_0$ is plotted as a function of the
  mole fraction $x$ in a mixture where $\rho_0=(1-x)\rho$ and
  $\rho_1=x\rho$, for $a=25$ and $\Delta a=5$ (one cannot have $\Delta
  a$ too large otherwise an underlying demixing transition causes HNC
  to fail).  The results are compared to data from Monte-Carlo (MC)
  trial particle identity swaps.  The agreement between HNC
  and MC is excellent, including the small deviation from linearity.
  The near-exact linear dependence can be represented by $\chi (1 -
  2x)$ where $\chi$ is the Flory $\chi$-parameter.
%
\item\verb+hm4_fig4-2_compare.py+ solves for the virial and
  compressibility route pressure for the MSA and HNC, for hard
  spheres, and compares with the accurate Carnahan-Starling equation of
  state, as in Fig~4.2 of \Refcite{HM06}.
%
\item\verb+hm4_fig10-3_compare.py+ solves the MSA and HNC for the
  RPM and compares with Monte-Carlo data, as in Fig~10.3 of
  \Refcite{HM06}.
%
\item\verb+attard_fig1_compare.py+ solves the HNC for the RPM to
  generate plots of total correlation functions like that shown in
  \Refcite{Att93}.
%
\item\verb+urpm_oneoff.py+\,, \verb+urpm_scan.py+ solve the ultrasoft
  restricted primitive model (URPM) \cite{WVA+13}, including the
  possibility of a neutral solvent species.  The first can be
  used for one-off calculations of structure and thermodynamics at a
  chosen state point.  The second is intended to scan a range
  of densities to locate the Kirkwood line (transition between pure
  exponential and damped oscillatory behaviour in the tails of
  $h_{\mu\nu}$).  Both routines are driven by command line options,
  via the \verb+argparse+ module.  To see the available options use
  \verb+--help+.
%
\item\verb+rpm_oneoff.py+\,, \verb+hs_oneoff.py+ are similar drivers
  for the RPM and hard sphere systems.  Use \verb+--show+ in these to
  get a graphical output.
%
\item \verb+driver1-4.f90+\,, \verb+driver3.py+ are retained
  for regression testing.
%
\item \verb+fftw_test.f90+ is test code for the fast Fourier transform library.
%
\end{itemize}

%\vspace{36pt}
\appendix
%
\section{Appendix}
%
\subsection{MSA charge-charge structure factor for RPM}
%
The MSA solution for the RPM is known analytically \cite{WL72}.  Some
results are summarised here, for benchmarking the code.  We
focus particularly on the charge-charge structure factor as the closed
form solution is fairly simple.  As in section \ref{sec:softRPM} we
suppose the potential is given by
%
\begin{equation}
  \beta v_{ij}=\Bigl\{\begin{array}{ll}
  \infty & r\le R\\[3pt]
  z_iz_j\lB/r & r>R
  \end{array}
\end{equation}
%
where $z_1=1$ and $z_2=-1$. In this section we shall use $R$ for the
hard core diameter, to match with literature notation.  We shall set
$\rho_1=\rho_2=\rho/2$.  Generically, because of symmetry
($h_{11}=h_{22}$ and $h_{12}=h_{21}$ \etc), the charge-charge
structure factor is given by
%
\begin{equation}
  S_{ZZ}=1+{\textstyle\frac{1}{2}}\rho(\tilde h_{11}-\tilde h_{12})\,.
\end{equation}
%
The symmetry also means that the OZ equation decouples,
%
\begin{equation}
  \tilde h_{11}-\tilde h_{12}=
  \tilde c_{11}-\tilde c_{12} + 
  {\textstyle\frac{1}{2}}\rho(\tilde h_{11}-\tilde h_{12})
  (\tilde c_{11}-\tilde c_{12})\,.
\end{equation}
%
This means that the charge-charge structure factor can be expressed
simply in terms of the Fourier-transformed direct correlation functions,
%
\begin{equation}
  S_{ZZ}=[1-{\textstyle\frac{1}{2}}\rho(\tilde c_{11}-\tilde c_{12})]^{-1}\,.
\end{equation}
%
Now we turn to the MSA solution of the RPM \cite{WL72}.  There are two
ways to calculate the charge-charge structure factor, both of which
lead to the same result.  The first is to Fourier transform the known
solution to the real space direct correlation function.  For this, we
introduce the inverse Debye length $\kappa=(4\pi\lB\rho)^{1/2}$, and
dimensionless version of the same, $x=\kappa R$.  In terms of the
latter we define $B=[1+x-(1+2x)^{1/2}]/x$.  Then the MSA solution for
the direct correlation function difference is
$c_{11}-c_{12}=-({4\lB}/{R})(B-{B^2 r}/{2R})$ for $r\le R$, and
$c_{11}-c_{12}= -{2\lB}/{r}$ for $r > R$ (the latter reflects the MSA
choice $c_{\mu\nu}=-\beta v_{\mu\nu}$ outside the hard core).  The
reason for focusing on the charge-charge structure factor is now
apparent as this difference expression does not involved the
Percus-Yevick solution for uncharged hard spheres, which would
otherwise accompany the individual direct correlation functions.

To do the Fourier transform we write this as
%
\begin{equation}
  \half(c_{11}-c_{12})=-\frac{\lB}{r}+\Bigl[\frac{\lB}{r}
    -\frac{2\lB}{R}\Bigl(B-\frac{B^2 r}{2R}\Bigr)\Bigr]\,\Theta(R-r)
  \label{eq:msa1}
\end{equation}
%
then, \cf\ \Eqref{eq:fFB},
%
\begin{equation}
  \half(\tilde c_{11}-\tilde c_{12})=-\frac{4\pi\lB}{k^2}+
  \frac{4\pi}{k}\int_0^R\!dr\,r\,\sin kr\,\Bigl[\frac{\lB}{r}
  -\frac{2\lB}{R}\Bigl(B-\frac{B^2 r}{2R}\Bigr)\Bigr]\,.
\end{equation}
%
The integral is easily evaluated and inserted in the expression
for the charge-charge structure factor to obtain, after some
simplification,
%
\begin{equation}
  S_{ZZ}(k)=\frac{k^4}{k^4+8q^4+4q^2(k^2-2q^2)\cos kR+8kq^3\sin kR}\,.
  \label{eq:szzrpm}
\end{equation}
%
In this, $q=[(1+2x)^{1/2}-1]/(2R)$.\footnote{Thus $x=2qR(1+qR)$, and
  $B=qR/(1+qR)$.}  \Eqref{eq:szzrpm} can be compared directly with the
charge-charge structure factor reported by \verb+rpm_oneoff.py+.
Complete agreement to within numerical accuracy should be observed.
Note that correctly $S_{ZZ}(k)=k^2/\kappa^2+O(k^4)$ as $k\to0$ (thus
the MSA satisfies the Stillinger-Lovett moment conditions), and
$S_{ZZ}(k)\to1$ as $k\to\infty$.

For completeness, the partner to \Eqref{eq:msa1} is the Percus-Yevick
hard sphere solution \cite{HM06, Wert63}
%
\begin{equation}
  \half(c_{11}+c_{12})=\Bigl\{\begin{array}{ll}
  -\lambda_1-6\eta\lambda_2x-\half\eta\lambda_1x^3 & x<1\,,\\
  0 & x > 1\,,\end{array}
  \label{eq:msa2}
\end{equation}
%
where $x=r/R$, $\eta=\pi\rho R^3/6$ is the packing fraction, and
%
\begin{equation}
\lambda_1=(1+2\eta)^2/(1-\eta)^4\,,\quad
\lambda_2=-(2+\eta)^2/4(1-\eta)^4\,.
\end{equation}
%

An alternative route to the charge-charge structure factor uses
analytic continuation on the Laplace-transformed direct correlation
function.  Specifically, Waisman and Lebowitz \cite{WL72} define $\bar
G(s)=\int_0^\infty dr\,e^{-sr} r\, [g_{11}(r)-g_{12}(r)]$ and give the
result (see also Stell and Sun \cite{SS75})
%
\begin{equation}
  \bar G(s)=-s(2q^2/\pi\rho)[(s^2+2qs+2q^2)e^{sR}-2q^2]^{-1}\,.
\end{equation}
%
From the definitions it is easy to infer
%
\begin{equation}
  S_{ZZ}(k)=1+\frac{2\pi\rho}{k}\imaginary \bar G(-ik)\,.
\end{equation}
%
If the indicated operations are carried out, one indeed recovers
\Eqref{eq:szzrpm}.

\thebibliography{}

\bibitem{WVA+13} P.~B.~Warren, A.~Vlasov, L.~Anton and A.~J.~Masters
  ``Screening properties of Gaussian electrolyte models, with
  application to dissipative particle dynamics'',
  J. Chem. Phys. {\bf138}, 204907 (2013).

\bibitem{Ng74} K.-C. Ng, ``Hypernetted chain solutions for the
  classical one-component plasma up to $\Gamma=7000$'',
  J. Chem. Phys. {\bf61}, 2680--2689 (1974).

\bibitem{HM06} J.-P.~Hansen and I.~R.~McDonald, ``Theory of simple
  liquids, 4th edition'' (Academic Press, Amsterdam, 2013).

\bibitem{KMP04} C.~T.~Kelley and B.~Montgomery Pettitt, ``A fast
  solver for the Ornstein-Zernike equations'',
  J. Comp. Phys. {\bf197}, 491--501 (2004).  Beware typos in this
  paper!  Note, $i$ and $j$ in the main text are $i-1$ and $j-1$ in
  this paper, and $n = N-1$.  This paper also suggests that $c$ and
  $e$ at $r = 0$ are evaluated by simple linear extrapolation, $c_0 =
  2 c_1 - c_2$ and $e_0 = 2 e_1 - e_2$; and $c$ and $e$ at $r =
  n\delta_r$ are zero, $c_n = e_n = 0$.

\bibitem{Vrbka09} L.~Vrbka, M.~Lund, I.~Kalcher, J.~Dzubiella,
  R.~R.~Netz, and W.~Kunz, ``Ion-specific thermodynamics of
  multicomponent electrolytes: A hybrid HNC/MD approach'',
  J.\ Chem.\ Phys.\ {\bf131}, 154109 (2009).

\bibitem{AL67} N.~W.~Ashcroft and D.~C.~Langreth, ``Structure of
  Binary Liquid Mixtures. II. Resistivity of Alloys and the Ion-Ion
  Interaction'', Phys.\ Rev.\ 159, 500 (1967).

\bibitem{Hir60} K.~Hiroike, ``A new approach to the theory of
  classical fluids, II'', Prog.\ Theor.\ Phys.\ {\bf24}, 317 (1960).

\bibitem{LB94} L.~Lue and D.~Blankschtein, ``Proper integral equations
  for interaction-site fluids: exact free-energy expressions'',
  J.\ Chem.\ Phys.\ {\bf100}, 3002 (1994).
  
\bibitem{WV14} P.~B.~Warren and A.~Vlasov ``Screening properties of
  four mesoscale smoothed charge models, with application to
  dissipative particle dynamics'', J.\ Chem.\ Phys.\ {\bf140}, 084904
  (2014).

\bibitem{Groot03}R.~D.~Groot, ``Electrostatic interactions in
  dissipative particle dynamics---simulation of polyelectrolytes and
  anionic surfactants'', J.\ Chem.\ Phys.\ {\bf118}, 11265 (2003).

\bibitem{GM+06}M.~Gonz\'alez-Melchor, E.~Mayoral, M. E. Vel\'azquez
  and J. Alejandre, ``Electrostatic interactions in dissipative
  particle dynamics using the Ewald sums'', J.\ Chem.\ Phys.\ {\bf125},
  224107 (2006).

\bibitem{NR92} W. H. Press, B. P. Flannery, S. A. Teukolsky,
  W. T. Vetterling, ``Numerical recipes in FORTRAN 77: 2nd edition''
  (CUP, Cambridge, 1992).

\bibitem{GW97} R.~D.~Groot and P.~B.~Warren, ``Dissipative particle
  dynamics: bridging the gap between atomistic and mesoscopic
  simulation'', J.\ Chem.\ Phys.\ {\bf107}, 4423 (1997).

\bibitem{WSG01} C.~M.~Wijmans, B.~Smit and R.~D.~Groot, ``Phase
  behavior of monomeric mixtures and polymer solutions with soft
  interaction potentials'', J.\ Chem.\ Phys.\ {\bf114}, 7644 (2001).

\bibitem{VLN13} A.~Vishnyakov, M.~T.~Lee and A.~V.~Neimark,
  ``Prediction of the critical micelle concentration of nonionic
  surfactants by dissipative particle dynamics simulations'',
  J.\ Phys.\ Chem.\ Lett. {\bf4}, 797 (2013).

\bibitem{Att93} P.~Attard, ``Asymptotic analysis of primitive model
  electrolytes and the electrical double layer'', Phys.\ Rev.\ E
  {\bf48}, 3604 (1993).

\bibitem{WL72} E.~Waisman and J.~L.~Lebowitz, ``Mean spherical model
  integral equation for charged hard spheres. II Results'',
  J.\ Chem.\ Phys.\ {\bf56}, 3093 (1972).

\bibitem{Wert63} M.~S.~Wertheim, ``Exact solution of the Percus-Yevick
  integral equation for hard spheres'', Phys.\ Rev.\ Lett.\ {\bf10},
  321 (1963).

\bibitem{SS75} G.~Stell and S.~F.~Sun, ``Generalised mean spherical
  approximation for charged hard spheres: The electrolyte regime'',
  J.\ Chem.\ Phys.\ {\bf63}, 5333 (1975).

\end{document}