# Free energy calculations of patchy particles

###### tags: `Post-doc projects`

---

### Owners (the only one with the permission to edit the main text, others can comment)

Giovanni, Laura, Frank

---

## Notation and definitions
##### Hard sphere ($V_\text{HS}$) and hard disk ($V_\text{HD}$) potential
\begin{equation}
V_\text{HS/HD}(\vec{r}_{ij}) =
\begin{cases}
\infty & \quad , \quad |\vec{r}_j-\vec{r}_i| \leq \sigma \\
0 & \quad , \quad |\vec{r}_j-\vec{r}_i| > \sigma
\end{cases}
\end{equation}
$\sigma$: core's diameter

##### Square shoulder potential
\begin{equation}
V_\text{SS}(\vec{r}_{ij}) =
\begin{cases}
\infty & \quad , \quad |\vec{r}_j-\vec{r}_i| \leq \sigma \\
\varepsilon & \quad , \quad \sigma < |\vec{r}_j-\vec{r}_i| \leq \delta \\
0 & \quad , \quad |\vec{r}_j-\vec{r}_i| > \delta
\end{cases}
\end{equation}
$\sigma$: core's diameter
$\delta$: potential well's diameter

##### Kern-Frenkel potential
\begin{equation}
V_\text{SS}(\vec{r}_{ij}) =
\begin{cases}
\infty & \quad , \quad |\vec{r}_j-\vec{r}_i| \leq \sigma \\
\varepsilon & \quad , \quad \sigma < |\vec{r}_j-\vec{r}_i| \leq \lambda_\text{p} \quad\land\quad
\begin{cases}
\hat{p}_{k_i}\cdot\hat{r}_{ij} > \cos\delta \\
\hat{p}_{k_j}\cdot\hat{r}_{ji} > \cos\delta \\
\end{cases} \quad\text{for some } (k_i,k_j)\\
0 & \quad , \quad |\vec{r}_j-\vec{r}_i| > \lambda_\text{p}
\end{cases}
\end{equation}
$\sigma$: hard core's diameter
$\lambda_\text{p}$: patch maximum interaction range
$\delta$: maximum angle for a bonded patch (in 2D corresponds to half the patch opening angle)
$k_i$: particle $i$'s patch index
$\hat{p}_{k_i}$: orientation versor of $k_i$-th patch of $i$-th particle

##### Ideal gas free energy
\begin{equation}
\frac{\beta F^\text{ex}(N,V)}{N} = \frac{\beta F(N,V)}{N} - \frac{\beta F^\text{id}(N,V)}{N}
\end{equation}
where $\beta=1/(k_BT)$ and $F^\text{id}(N,V)$ is the free energy of an ideal monoatomic gas:
\begin{equation}
\frac{\beta F^\text{id}(N,V)}{N} = \log(\rho\Lambda^D)-1+\frac{1}{2N}\log(2\pi N)
\end{equation}
where $D$ is the dimension of the space, $\rho=N/V$ the particles' density, $\Lambda= \sqrt{h^2/(2\pi m k_BT)}$ is the thermal wavelength, $h$ the Planck's constant.

##### Simulation units
$\sigma$: length
$\varepsilon$: energy (but usually expressed in $k_BT$ units)
$\mu$: mass
$\tau = \sqrt{\frac{\varepsilon}{\mu\sigma^2}}$: time

## Free energy of solids

To calculate $\beta F(V,T)$ I took as a reference point the corresponding ideal Einstein solid and, by following the Frenkel-Ladd scheme, integrated through a thermodynamic path transforming the hamiltonian of the system through the following steps:
\begin{eqnarray}
H_1 &=& H_\text{Ein}^* = P + V_\text{Ein} = \sum_{i}^{N}\frac{p_i^2}{2m} + \lambda_\text{max}\sum_{i}^{N}\frac{(\vec{r}_i-\vec{R}_i)^2}{\sigma^2} \\
H_2 &=& H_{\lambda_\text{max}}^* = P + V_\text{Ein} + V = P + V_\text{Ein} + \sum_{i<j}^{N}V_\text{SS}(\vec{r}_{ij}) \\
H_3 &=& H^* = P + V \\
H_4 &=& H = P + V
\end{eqnarray}
Here I indicate with an asterisk ($*$) the systems in which the center of mass is constrained. $\vec{R}_i$ are the lattice positions of the ideal Einstein solid.

