# Nested Sampling for Exploring Lennard-Jones Clusters

Lune Maillard, Fabio Finocchi, César Godinho, Martino Trassinelli

► To cite this version:

Lune Maillard, Fabio Finocchi, César Godinho, Martino Trassinelli. Nested Sampling for Exploring Lennard-Jones Clusters. Physical Sciences Forum, 2025, 12, pp.8. 10.3390/psf2025012008. hal-05377063

HAL Id: hal-05377063

https://hal.science/hal-05377063v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Proceeding Paper

# Nested Sampling for Exploring Lennard-Jones Clusters †

Lune Maillard <fnref n="1" />*, Fabio Finocchi <fnref n="1" />, César Godinho <fnref n="2" /> ID and Martino Trassinelli <fnref n="1" />,*, ID

## Abstract

Lennard-Jones clusters, while an easy system, have a significant number of non equivalent configurations that increases rapidly with the number of atoms in the cluster. Here, we aim at determining the cluster partition function; we use the nested sampling algorithm, which transforms the multidimensional integral into a one-dimensional one, to perform this task. In particular, we use the *nested_fit* program, which implements slice sampling as search algorithm. We study here the 7-atom and 36-atom clusters to benchmark *nested_fit* for the exploration of potential energy surfaces. We find that *nested_fit* is able to recover phase transitions and find different stable configurations of the cluster. Furthermore, the implementation of the slice sampling algorithm has a clear impact on the computational cost.

**Keywords:** Lennard-Jones; clusters; nested sampling

## 1. Introduction

Lennard-Jones clusters are relevant models for rare gas atom clusters [1] but the number of configurations increases exponentially with the number of degrees of freedom [2]. They have been widely studied using a variety of methods such as basin-hopping [3], Metropolis Monte Carlo [4], Molecular Dynamics [5] or parallel tempering [6,7].

In this work, we use the nested sampling algorithm [8] to study those clusters. This algorithm is mainly employed to analyse data in fields such as cosmology and astrophysics [9] but is applied more and more in materials science for the exploration of potential energy surfaces, including Lennard-Jones clusters [10]. Here, we utilise *nested_fit* [11–14], which has a different implementation of nested sampling than the one used in Ref. [10]. Indeed, we want to benchmark *nested_fit* for the exploration of potential energy surfaces for classical and, in the future, quantum systems.

We first present the implementation of the nested sampling algorithm for the specific case of Lennard-Jones cluster in Section 2. In Section 3, we study, in particular, clusters of two sizes. In Section 4, we look at the computational cost of our implementation. Finally, we conclude in Section 5.

## 2. Nested Sampling for Lennard-Jones Clusters

For a cluster of *N* atoms, the truncated and shifted Lennard-Jones potential takes the following form [15,16]:

<figure><img src="image_2.png" /></figure>

check for
updates

Academic Editor: Geert Verdoolaege

Published: 13 October 2025

**Citation:** Maillard, L.; Finocchi, F.; Godinho, C.; Trassinelli, M. Nested Sampling for Exploring Lennard-Jones Clusters. *Phys. Sci. Forum* **2025**, *12*, 8. https://doi.org/10.3390/psf2025012008

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Institut des Nanosciences de Paris, Sorbonne Université, CNRS, INSP, F-75005 Paris, France;
fabio@insp.jussieu.fr</footnote>

<footnote n="2">² Laboratory of Instrumentation, Biomedical Engineering and Radiation Physics (LIBPhys-UNL), Department of Physics, NOVA School of Science and Technology, NOVA University Lisbon, 2892-516 Caparica, Portugal; c.godinho@campus.fct.unl.pt</footnote>

<footnote>* Correspondence: lune.maillard@insp.upmc.fr (L.M.); martino.trassinelli@insp.jussieu.fr (M.T.)</footnote>

