Hyper-netted chain (HNC) solver for Ornstein-Zernike (OZ) equation

Current version: v1.0 - initial working version

Summary

Implements a python module pyHNC for solving the Ornstein-Zernike (OZ) equation using the hypernetted-chain (HNC) closure, for single-component systems, with soft potentials (no hard cores) such as dissipative particle dynamics (DPD). It uses the FFTW library to do the Fourier transforms, accessed via the pyFFTW wrapper, together with basic NumPy function calls.

The code is intended for rapid prototyping, but is also partly pedagogical with the intent of attempting to capture some of the tacit knowledge involved in undertaking this kind of calculation.

Basic codes in the repository include:

• pyHNC.py : python module implementing the functionality;
• fftw_demo.py : test FFTW for Fourier-Bessel transforms;
• dpd_demo.py : demonstrate the capabilities for standard DPD;
• dpd_eos.py : calculate data for standard DPD equation of state (EoS);
• dpd_gw_compare.py : compare to Groot and Warren, J. Chem. Phys. 107, 4423 (1997).

For more details see extensive comments in the codes, and also the documentation for the parallel SunlightHNC project. The book "Theory of Simple Liquids" by Jean-Pierre Hansen and Ian R. McDonald is foundational -- either the 3rd edition (2006) or the 4th edition (2013).

Simplifications compared to the (faster) SunlightHNC include the fact that hard cores are not implemented, only a single component is assumed, and simple Picard iteration is used rather than the Ng accelerator. The present code is also implemented entirely in python, rather than SunlightHNC which is mostly implemented in FORTRAN 90.

The --sunlight option for dpd_demo.py compares the present implementation to results from SunlightHNC. For this to work, the compiled python interface (*.pyf) and shared object (*.so) dynamically-linked library from SunlightHNC should be made acessible to dpd_demo.py. This can be done for example by copying oz.pyf and oz.*.so from SunlightHNC (available after running make) to the directory containing dpd_demo.py.

Codes related to the nDPD model in Sokhan et al., Soft Matter 19, 5824 (2023);

• ndpd_demo.py : vanilla nDPD, as in the above paper;
• ndpd_rpa.py : implement the RPA and EXP approximations for the EoS;
• ndpd_wca.py : implement Weeks-Chandler-Andersen (WCA) high-temperature approx for EoS;
• ndpd_liquidus.py : estimate the liquidus as the point where p = 0; condor-enabled;

Other codes related to many-body DPD (MBDPD).

• mdpd_hnc.py : various HNC variants for MBDPD; condor-enabled;
• mdpd_dft.py : 'vanilla' DFT for MBDPD;
• mdpd_percus.py : Percus-like DFT for MBDPD;

Other files herein:

• timing.py : for measuring performance of condor DAGMan jobs.
• solute_demo.ipynb : jupyter notebook illustrating solute calculations;
• mu_dpd_hendrikse2025.csv : solute excess chemical potential data digitised from Fig. 1 of Hendrikse et al., Phys. Chem. Chem. Phys. 27, 1554-66 (2025).

What's being solved here?

HNC closure of the OZ equation

What's being solved here is the Ornstein-Zernike (OZ) equation in reciprocal space in the form

$$h(q) = c(q) + \rho\, h(q)\, c(q)$$

in combination with the hypernetted-chain (HNC) closure in real space as

$$g(r)=\exp[-v(r)+h(r)-c(r)]$$

using Picard iteration.

Here $\rho$ is the number density, $v(r)$ is the potential in units of $k_\text{B}T$, $g(r)$ is the pair correlation function, $h(r) = g(r) - 1$ is the total correlation function, and $c(r)$ is the direct correlation function which is defined by the OZ equation. In practice the OZ equation and the HNC closure are written and solved iteratively (see next) in terms of the indirect correlation function $e(r) = h(r) - c(r)$.

Algorithm

Given an initial guess $c(r)$, the solver implements the following scheme (cf SunlightHNC):

• Fourier-Bessel forward transform $c(r) \to c(q)$ ;
• solve the OZ equation for $e(q) = c(q) / [1-\rho\, c(q)]-c(q)$ ;
• Fourier-Bessel back transform $e(q) \to e(r)$ ;
• implement the HNC closure as $c\prime(r)=\exp[-v(r)+e(r)]-e(r)-1$ ;
• replace $c(r)$ by $\alpha\,c\prime(r)+(1-\alpha)\,c(r)$ (Picard mixing step);
• check for convergence by comparing $c(r)$ and $c\prime(r)$ ;
• if not converged, repeat.