#### Solid particles

In the case of solid particles, like patchy particles, that also have rotational degrees of freedom ($\{\vec{\Omega}_i\}_N$) on top of translational ones, an extra-term needs to be added to (i) the kinetic energy ($P_\Omega$, depending on the angular and inertia momenta), to (ii) the Einstein potential energy (to also constrain the orientation of particles in the solid), and (iii) the energy of the system $V$ actually depends on both the positions ($\{\vec{r}_i\}_N$) and orientations ($\{\vec{\Omega}_i\}_N$) of the particles.

For 2-dimensional identical patchy particles, hence we have only two extra degrees of freedom per particle, the orientation angle of the first patch versor, $\theta_i$, with an associated angular momentum $p_{\theta_i} = I\omega_i = I\dot{\theta_i}$. This ends up in a rotational kinetic energy term:
\begin{equation}
P_\Omega = \sum_{i}^{N}\frac{p_{\theta_i}^2}{2I}
\end{equation}
and in order to obtain the rotational partition function we have to integrate by considering the integral measure:
\begin{equation}
\int\prod_{i}^{N}\frac{\mathrm{d}p_{\theta_i}\mathrm{d}\theta_i}{hn}
\end{equation}
where n is a rotational symmetry factor of the molecule (we integrate over the full $[0,2\pi)$ interval), that in our case is equal to the number of patches. The integral in the space of momenta leads to the free energy term $\beta F_\Omega^{\text{id}} = N\log(2\pi\tilde{\Lambda}/n)$, where we defined $\tilde{\Lambda} = n\sqrt{\hbar^2/(2\pi I k_BT)}$. Since usually in such a kind of calculations, concerning only phase equilibria, temperature is just a constant and we compare everything at the same $T$, we can purposely choose a suitable value for $I$ making $\tilde{\Lambda}=1$. The current definition of $\tilde{\Lambda}$ also allows for the presence of a normalization factor $2\pi$ in the angular configurational integrals, for which the measure becomes:
\begin{equation}
\int_{-\pi}^{\pi}\prod_{i}^{N}\frac{\mathrm{d}\theta_i}{2\pi} \qquad .
\end{equation}

### Free energy of the ideal Einstein solid

With the constraint of the centre of mass' position:
\begin{eqnarray}
\beta F_1 &=& -\log\int \frac{d\vec{r}^Nd\vec{p}^N}{h^{(N-1)D}}\delta(\vec{R}_\text{cm})\delta(\vec{P}_\text{cm})e^{-\beta H_\text{Ein}^*} \nonumber\\
&=& -\frac{1}{2}(N-1)D\log\left( \frac{\pi k_BT}{\lambda_\text{max}}\frac{\sigma^2}{\Lambda^2} \right)
\end{eqnarray}

#### Solid particles

To constrain the orientation of particles I considered two equivalent choices for the angular part of the Einstein potential, one harmonic and the other sinusoidal:
\begin{eqnarray}
V_h(\theta_i) &=& \lambda_\text{max}\left(\frac{n\bar{\theta_i}}{\pi}\right)^2 \qquad,\qquad \bar{\theta_i}=(\theta_i-\theta_i^0 + \frac{\pi}{n})\%\left(\frac{2\pi}{n}\right) - \frac{\pi}{n} \\
V_s(\theta_i) &=& \lambda_\text{max}\sin^2\left(\frac{n(\theta_i-\theta_i^0)}{2}\right)
\end{eqnarray}
where $\theta_i^0$ are the equilibrium orientations in the Einstein crystal[^1].