<footnote>† Presented at the 43rd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, Ghent, Belgium, 1–5 July 2024.</footnote>
$$V = \sum_{\substack{1 \le k < j \le N \\ r_{ij} < r_c}} 4\epsilon \left( \left( \frac{r_0}{r_{ij}} \right)^{12} - \left( \frac{r_0}{r_{ij}} \right)^6 - \left( \frac{r_0}{r_c} \right)^{12} + \left( \frac{r_0}{r_c} \right)^6 \right).$$

We have that $r_{ij}$ is the distance between atoms i and j and that $r_c$ is the cutoff radius, used to remove the interaction at infinite range. The partition function for this system can be written as

$$Z(\beta) = \frac{1}{N! h^{3N}} \int \exp\left(-\beta\left(\frac{1}{2}\sum_{i=1}^{3N} \frac{p_i^2}{m} + V(\mathbf{x})\right)\right) dxd\mathbf{p} \\ = \frac{1}{N! h^{3N}} \int \exp\left(-\beta\left(\frac{1}{2}\sum_{i=1}^{3N} \frac{p_i^2}{m}\right)\right) d\mathbf{p} \int \exp(-\beta(V(\mathbf{x})))d\mathbf{x} = Z_k * Z_c,$$

where $\mathbf{p} = (p_1, \dots, p_{3N})$ is the vector of momenta, $\mathbf{x}$ the vector of all atomic positions, $\beta = \frac{1}{k_B T}$ the inverse temperature and $h$ the Planck constant. The separation in the equation above is possible as the potential $V$ only depends on the atom positions and not on the momenta. The term $Z_k = \frac{1}{N! h^{3N}} \int \exp\left(-\beta\left(\frac{1}{2}\sum_{i=1}^{3N} \frac{p_i^2}{m}\right)\right) d\mathbf{p}$ is easy to compute and gives $(\sqrt{2\pi m/(\beta h^2)})^{3N}/N!$. The term $Z_c = \int \exp(-\beta(V(\mathbf{x})))d\mathbf{x}$ can be rewritten as $Z_c = \int \rho(E) \exp(-\beta E)dE$ where $\rho(E)$ is the density of states [17] and does not in general have an analytical form. We compute $Z_c$ by using a Monte Carlo sampling and, more precisely, the nested sampling algorithm which works in the following way [8,10]:

1. K points, called live points, are uniformly sampled from the entire space.

2. At each iteration $i$, the point $\mathbf{x}_{old}$ associated to the highest energy is removed and replaced by a point $\mathbf{x}_{new}$ with an energy that is strictly lower: $V(\mathbf{x}_{new}) < V(\mathbf{x}_{old})$. The method to find this new point will be presented in Section 2.1. Denoting $E_i = V(\mathbf{x}_{old})$, this point will contribute to the partition function as

$$c_i = w_i e^{-\beta E_i} \text{ with } w_i = \frac{1}{2} \left( \left( \frac{K}{K+1} \right)^{m-1} - \left( \frac{K}{K+1} \right)^{m+1} \right). \quad (1)$$

The term $w_i$ is the approximation of the density of states (DOS) that is evaluated by statistical considerations.

3. This procedure is repeated until the current contribution is small compared to previous contributions: $\log(c_i) - \log(c_{max}) < \delta$ with $c_{max} = \max_{m \le i}(c_m)$ at an inverse temperature $\beta$ chosen by the user. Here, we use $\delta = -10$, the value used in Ref. [14] that gave good results for the harmonic potential.

As the interaction potential does not depend on the temperature, we can compute the partition function at all temperatures by estimating pairs of $(E_i, w_i)$ values during one single exploration of the potential [10] as

$$Z_c(\beta) \approx \sum_i c_i.$$

From the partition function, other properties of the system can be calculated, such as the internal energy

$$U = -\frac{\partial \log(Z)}{\partial \beta} = -\frac{\partial \log(Z_k)}{\partial \beta} - \frac{\partial \log(Z_c)}{\partial \beta} = \frac{3N}{2} k_B T - \frac{\partial \log(Z_c)}{\partial \beta}$$

and the heat capacity