Typically this works for a Picard mixing fraction $\alpha = 0.2$, and for standard DPD for example convergence to an accuracy of $10^{-12}$ for a grid size $N_g = 2^{12} = 8192$ with a grid spacing $\Delta r = 0.02$ is achieved with a few hundred iterations (a fraction of a second CPU time).

Initial guess

A suitable initial guess for soft potentials is $c(r) = -v(r)$ ; this is the random-phase approximation (RPA), which for systems without hard cores is equivalent to the mean spherical approximation (MSA). If a similar problem has already been solved then it is often convenient to use the previous solution as the starting point (the solver has been 'warmed up' as it were). This is the default approach.

Structural properties

Once converged, the pair correlation function and static structure factor can be found from:

$$\begin{align} &g(r) = 1 + h(r) = 1 + e(r) + c(r)\,,\\\ &S(q) = 1 + \rho h(q) = 1 + \rho[e(q) + c(q)]\,. \end{align}$$

Thermodynamics

Energy density and virial pressure

Thermodynamic quantities can also now be computed, for example the excess energy density and virial pressure follow from Eqs. (2.5.20) and (2.5.22) in Hansen and McDonald, "Theory of Simple Liquids" (3rd edition) as:

$$\begin{align} &e = \frac{3}{2}\,\rho + 2\pi\rho^2 \int_0^\infty \text{d}r\,r^2\, v(r)\,g(r)\,,\\\ &p = \rho + \frac{2\pi\rho^2}{3} \int_0^\infty \text{d}r\,r^3 f(r)\,g(r) \end{align}$$

where $f(r) = -\text{d}v/\text{d}r$. The first terms here are the ideal contributions, in units of $k_\text{B} T$.

In practice these should usually be calculated with $h(r) = g(r) - 1$, since the mean-field contributions (i.e. the above with $g(r) = 1$) can usually be calculated analytically. Note that in this case an integration by parts shows that the two integrals are actually the same, and are essentially equal to the value of the potential at the origin in reciprocal space:

$$\frac{2\pi\rho^2}{3}\int_0^\infty \text{d}r\,r^3 f(r) =2\pi\rho^2\int_0^\infty \text{d}r\, r^2\, v(r) =\frac{\rho^2}{2}\int \text{d}^3\mathbf{r}\,v(r) =\frac{\rho^2}{2}\,v(q=0)\,.$$

Compressibility

Eq. (2.6.12) in Hansen and McDonald shows that in units of $k_\text{B} T$ the isothermal compressibility satisfies

$$\rho\chi_\text{T} = 1 +4\pi\rho\int_0^\infty \!\!\text{d}r\, r^2\,h(r)$$

where $\chi_\text{T} = - (1/V)\, \partial V/\partial p$. In terms of the EoS, this last expression can be written as $\chi_\text{T}^{-1} = \rho\,\text{d}p/\text{d}\rho$. Further, in reciprocal space the OZ equation (above) can be written as

$$[1+\rho\,h(q)]\,[1-\rho\,c(q)] = 1\,.$$

Employing this at $q = 0$, one therefore obtains

$$\frac{\text{d}p}{\text{d}\rho} = [\rho\chi_\text{T}]^{-1} =1-4\pi\rho\int_0^\infty \!\!\text{d}r\, r^2\, c(r)\,.$$

Given $c(r)$ as a function of density, this can be integrated to find $p(\rho)$. This is known as the compressibility route to the EoS.

Chemical potential

A result peculiar to the HNC is the closed-form expression for the chemical potential given in Eq. (4.3.21) in Hansen and McDonald (I thank Andrew Masters for drawing my attention to this),

$$\frac{\mu}{k_\text{B}T}=\ln\,\rho +4\pi\rho\int_0^\infty\!\!\text{d}r\, r^2\Big[\frac{1}{2}h(h-c)-c\Bigr]\,.$$

Here the reference standard state corresponds to $\rho = 1$. Since the Gibbs-Duhem relation in the form $\text{d}p=\rho\,\text{d}\mu$ can be integrated to find the pressure, this affords another route to the EoS: the chemical potential route.

Free energy and coupling constant integration