[^1]: In the definition of the harmonic angular potential $V_h$ I rescaled the angle of a factor $\pi/n$ to make both the considered function shapes to have the same maximum potential barrier, for test purposes while writing the code.

$\DeclareMathOperator\erf{erf}$
The corresponding partition functions are:
\begin{aligned}
\mathcal{Z}_\Omega^h &= \left[\frac{1}{2\pi}\int_{-\pi}^{\pi}d\theta e^{-\beta\lambda_\text{max}\left(\frac{n\bar{\theta}}{\pi}\right)^2}\right]^N = \left[\frac{1}{2}\sqrt{\frac{\pi k_BT}{\lambda_\text{max}}}\erf\left(\sqrt{\frac{\lambda_\text{max}}{k_BT}}\right)\right]^N \\
\mathcal{Z}_\Omega^s &= \left[\frac{1}{2\pi}\int_{-\pi}^{\pi}d\theta e^{-\beta\lambda_\text{max}\sin^2\left(\frac{n\theta}{2}\right)}\right]^N = e^{-N\frac{\beta\lambda_\text{max}}{2}}\left[\mathcal{I}_0\left(\frac{\beta\lambda_\text{max}}{2}\right)\right]^N
\end{aligned}
where $\erf (x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^2}dt$ is the error function, and $\mathcal{I}_k(z)=\frac{1}{\pi}\int_{0}^{\pi}e^{z\cos\phi}\cos(k\phi)d\phi$ are the modified Bessel functions of the first kind (for integer $k$).

After a quick test I realized that there seems to be no substantial difference in the computational efficiency between the two of them.

#### Asymptotic behavior

It is always useful, as a consistency check, to look at the shape of the derivative of the average energy as a function of $\lambda$ and at its asymptotic behavior for very large values of $\beta\lambda$.
In these conditions, the confinement imposed on the system should be very strong, and the only important features are in the end the global minima of the potential energy and their local properties, mainly the range of the confinement, given by the second order positional derivatives.

In that limit we can easily approximate almost all the integrals as gaussian (do not remember the name$\ldots$), it can be shown by naively expanding all the involved functions as Taylor series around the minima of the potential $u(x)$:
\begin{equation}
\int_{x_i}^{x_f}a(x)e^{-\beta\lambda u(x)}\mathrm{d}x \qquad \xrightarrow{\beta\lambda\gg1} \qquad \frac{a^{(i)}(0)}{i!}e^{-\beta\lambda u(0)} \int_\mathbb{R}x^ie^{-\frac{\beta\lambda}{2} u^{\prime\prime}(0)}\mathrm{d}x
\end{equation}
where $i$ is the lower non-zero term in the polynomial expansion of $a(x)$.
By applying this approximation to the calculation of the rotational partition functions we see that it will always be proportional to $1/\sqrt{\beta\lambda u^{\prime\prime}}$. Since for both $V_h(\theta)$ and $V_c(\theta)$ the second derivative is proportional to $n^2$[^2], that term cancel the $n$ term arising from the fact that they both have a period of $2\pi/n$.

[^2]:This can always be changed by rescaling of a suitable factor the potential energy.

In the limit $\beta\lambda\gg1$, the rotational partition functions become:
\begin{aligned}
\mathcal{Z}^{h}_\Omega &\simeq \left[\frac{1}{2}\sqrt{\frac{\pi}{\beta\lambda}}\right]^N \\
\mathcal{Z}^{c}_\Omega &\simeq \left[\frac{1}{\sqrt{\pi\beta\lambda}}\right]^N
\end{aligned}
It is also easy to see that in the harmonic approximation, whatever the shape of the potential, we will always have
\begin{equation}
\left\langle\frac{1}{N}\frac{\mathrm{d}V_\Omega}{\mathrm{d}\lambda}\right\rangle_\lambda \simeq \frac{1}{2\beta\lambda}
\end{equation}
like for the translational degrees of freedom.