$$C_v = \frac{\partial U}{\partial T} = \frac{3N}{2} k_B - \frac{\partial}{\partial T} \frac{\partial \log(Z_c)}{\partial \beta}.$$
In this paper, we use, in particular, the *nested_fit* program. In the next sections, we present some of the specific implementation of nested sampling in *nested_fit*.

## 2.1. Slice Sampling Real and Slice Sampling Transformed

The method used to find the new point $x_{new}$ in *nested_fit* is slice sampling [18,19]. This method consists of uniformly choosing new exploration points on a slice of the volume defined by the constraint $V(x) < V(x_{old})$. In one dimension, one of the live points is randomly chosen and a slice is built around it until the end points of the slice have an energy higher than $x_{old}$ or are out of the sampling space. A point is then sampled from within the slice and accepted if it verifies the constraint and rejected otherwise. In that case, another point is sampled until an acceptable point is found. In a multidimensional setting, a change in coordinates is first performed to efficiently explore all the parameter spaces, even in the presence of strong correlation. This transformation is performed via the Cholesky transformation of the live points' covariance matrix at the considered step to transform the point coordinates into new coordinates with dimensions ~ $O(1)$ in all directions [19]. The one-dimensional algorithm is then applied recursively to the vectors of $n_{bases}$ randomly generated orthonormal bases (see Ref. [14]). Hence, to compute an update candidate, one live point, corresponding to one N-atom cluster, is moved for $n_{bases} \times 3N$ steps, each one moving all atoms in the cluster. Increasing the value of $n_{bases}$ implies that there is a better decorrelation of samples but it also results in a higher computational cost.

There are two ways of implementing this algorithm:

1. First, the steps can be performed in the transformed space. This will be referred to as *slice sampling transformed* and was used in Ref. [14]. However, the sampled points need to be transformed back to the real space to compute the energy and check that the points are within the bounds. This brings a significant computational overhead to the computation.

2. Second, the steps can be performed in the real space and only the slices are chosen in the transformed space. This will be referred to as *slice sampling real* and was used in Refs. [19,20]. In this case, only the orthonormal bases need to be transformed to the real space.

## 2.2. Parallelisation

The nested sampling algorithm can be easily parallelised. One method that is used to make nested sampling parallel is, instead of searching and substituting one point per iteration, to search *r* new points in parallel and then to replace *r* points at once [21]: the *r* points with the highest energy are removed and replaced by *r* new points with an energy lower than the energies of all removed points. In that case, the compression factor is $\frac{K-r+1}{K}$ instead of $\frac{K}{K+1}$ (see Equation (1)) [17]. However, this parallelisation increases the variance of the results [17]. A different approach is used in Polychord using MPI protocols [19,20]: one primary process is in charge of the sequential replacing of the live points while all other secondary processes are in charge of finding new live points continuously, with the threshold criterion changed by the primary process. The sequential replacing is thus performed by the primary process with the points that are found by the secondary processes. A point found for an iteration with a specific energy constraint is valid for a further iteration if it verifies the new constraint [19]. In this case, the compression factor is still $\frac{K}{K+1}$.

The method we use is very similar to the one used in Polychord. First, the point with the highest energy $x_{max}$ is identified; we note its associated energy $V_{max}$. Then, *r* searches are performed in parallel to find a new point with energy lower than $V_{max}$. We thus obtain a list of *r* new points each with its associated energy. Finally, the points are sequentially
added (without parallelisation) into the set of live points by comparing their energy with $V_{max}$, which evolves with the addition of the new points. There are two possibilities:

* If the energy is lower than $V_{max}$, the point is added to the live points, $x_{max}$ is removed and $V_{max}$ is updated to the new maximum energy of the set of live points. The new value will be used for new points not yet added or rejected.

* If the energy is higher than $V_{max}$, the new point is rejected and $V_{max}$ is not updated.

Hence, for each set of searches performed, between one and *r* points are replaced in our set of live points. Indeed, by definition, all the searches have an energy lower than the initial $V_{max}$ so that at least one point is added.