It follows from the basic definition of the free energy $F=-\ln\int\text{d}^N\{\mathbf{r}\}\,\exp(-U)$ that $\partial F/\partial\lambda = \langle\partial U/\partial\lambda\rangle$ where $\lambda$ is a parameter in the potential function $U$. We can therefore calculate the free energy from $F = F_0 + \int_0^1\text{d}\lambda\, \langle\partial U/\partial\lambda\rangle_\lambda$. If $\lambda$ is simply a multiplicative scaling, then $\partial U/\partial\lambda = U$ and we have a coupling constant integration scheme, $F = F_0 + \int_0^1\,\text{d}\lambda\,\langle U\rangle_\lambda$ where the indicated average should be taken with the potential energy scaled by a factor $\lambda$. In this scheme $F_0$ is just the free energy of an ideal gas of non-interacting particles since $\lambda\to0$ switches off the interactions.

Since the free energy can be differentiated to find the pressure, this is the basis for the so-called energy route to the EoS. For example, if the free energy density is available as a function of density, $f(\rho)$, the pressure follows from $p = -\partial F/\partial V = \rho\,\text{d}f/\text{d}\rho-f$. It also follows that $[\rho\chi_\text{T}]^{-1} = \text{d}p/\text{d}\rho = \rho\,\text{d}^2f/\text{d}\rho^2$.

The mean-field contribution to the free energy can be calculated immediately since the contribution to the energy density $2\pi\rho^2\int_0^\infty\text{d}r\,r^2\,v(r)$ is independent of $\lambda$ and therefore $\int_0^1\,\text{d}\lambda$ applied to this term trivially evaluates to the same. Furthermore, since this term is $\propto\rho^2$, following the indicated route to the pressure shows that this exact same term appears there too. So the mean-field contribution to the pressure by coupling constant integration is the same as the virial route mean-field pressure.

For the non-mean-field correlation contribution we sketch the algorithm:

• solve the HNC closure of OZ equation for the scaled pair potential $\lambda\,v(r)$ to get $h(r;\lambda)$ ;
• calculate th excess energy $\Delta e(\lambda)=2\pi\rho^2\int_0^\infty\text{d}r\, r^2\,v(r)\,h(r;\lambda)$ with the unscaled pair potential;
• the excess correlation free energy is then the integral $\Delta f=\int_0^1\text{d}\lambda\,\Delta e(\lambda)$
• the excess correlation pressure then follows from $\Delta p = \rho^2\,\text{d}(\Delta f/\rho)/\text{d}\rho$. This should be added to the mean-field contribution to obtain the excess pressure, and the whole added to the ideal contribution to find the total pressure.

In practice a basic trapezium rule can often suffice for the integration step in the above. The derivative with respect to density would usually be computed numerically too.

Solutes

The above methodology can be repurposed to solve also the case of an infinitely dilute solute inside a solvent. To do this we start from the OZ equations for a two-component mixture and specialise to the case where the density of the second component vanishes. In this limit the OZ equations partially decouple in the sense that the solvent case reduces to the above one-component problem, which can be solved as already indicated. The off-diagonal OZ relations become

$$\begin{align} &h_{01}(q) = c_{01}(q)+\rho_{0}\,h_{01}(q)\,c_{00}(q)\,,\\\ &h_{01}(q) = c_{01}(q)+\rho_{0}\,h_{00}(q)\,c_{01}(q)\,. \end{align}$$

The equivalence between the two can be proven from the OZ relation for the solvent. These should be supplemented by the HNC closure for the off-diagonal case

$$g_{01}(r) = \exp[-v_{01}(r)+h_{01}(r)-c_{01}(r)]\,.$$

The second off-diagonal OZ relation can be written as

$$h_{01}(q) = [1+\rho_{0}\,h_{00}(q)]\,c_{01}(q) = S_{00}(q)\,c_{01}(q)\,,$$

where $S_{00}(q)$ is the solvent structure factor. It follows that the solute problem can be solved be re-purposing the exact same algorithm as for the one-component system, replacing the OZ equation step by the assignment,

$$e_{01}(q) = [S_{00}(q)-1]\,c_{01}(q)\,.$$

Applications of this infinitely-dilute solute limit are in the process of being investigated.

FFTW and Fourier-Bessel transforms

The code illustrates how to implement three-dimensional Fourier transforms using FFTW. The starting point is the Fourier transform pair

$$\begin{align} &g(\mathbf{q}) = \int\!\text{d}^3\mathbf{r}\, e^{-i\mathbf{q}\cdot\mathbf{r}}\,f(\mathbf{r})\,,\\\ &f(\mathbf{r}) = \frac{1}{(2\pi)^3}\int\!\text{d}^3\mathbf{q}\, e^{i\mathbf{q}\cdot\mathbf{r}}\,g(\mathbf{q})\,,\\\ \end{align}$$

If the functions have radial symmetry, these reduce to the forward and backward Fourier-Bessel transforms