### Free energy of the Einstein solid with interacting particles

\begin{equation}
\beta F_2 = \beta F_1 -\log\left\langle e^{-\beta V}\right\rangle_{H_1, \lambda_\text{max}} \approx \beta F_1 + \beta V(\vec{R}^N)
\end{equation}
For a more precise evaluation at finite $\lambda_\text{max}$ (actually for stepwise constant potentials the previous formula can be safely used in any case) I used the following formula, to numerically handle the sum of exponentials inside the logarithm:
\begin{equation}
\beta F_2 = \beta F_1 +\beta\left\langle V\right\rangle -\log\left\langle e^{-\beta (V-\left\langle V\right\rangle)}\right\rangle_{H_1, \lambda_\text{max}}
\end{equation}

### Free energy of the target system with constrained c.o.m.

\begin{equation}
\beta F_3 = \beta F_2 -\beta\int_{0}^{\lambda_\text{max}}d\lambda\left\langle V_\text{Ein}\right\rangle_\lambda
\end{equation}
To evaluate the integral above in a more accurate way I used the Gauss-Legendre quadrature, which is applicable to definite integrals of function well approximated by a linear combination of Legendre polynomials over a symmetric interval $[-1;1]$:
\begin{equation}
\int_{-1}^{1}dx f(x) = \sum_{i=1}^{n}f(x_i)w_i
\end{equation}
where the weights $w_i$ and the nodes $x_i$ are calculated such that the formula is exact for arbitrary polynomials of degree $n_\text{max}\leq2n-1$.
To compute integrals on arbitrary intervals $[a;b]$ a change of variable leads to:
\begin{aligned}
\int_{a}^{b}dy g(y) &= \frac{b-a}{2}\int_{-1}^{1}dx g\left(y(x)\right) = \sum_{i=1}^{n}g\left(\frac{b-a}{2}x_i+\frac{a+b}{2}\right)w_i \\
x &= \frac{2}{b-a}\left(y-\frac{a+b}{2}\right)
\end{aligned}
where $x_i$ and $w_i$ are the same as above. To calculate the free energy integral we conveniently perform the change of variable $\lambda \to \log(\lambda + c)$, where $c$ is a constant suitably chosen, according to the value of $n$, such that the limit values of the derivative of the energy for small $\lambda$ values is similar to the same quantity calculated in the unperturbed solid.
\begin{eqnarray}
\int_{0}^{\lambda_\text{max}}\mathrm{d}\lambda\left\langle V_\text{Ein}\right\rangle_\lambda &=& \int_{\log c}^{\log(\lambda_\text{max}+c)}\mathrm{d}[\log(\lambda+c)](\lambda+c)\left\langle V_\text{Ein}\right\rangle_\lambda = \nonumber\\
&&= \sum_{i=1}^{n} w_i (\lambda_i+c)\left\langle V_\text{Ein}\right\rangle_{\lambda_i} \\
\lambda_i &=& \exp\left[\frac{1}{2}\log\left(\frac{\lambda_\text{max}+c}{c}\right)x_i+\frac{1}{2}\log\left(c\left(\lambda_\text{max}+c\right)\right)\right] - c \nonumber
\end{eqnarray}

#### Solid particles

In the case of 2D patchy particles, I sample the orientation in the interval $[-\pi,\pi)$. An alternative way to do that is to embed in the simulation the symmetry of particles, letting $\theta$ to span over the interval $[-\frac{\pi}{n},\frac{\pi}{n})$. In that case, in the calculation of the Einstein free energy in the following section one should avoid to consider the symmetry factor in the integral. In my case (the first one), to assure the correct sampling of the whole phase space in simulations, which would be hindered for large values of $\beta\lambda$, I added a Monte Carlo move that every $\sim 10$ proposed rotations also let the angle to jump of $\pm\frac{2\pi}{n}$.