## 2.3. *Computing the Covariance Matrix*

Another aspect of *nested_fit* that brings a computational overhead is the computation of the covariance matrix and its Cholesky transform at every iteration. To accelerate the program, we have chosen to compute those two matrices only every 0.05 K iterations. This value was chosen as a trade-off between accuracy and computational cost.

# 3. Lennard-Jones Clusters

We now study Lennard-Jones clusters of two sizes: 7 atoms and 36 atoms. In the case of the 7-atom cluster, we compare our results with those found in Ref. [10]. In this work, we use reduced units, i.e., the temperature is in units of $\frac{k_B T}{\epsilon}$ with $\epsilon = 1$. In reduced units, the stopping temperature is 0.01 for the 7-atom cluster and 0.005 for the 36-atom cluster. For the simulation, the atoms are confined in a cubic box of side $L$. To fix the size of the box, we will look at the value of the density $\rho$ in unit of $r_0^{-3}$: $\rho = \frac{N}{(L/r_0)^3}r_0^{-3}$. In all our cases, we take $L = 6$ and use $r_0$ to tune the density to a given desired value. Within reduced units, the results are valid for all atomic species interacting via the Lennard-Jones potential: the specific $(r_0, \epsilon)$ values, corresponding to the atomic species considered, need to be reintroduced instead of the reduced ones. Hence, the specific value of $r_0$ is not important. The density characterises the size of the space the particles evolve in, relative to the number of atoms $N$. We use slice sampling real with the covariance matrix computed every 0.05 K iterations; the procedure is repeated eight times to estimate the uncertainties.

Here, we aim to compute the partition function and derivatives of the clusters, aiming at detecting the possible phase transitions. In an infinite system, a phase transition is indicated by a discontinuity in the heat capacity curve, either in the form of a singularity (delta function) or a critical exponent [22]. However, here, we work in a finite system composed of $N$ particles, in which case a phase transition is not shown by a discontinuity but rather by a peak (or a shoulder) in the heat capacity curve [23].

## 3.1. *Seven Atoms*

First, we consider a cluster made of $N = 7$ atoms with a density of $\rho = 9.27 \times 10^{-3} r_0^{-3}$ (corresponding to $r_0 = 0.659$). The size of the cluster and the density are chosen in order to compare our results with those obtained in Ref. [10]. We use $K = 1000$ live points and $n_{bases} = 5$ so that the samples are sufficiently decorrelated for the comparison. Ref. [10], through the *pymatnest* [10,17,24,25] program, moves the atoms one at a time using rejection Gibbs sampling with $K = 500$ live points. We both use a cutoff radius $r_c$ of $3r_0$. Both heat capacity curves obtained are presented in Figure 1a with examples of solid, liquid and gas configurations. We see that a good agreement with Ref. [10] is obtained. More importantly, *nested_fit* is able to recover both the evaporation and melting phase transitions. We also ran *nested_fit* with $K = 500$ live points but obtained an evaporation peak that was slightly shifted towards higher temperatures compared to the curve from Ref. [10], indicating that there are too few points to obtain convergence. The difference in search method may
explain why we need to double the points compared to Ref. [10]. Indeed, *pymatnest* was developed specifically for the exploration of energy surfaces while *nested_fit* has applications in both data analysis and materials science. The search method used by the latter thus has to work in both cases. Furthermore, in Ref. [10], the Lennard-Jones function was called $2 \times 10^7$ times, while here, it was called $5.4 \times 10^7$ times. There are two main reasons for this increase: first, we use twice as many live points, which should more or less double the number of calls, and second, in slice sampling, the function is called not only to move the points but also to extend the slice. In Figure 1b, we see that the acceptance ratio of a slice sampling search decreases when the volume decreases, showing that it is then more difficult to find new sampling points. The explorations were performed on a 64 bits computer with four CPUs of frequency 3.3 GHz and took around one minute per exploration.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Heat capacity for the Lennard-Jones cluster with 7 atoms. Comparison with the results in Ref. [10]. The configurations shown were obtained with XCrySDen [26]. (b) Evolution of slice sampling acceptance ratio during one of the eight runs.</figcaption></figure>