$$\begin{align} &g(q)=\frac{4\pi}{q}\int_0^\infty\!\! \text{d}r\,r\,f(r)\,\sin\,qr \,,\\\ &f(r)=\frac{1}{2\pi^2\,r}\int_0^\infty\!\! \text{d}q\,q\,f(q)\,\sin\,qr \,. \end{align}$$

From the FFTW documentation, RODFT00 implements

$$Y_k=2\sum_{j=0}^{n-1} X_j\,\sin\Bigl[\frac{\pi(j+1)(k+1)}{n+1}\Bigr]\,,$$

where $n$ is the common length of the arrays $X_j$ and $Y_k$.

To cast this into the right form, set $\Delta r\times\Delta q=\pi/(n+1)$ and assign $r_j=(j+1)\times\Delta r$ for $j=0$ to $n-1$, and likewise $q_k=(k+1)\times\Delta q$ for $k=0$ to $n-1$, so that

$$Y_k=2\sum_{j=0}^{n-1}X_j\,\sin\,q_kr_j\,.$$

For the desired Fourier-Bessel forward transform we can then write

$$g(q_k) = \frac{2\pi\Delta r}{q_k}\times 2\sum_{j=0}^{n-1}r_j\,f(r_j)\,\sin\,q_kr_j\,,$$

with the factor after the multiplication sign being calculated by RODFT00.

The Fourier-Bessel back transform is handled similarly.

On FFTW efficiency

Timing tests (below) indicate that FFTW is very fast when the array length $n$ in the above is a power of two minus one, which doesn't quite seem to fit with the documentation.

Hence, the grid size $N_g=n+1$ in pyHNC is typically a power of two, but the arrays passed to FFTW are shortened to $N_g-1$. Some typical timing results on a moderately fast Intel^® NUC11TZi7 with an 11th Gen Intel^® Core™ i7-1165G7 processor (up to 4.70GHz) support this. For example with 4.2 million grid points

$ time ./fftw_demo.py --ng=2^22 --deltar=1e-3
ng, Δr, Δq, iters = 4194304 0.001 0.0007490140565847857 10
FFTW array sizes = 4194303
real	0m4.087s
user	0m3.928s
sys	0m0.822s

$ time ./fftw_demo.py --ng=2^22+1 --deltar=1e-3
ng, Δr, Δq, iters = 4194305 0.001 0.0007490138780059611 10
FFTW array sizes = 4194304
real	0m10.682s
user	0m9.840s
sys	0m1.505s

$ time ./fftw_demo.py --ng=2^22-1 --deltar=1e-3
ng, Δr, Δq, iters = 4194303 0.001 0.0007490142351636954 10
FFTW array sizes = 4194302
real	0m14.539s
user	0m14.079s
sys	0m1.121s

In the code, FFTW is set up with the most basic FFTW_ESTIMATE planner flag. This may make a difference in the end, but timing tests indicate that with a power of two as used here, it takes much longer for FFTW to find an optimized plan, than it does if it just uses the simple heuristic implied by FFTW_ESTIMATE. Obviously some further investigations could be undertaken into this aspect.

The TL ; DR take-home message here is set $N_g$ to a power of two !

Choice of grid size

From above $\Delta r\times\Delta q=\pi/N_g$ can be inverted to suggest $N_g=\pi/\Delta r\,\Delta q$. Since presumably we want the grid resolution in real space and reciprocal space to be comparable, $\Delta r\simeq\Delta q$, and we want $N_g = 2^r$, this suggests the following table (where $\Delta q$ is computed from $\Delta r$ and $N_g$,

--deltar=0.05 --ng=2^11 (ng=2048 ⇒ Δq ≈ 0.031  )
--deltar=0.02 --ng=2^13 (ng=8192 ⇒ Δq ≈ 0.019  )
--deltar=0.01 --ng=2^15 (ng=32768 ⇒ Δq ≈ 0.0096 )
--deltar=5e-3 --ng=2^17 (ng=131072 ⇒ Δq ≈ 4.79e-3)
--deltar=2e-3 --ng=2^20 (ng=1048576 ⇒ Δq ≈ 1.50e-3)
--deltar=1e-3 --ng=2^22 (ng=4194304 ⇒ Δq ≈ 0.749e-3)

Copying

This program 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.

This program 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 this program. If not, see http://www.gnu.org/licenses/.

Copyright

This program is copyright © 2023-2026 Patrick B Warren (STFC).
Additional modifications copyright © 2025 Joshua F Robinson (STFC).

Contact

Send email to patrick.warren{at}stfc.ac.uk.