### Free energy of the system with unconstrained c.o.m.

\begin{equation}
\beta F_4 = \beta F_3 -\log\left[ \frac{N^{D/2}}{\rho\Lambda^D} \right]
\end{equation}
since dealing with solids of identical particles with single-particle unit cells.

### Total free energy

By summing up all the terms the final free energy per particle is obtained:
\begin{eqnarray}
\frac{\beta F}{N} &=& \frac{1-N}{2N}D\log\left(\frac{\pi k_BT}{\lambda_\text{max}}\right) + \frac{1}{N}\log\left(\rho\sigma^D\right) - \frac{D}{2N}\log N ~ ~ + \\
&-& D\log\left(\frac{\sigma}{\Lambda}\right) + \beta \frac{V(\vec{R}^N)}{N} - \beta\int_{0}^{\lambda_\text{max}}d\lambda\left\langle V_\text{Ein}\right\rangle_\lambda
\end{eqnarray}
In the following I consider $\Lambda = \sigma$.

Then, by subtracting the ideal gas part the excess free energy per particle is obtained:
\begin{eqnarray}
\frac{\beta F^\text{ex}}{N} &=& \frac{1-N}{2N}D\log\left(\frac{\pi k_BT}{\lambda_\text{max}}\right) - \frac{N-1}{N}\log\left(\rho\sigma^D\right) ~ ~ +\\
&-& \frac{D+1}{2N}\log N + 1 -\frac{1}{2N}\log\left(2\pi\right) ~ ~ +\\
&+& \beta \frac{V(\vec{R}^N)}{N} - \beta\int_{0}^{\lambda_\text{max}}d\lambda\left\langle V_\text{Ein}\right\rangle_\lambda
\end{eqnarray}

#### Solid particles ($D=2$)

For patchy particles the additional contribution from the orientational part of the Einstein crystal have to be considered:
\begin{eqnarray}
\frac{\beta F_{\Omega,h}}{N} &=& \log\tilde{\Lambda} + \frac{1}{2}\log\left(\frac{4\lambda_\text{max}}{\pi k_BT}\right) - \log\left[\erf\left(\sqrt{\frac{\lambda_\text{max}}{k_BT}}\right)\right] \\
\frac{\beta F_{\Omega,s}}{N} &=& \log\tilde{\Lambda} + \frac{1}{2}\frac{\lambda_\text{max}}{k_BT} - \log\left[\mathcal{I}_0\left(\frac{\lambda_\text{max}}{2k_BT}\right)\right] =\\
& & \log\tilde{\Lambda} -\log\left[e^{-\frac{\lambda_\text{max}}{2k_BT}}\mathcal{I}_0\left(\frac{\lambda_\text{max}}{2k_BT}\right)\right] \nonumber
\end{eqnarray}
by looking at the second formula, a naive evaluation would lead to a difference among two numbers that can actually be very big, resulting in large numerical errors (*e.g.* the relative error can easily reach $o(10^2)$ for $\beta\lambda_\text{max}\sim200$). Then it is useful to recast the expression, in such a way to have as an argument of the logarithm a numerically manageable number, in this case an exponentially scaled Bessel function (it can be practically evaluated by using several math packages, I used <tt>scipy.special.ive, scipy.special.i0e</tt> in my python script).

#### Finite size corrections

To obtain the free energy in the thermodynamic limit of $N\to\infty$, I need to extrapolate it by considering that the leading terms due to finite size effects are proportional to $\log N/N$ and $N^{-1}$.

[^3]: W. G. Hoover, *J. Chem. Phys.* **49**, 1981–1982 (1968)