## 3.2. Thirty-Six Atoms

Second, we consider a cluster composed of $N = 36$ atoms with a density of $\rho = 1.79 \times 10^{-2}r_0^{-3}$ (corresponding to $r_0 = 0.475$). This example is quite challenging as a low-temperature peak is observed, associated with a Mackay-anti-Mackay phase transition [27,28]. Here, we use $n_{bases} = 3$ as the energy function is more costly than in the previous case. The heat capacity obtained with $K = 70,000$ is shown in Figure 2. Again, we can see the presence of a peak corresponding to evaporation and a shoulder corresponding to melting. However, in this case, there is also a peak at lower temperature (around $T = 0.135$) that corresponds to a solid-solid phase transition. From the standard deviation of the heat capacity curve, we can see that the position of this peak is more unstable from one run to the other than the position of the evaporation peak. Indeed, this low-temperature peak is harder to obtain—and was not seen in runs with fewer live points—as we need enough points to fall in each basin so that it is correctly explored. A variation in the proportion of samples in each basin could explain this difference between the runs. We found the positions of the peaks to be similar when more live points—up to $K = 110,000$—were used, indicating that $K = 70,000$ live points are enough to attain convergence. The two solid configurations are shown in Figure 2. The low-temperature peak was also found in Refs. [10,27,28] at around $T = 0.145$, which is in the interval [0.1, 0.15] given by our eight runs. The research in Refs. [10,27,28] was performed at different densities, but this does not affect the position of the solid-solid transition peak. We also tried using $n_{bases} = 1$ but we observed a shift in the evaporation peak: the number of bases used is most likely too small to decorrelate the samples. The explorations were performed on a 64 bits computer with 64 CPUs of frequency 2.0 GHz and took around 43 hours per exploration.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Heat capacity for the Lennard-Jones cluster with 36 atoms. The configurations shown were obtained with XCrySDen [26]. The solid line is the mean obtained from the eight runs and the dashed curves correspond to the first standard deviation.</figcaption></figure>

**4. Program Profiling and Optimisation**

In this section, we look at the impact of the different developments made to `nested_fit` on the computational resources needed. For that, we take the 36-atom cluster for its high number of parameters (108). First, we look at the result of profiling—performed with `valgrind` [29,30] and `gprof2dot` [31]—to see which parts of the algorithm are the most computationally expensive. We use $K = 500$ and $n_{bases} = 1$. We consider three cases: slice sampling transformed with the calculation of the covariance matrix at each iteration, slice sampling transformed with the calculation of the covariance matrix every 0.05 K iterations and slice sampling real with the calculation of the covariance matrix every 0.05 K iterations. The results of the profiling are presented in Figures 3, 4 and 5, respectively. For the first case, version 4.5.5 of `nested_fit` was used and version 4.6.0 for the two other cases.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Profiling for slice sampling transformed with the calculation of the covariance matrix at each iteration.</figcaption></figure>

The purpose of each function is explained in Table 1. Only the relevant boxes are represented. The last three lines in a box correspond to (1) the percentage of time spent in this function and its children function, (2) the percentage of time spent in this function alone and (3) the number of times the function was called. The arrows linking one box to another represent the first (parent) function calling the second (child) with the two numbers being the percentage of time the child function transfers to its parent and the number of times the parent function calls the child function. A box is only represented if at least 5% of the time is spent in the function.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Profiling for slice sampling transformed with the calculation of the covariance matrix every 0.05K iteration.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Profiling for slice sampling real with the calculation of the covariance matrix every 0.05 K iteration.</figcaption></figure>

First, looking at Figures 3 and 4, we can see that, in the case of slice sampling transformed, over 50% of the time is spent returning to the real space to calculate the likelihood and check the boundaries (i.e. in the part_like_sub and test_bnd_sub functions), which is not needed in slice sampling real. Thus, more time is spent changing space than calculating the function studied, leading to a computational overhead for slice sampling transformed. Furthermore, we can see that when the covariance function is computed at each iteration (Figure 3), this calculation takes around 20% of the time whereas when it is only computed every 0.05K iterations (here every 25 iterations) (Figures 4 and 5), the box does not appear, which means that less than 5% of the time is spent on this particular task. Finally, we see that, for slice sampling transformed with the covariance matrix computed at each iteration (Figure 3), around 12% of the time is spent calculating the function studied. This number rises up to around 55% (4.4 times more) when using slice sampling real and computing the covariance function every 0.05K iteration (Figure 5). In both cases, the number of times the energy function was called is roughly the same (around 580,000 times). This shows that, by performing slice sampling in the real space rather than in the transformed space and not computing the covariance matrix at every iteration, we removed calculations that are not needed and that bring a significant computational overhead: the computational time is reduced by a factor of 2.8, as can be deduced from Table 2, which gives the computational time of all three cases for 20,000 iterations.

For the impact of the parallelisation on the computation time, we also look at the time required to run *nested_fit* with and without parallelisation using slice sampling transformed with $K = 1000$ and $n_{bases} = 1$. The covariance matrix is computed every iteration. The parallelised case over 64 cores is around 21 times faster than the non-
parallelised case (the non-parallelised case takes 85 min while the parallelised case takes 4 min). A possible explanation for this dividing factor being smaller than the number of cores used is that the time gained by simultaneously searching 64 new points may be partially counterbalanced by the time taken for the sequential adding of all the points, some of which may be rejected due to the constraint changing with each added point.

<table><caption>**Table 1.** Purpose of the functions shown in Figures 3–5.</caption><tr><td><code>_mod_search_new_point_MOD_slice王</code></td><td>Performs one iteration of nested sampling with slice sampling real</td></tr><tr><td><code>_mod_search_new_point_MOD slice</code></td><td>Performs one iteration of nested sampling with slice sampling transformed</td></tr><tr><td><code>_mod_ likelihood_MOD_loglikelihood</code></td><td>Calls the function to calculate (here the Lennard-Jones potential)</td></tr><tr><td><code>_mod_search_new_point_MOD mat cov</code></td><td>Computes the covariance matrix</td></tr><tr><td><code>_mod_search_new_point_MOD_base o_n</code></td><td>Creates an orthonormal basis</td></tr><tr><td><code>_mod_search_new_point_MOD_part like sub</code></td><td>Calculates the function for a point in the transformed space</td></tr><tr><td><code>_mod_search_new_point_MOD_test bnd sub</code></td><td>Checks if a point in the transformed space is inside the sampling space</td></tr></table>

<table><caption>**Table 2.** Time taken by all three cases to perform 20,000 iterations without parallelisation.</caption><thead><tr><th></th><th>Version 4.5.5-Slice<br>sampling transformed</th><th>Version 4.6.0-Slice<br>sampling transformed</th><th>Version 4.6.0-Slice<br>sampling real</th></tr></thead><tbody><tr><th>Time (s)</th><td>442.19</td><td>367.57</td><td>158.13</td></tr></tbody></table>

## 5. Conclusions

In this work, we have applied the nested sampling algorithm to study classical Lennard-Jones clusters.

First, we have studied clusters of two sizes. As for the 7-atom cluster, we have seen that we were able to recover phase transitions by comparing our work with that in Ref. [10]. We saw that we required more live points and more likelihood calls than those in Ref. [10], likely due to the different choices of method to find a new point, which is an inconvenience for a more computationally expensive energy function. We then studied the 36-atom cluster, which presents a solid-solid phase transition at low temperature that *nested_fit* was able to recover provided that enough live points were used.