From ref.[^3] we can assume that the value of the integral over $\lambda$ in the limit $N\to\infty$ goes like:
\begin{equation}
\frac{\beta}{N}(F_3-F_2) = C + \frac{\log N}{N} + \mathcal{O}\left(N^{-1}\right)
\end{equation}
and considering the other terms $\mathcal{O}\left(\log N/N\right)$ stemming from the center of mass constraint and the ideal gas finite size corrections we obtain that:
\begin{eqnarray}
\frac{\beta F}{N} + \frac{D-2}{2}\frac{\log N}{N} &=& f(\rho,T) + \mathcal{O}\left(N^{-1}\right) \\
\frac{\beta F^\text{ex}}{N} + \frac{D-1}{2}\frac{\log N}{N} &=& f^\text{ex}(\rho,T) + \mathcal{O}\left(N^{-1}\right) \\
\end{eqnarray}
that for sufficiently large $N$ should behave as linear functions of $N^{-1}$.

## Free energy of the liquid phase

To calculate the free energy of the fluid phase I picked up the ideal gas as a reference system, first integrating along a path from infinite to a reference density, then by integrating over temperature, which in the case of Kern-Frenkel potential ($V_\text{KF}$) is equivalent to a continuous variation of the hamiltonian, since the partition function only depends on the factor $\beta\epsilon$.

### Integration along isochores

To find the free energy of the hard disks fluid I started from the known ideal gas free energy and integrated along an isochoric path to the target density, from the 0-density limit:
\begin{equation}
\frac{\beta F(V,T)}{N} = -\frac{\beta}{N}\int_{\infty}^{V}P\mathrm{d}V^\prime = \frac{\beta F^\text{id}(V,T)}{N} + \beta\int_{0}^{\rho}\frac{P_\text{ex}}{{\rho^\prime}^2}\mathrm{d}\rho^\prime
\end{equation}
In this case I independently calculated through Gauss-Legendre quadrature for three values of the density (with this procedure I could obtain a better estimate of the statistical error, probably), and also fitted the excess compressibility factor $\beta P_\text{ex}/\rho$ from the equation of state as a virial expansion (I used a polynomial of degree $20$ and constrained the first three coefficients, known analytically, up to a maximum density $\rho_M = 0.8538$, on 82 points for which I averaged the pressure over 10 independent runs long $t_\text{max}>10^6\tau$; I also checked the values of the first three coefficients by an unconstrained fit, and they agree quite well with the analytical ones). The error on the free energy obtained by the virial fit is tough to estimate, since the covariance matrix of the fit is singular, but they lie within the statistical error bars of the ones obtained by Gauss-Legendre quadrature.

### Integration along isotherms

At a fixed particles density the free energy difference between the target system (2) and one for which the free energy is known (1) (usually hard disks/spheres fluid) is evaluated:
\begin{eqnarray}
H_1 &=& V_\text{HD} \\
H_2 &=& V_\text{HD} + \epsilon W
\end{eqnarray}
then integrating over a parameter continuously turning on $W$:
\begin{equation}
\frac{\beta\Delta F_{12}}{N} = \int_{0}^{\epsilon}\left\langle\frac{\beta}{N}W\right\rangle_\lambda\mathrm{d}\lambda
\end{equation}

#### patchy particles

For patchy particles ($W=V_\text{KF}$) we would have $\beta F_1 = \beta F_\text{HD} + \log\tilde{\Lambda}$, to account for the angular degrees of freedom in the non-interacting patchy particles system. As mentioned in the case of solids, we can choose $I$, according to $T$ and $n$, in such a way to cancel the second addendum. We will have:
\begin{equation}
\frac{\beta\Delta F_{12}}{N} = \int_{0}^{\beta\epsilon}\left\langle\frac{V^\prime_\text{KF}}{N}\right\rangle_{\beta\lambda}\mathrm{d}(\beta\lambda)
\end{equation}
where $V^\prime_\text{KF}$ is the Kern-Frenkel potential with a unitary potential energy well (this average substantially becomes a bonds counting with a negative sign).

## Fluid -- solid coexistence