We then looked at the impact of different implementation types on the computational resources. We saw that using slice sampling transformed brought a consequent computational overhead through the need to transform back to the real space to compute the energy. Slice sampling real, which works directly in the real space, does not need to perform this task and therefore spends more than triple the time calculating the energy compared to slice sampling transformed, even though the function was called around the same number of times. We also saw that, using parallelisation on 64 cores, we gained a factor of 21 in the computation time. All of these developments allowed us to study bigger and more complex clusters that require a higher number of live points. For example, one run for the 36-atom Lennard-Jones clusters with $K = 70,000$ took around 43 h to run with the parallelisation.

In the future, we would like to study quantum Lennard-Jones clusters which are by about an order of magnitude more computationally expensive than the classical ones.
**Author Contributions:** Conceptualisation, L.M., M.T. and F.F.; methodology, L.M., M.T. and F.F.; software, M.T., L.M. and C.G.; validation, L.M.; formal analysis, L.M.; investigation, L.M.; resources, F.F.; data curation, L.M. and M.T.; writing—original draft preparation, L.M.; writing—review and editing, F.F. and M.T.; visualisation, L.M.; supervision, M.T. and F.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Raw data can be provided on demand.

**Acknowledgments:** Work realised with the support of the Sorbonne Center for Artificial Intelligence 579—Sorbonne University—IDEX SUPER 11-IDEX-0004. The authors thank Philipe Depondt, Simon Huppert and Julien Salomon for helpful discussions.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Ashcroft, N.W.; Mermin, N.D. *Solid State Physics*; Holt- Rinehart and Winston: New York, NY, USA, 1976.

2. Doye, J.P.K.; Miller, M.A.; Wales, D.J. Evolution of the potential energy surface with size for Lennard-Jones clusters. *J. Chem. Phys.* **1999**, *111*, 8417–8428. [CrossRef]

3. Wales, D.J.; Doye, J.P.K. Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms. *J. Phys. Chem.* **1997**, *101*, 5111–5116. [CrossRef]

4. Quirke, N.; Sheng, P. The melting behavior of small clusters of atoms. *Chem. Phys. Lett.* **1984**, *110*, 63–66. [CrossRef]

5. Beck, T.L.; Berry, R.S. The interplay of structure and dynamics in the melting of small clusters. *J. Chem. Phys.* **1988**, *88*, 3910–3922. [CrossRef]

6. Neirotti, J.P.; Calvo, F.; Freeman, D.L.; Doll, J.D. Phase changes in 38-atom Lennard-Jones clusters. I. A parallel tempering study in the canonical ensemble. *J. Chem. Phys.* **2000**, *112*, 10340–10349. [CrossRef]

7. Calvo, F.; Neirotti, J.P.; Freeman, D.L.; Doll, J.D. Phase changes in 38-atom Lennard-Jones clusters. II. A parallel tempering study of equilibrium and dynamic properties in the molecular dynamics and microcanonical ensembles. *J. Chem. Phys.* **2000**, *112*, 10350–10357. [CrossRef]

8. Skilling, J. Nested Sampling. *AIP Conf. Proc.* **2004**, *735*, 395–405. [CrossRef]

9. Ashton, G.; Bernstein, N.; Buchner, J.; Chen, X.; Csányi, G.; Fowlie, A.; Feroz, F.; Griffiths, M.; Handley, W.; Habeck, M.; et al. Nested sampling for physical scientists. *Nat. Rev. Methods Prim.* **2022**, *2*, 39. [CrossRef]

10. Pártay, L.B.; Bartók, A.P.; Csányi, G. Efficient Sampling of Atomic Configurational Spaces. *J. Phys. Chem.* **2010**, *114*, 10502–10512. [CrossRef]

11. Trassinelli, M. The Nested_fit Data Analysis Program. *Proceedings* **2019**, *33*, 14. [CrossRef]

12. Trassinelli, M. Bayesian data analysis tools for atomic physics. *Nucl. Instruments Methods Phys. Res. Sect. Beam Interact. Mater. Atoms* **2017**, *408*, 301–312. [CrossRef]

13. Trassinelli, M.; Ciccadicola, P. Mean Shift Cluster Recognition Method Implementation in the Nested Sampling Algorithm. *Entropy* **2020**, *22*, 185. [CrossRef] [PubMed]

14. Maillard, L.; Finocchi, F.; Trassinelli, M. Assessing Search and Unsupervised Clustering Algorithms in Nested Sampling. *Entropy* **2023**, *25*, 347. [CrossRef]

15. Shi, W.; Johnson, J. Histogram reweighting and finite-size scaling study of the Lennard-Jones fluids. *Fluid Phase Equilibria* **2001**, *187-188*, 171–191. [CrossRef]

16. Pártay, L.B.; Ortner, C.; Bartók, A.P.; Pickard, C.J.; Csányi, G. Polytypism in the ground state structure of the Lennard-Jonesium. *Phys. Chem. Chem. Phys.* **2017**, *19*, 19369–19376. [CrossRef]

17. Pártay, L.B.; Csányi, G.; Bernstein, N. Nested sampling for materials. *Eur. Phys. J.* **2021**, *94*, 159. [CrossRef]

18. Neal, R.M. Slice sampling. *Ann. Stat.* **2003**, *31*. 705–767. [CrossRef]

19. Handley, W.J.; Hobson, M.P.; Lasenby, A.N. PolyChord: Nested sampling for cosmology. *Mon. Not. R. Astron. Soc. Lett.* **2015**, *450*, L61–L65. [CrossRef]

20. Handley, W.J.; Hobson, M.P.; Lasenby, A.N. polychord: Next-generation nested sampling. *Mon. Not. R. Astron. Soc.* **2015**, *453*, 4385–4399. [CrossRef]

21. Henderson, R.W.; Goggans, P.M. Parallelized nested sampling. *Aip Conf. Proc.* **2014**, 1636, 100–105. [CrossRef]

22. Hüller, A. First order phase transitions in the canonical and the microcanonical ensemble. *Z. Für Phys. Condens. Matter* **1994**, *93*, 401–405. [CrossRef]

23. Medved', I.; Trník, A.; Vozár, L. Modeling of heat capacity peaks and enthalpy jumps of phase-change materials used for thermal energy storage. *Int. J. Heat Mass Transf.* **2017**, *107*, 123–132. [CrossRef]
24. Baldock, R.J.N.; Bernstein, N.; Salerno, K.M.; Pártay, L.B.; Csányi, G. Constant-pressure nested sampling with atomistic dynamics. *Phys. Rev. E* **2017**, *96*, 043311. [CrossRef]

25. Baldock, R.J.N.; Pártay, L.B.; Bartók, A.P.; Payne, M.C.; Csányi, G. Determining pressure-temperature phase diagrams of materials. *Phys. Rev. B* **2016**, *93*, 174108. [CrossRef]

26. Kokalj, A. XCrySDen—a new program for displaying crystalline structures and electron densities. *J. Mol. Graph. Model.* **1999**, *17*, 176–179. [CrossRef]

27. Mandelshtam, V.A.; Frantsuzov, P.A. Multiple structural transformations in Lennard-Jones clusters: Generic versus size-specific behavior. *J. Chem. Phys.* **2006**, *124*, 204511. [CrossRef] [PubMed]

28. Frantsuzov, P.A.; Mandelshtam, V.A. Size-temperature phase diagram for small Lennard-Jones clusters. *Phys. Rev. E* **2005**, *72*, 037102. [CrossRef]

29. Armour-Brown, C.; Borntraeger, C.; Fitzhardinge, J.; Floyd, P.; Hughes, T.; Jovanovic, P.; Jevtic, D.; Krohm, F.; Johnson, M.; Love, C.; et al. Valgrind. Available online: https://valgrind.org/ (accessed on 17 June 2025).

30. Nethercote, N.; Seward, J. Valgrind: A framework for heavyweight dynamic binary instrumentation. *Acm Sigplan Not.* **2007**, *42*, 89–100. [CrossRef]

31. Fonseca, J. gprof2dot. Available online: https://github.com/jrfonseca/gprof2dot (accessed on 17 June 2025).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.