{'fraction_non_alphanumeric': 0.024172591321402304, 'fraction_numerical': 0.007403775435155676, 'mean_word_length': 5.186229905975129, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite temperature properties of clusters by replica exchange metadynamics: the water nonamer

We introduce an approach for the accurate calculation of thermal properties of classical nanoclusters. Based on a recently developed enhanced sampling technique, replica exchange metadynamics, the method yields the true free energy of each relevant cluster structure, directly sampling its basin and measuring its occupancy in full equilibrium. All entropy sources, whether vibrational, rotational anharmonic and especially configurational -- the latter often forgotten in many cluster studies -- are automatically included. For the present demonstration we choose the water nonamer (H2O)9, an extremely simple cluster which nonetheless displays a sufficient complexity and interesting physics in its relevant structure spectrum. Within a standard TIP4P potential description of water, we find that the nonamer second relevant structure possesses a higher configurational entropy than the first, so that the two free energies surprisingly cross for increasing temperature.


Introduction
The physics of clusters formed by a small number of atoms or molecules is basic to many branches of science. 1 Clusters are the primary nanosystem where matter behaves differently from its usual bulk aggregation. In condensed matter physics it has long been of interest how collective features of macroscopic systems, such as ordering phenomena and phase transitions, emerge in a nutshell in small size clusters. 2,3 As is well known, it often proves harder to predict the properties of a small aggregate than those of an infinitely large system, where one can take advantage of regularity. In computational physics, identifying and characterizing cluster structures provide a prototypically difficult optimization problem. 4 The potential energy landscape as a function of the atomic or molecular coordinates is generally complex, with many local minima separated by large barriers even in relatively small clusters, where the simple identification of the lowest energy structure may be a serious computational challenge. State-of-the-art optimization techniques are for example sometimes benchmarked on their capability to find the global energy minimum of Lennard-Jones clusters. [5][6][7] Simulation studies have been devised, using enhanced sampling methods, for example aggregation-volume-bias Monte Carlo and umbrella sampling, to obtain the free energy and density of state of clusters or crystalline nuclei in liquids. [8][9][10][11] The recent boost of nanoscience increasingly demands a better and better ability to characterize, calculate, and understand theoretically the peculiar properties of matter at the nanoscale, a subject of broad technological, practical as well as conceptual implications.
The traditional manner to approach cluster simulation is quite simple. As a first step one attempts to identify all the locally stable low energy structures, usually by some global optimization technique, such as simulated annealing, genetic algorithms, etc . The computational cost of these techniques is largely dominated by repeated potential energy evaluation, and one must carefully choose a level of description of the mechanics of atomic coordinates (quantum or semiempirical or classical) and of the electronic degrees of freedom allowing at the same time an exhaustive search and a sufficiently accurate estimate of the energy. Sometimes a preliminary search with a less accurate empirical intermolecular potential can be followed by a subsequent refinement leading to a more accurate energy estimate and thus to the identification of the relevant structures with a higher level description. If one is interested in the properties of the system at a very low temperature, what matters is the global energy minimum. When however temperature grows so that k B T is comparable to the energy gap ∆ 01 above the global minimum, higher energy states become relevant and must be identified. It is traditional to consider a relatively small number of "relevant states", analogous to inherent states in glass physics, crudely consisting of locally stable atomic configurations or cluster geometries, that identify local energy minima in the space of atomic coordinates.
A finite set of relevant states can be obtained for example by (real or conceptual) annealing of the system from a much larger number and variety of finite temperature instantaneous molecular configurations, cooled down to reach the closest local energy minimum at T = 0. At finite temperature, each relevant state s represents the center of a basin in configuration space. The logarithm of the volume of the relevant state basin volume in configuration space represents its entropy S s .
When the thermal occupancy of each relevant structure s is explicitly calculable, this is equivalent to specifying its free energy F s = E s − T S s . A crucial question that needs to be addressed is, what contributes to the entropy of the relevant state, and how it can be properly calculated. The simplest and commonest approach is to consider free rotations and harmonic vibrations around the T = 0 local energy minimum of the relevant state, and calculate the vibrational and rotational entropy involved, in the rigid-rotor-harmonic approximations. 8,12,13 This is however not generally adequate, and in fact the rot-vibrational approximation may yield free energies of relevant structures in serious error. The full entropy of a state, measuring the size of its stability basin should include not just vibrations and rotations, but all possible coordinate transformations leading to essentially equivalent configurations. Proper determination of configurational entropy is a factor which critically reflects the abundance of a relevant state in the cluster's thermal equilibrium. Although these concepts are natural, and well rooted in classic cluster literature, 2,3 they are otherwise not always properly implemented in everyday practice where configurational entropy, difficult to define and evaluate, is often overlooked. A practical and accurate calculation approach including automatically all sources of entropy, especially for the more complex molecular clusters, is highly desirable.
In this paper we introduce and demonstrate a procedure leading to a computationally convenient and conceptually correct extraction of the free energy of the relevant states of a cluster at finite temperature, including all non-electronic entropy sources contributions. In the simple case we use for our demonstration, a small water cluster, the entropy is found to differ drastically between one relevant state and another. This working case also shows explicitly that the neglect of configurational contributions leads to large errors, totally unacceptable practically no less than conceptually.
Our method simultaneously addresses both steps normally followed in traditional cluster study techniques: a) the search for the relevant structures, and b) the study of their finite-temperature properties. The approach, based on a recently developed enhanced sampling technique, replica exchange metadynamics, 14 allows the computation of the true free energy of each relevant structure, directly sampling its basin and measuring its occupancy (overall abundance) in full thermal equilibrium, thus including automatically all entropy sources.
For our demonstration purposes we chose the water nonamer (H 2 O) 9 , a very simple dielectric cluster which nonetheless displays a sufficient complexity and interesting physics in its relevant structure spectrum. Although not directly addressed in this study, well characterized spectroscopic properties of this particular cluster are available. 15 We studied the cluster properties covering a temperature range of potential interest for atmospheric science (100-250 K), and described the forces between water molecules with a simple but realistic and widely tested empirical force scheme, the so-called TIP4P potential. 16 Since the electronic shells of this system are closed and widely gapped, the use of such effective interatomic forces is not unreasonable, although clearly not definitive. The cluster thermal evolution is simulated by replica exchange metadynamics (RE-META), 14 an approach which requires the simultaneous running of several molecular dynamics simulations at different temperatures. Each simulation ("replica") is biased by an additional history-dependent "metadynamics potential". 17 This artificial bias energy term acts on an appropriately chosen small set of collective variables (CV), global or local order parameter-like quantities representing some physically motivated function of molecular coordinates. At fixed time intervals, the replicas are forced to attempt an exchange of their coordinate configurations.
The exchange move is accepted or rejected according to the usual Metropolis criterion, 18,19 including the change of the history dependent CV-based bias potential between the two replicas.
The configurations visited in the RE-META simulation correspond to visiting the free energy landscape of the cluster -as a function of the CV variables -in parallel at all temperatures. Combining replica exchange and metadynamics leads to an overall kinetics which, even at low temperatures, rapidly becomes diffusive in collective variable space, resulting in an extremely efficient configuration sampling. 14 In RE-META, replica exchange is improved because metadynamics explores high free energy configurations within each replica, ultimately leading to faster decorrelation. Metadynamics is in turn strengthened, since replica exchange improves the sampling of the degrees of freedom that are not explicitly biased -thus mitigating one major problem of metadynamics. 14 The water nonamer example confirms that the free energies of the lowest relevant states are indeed heavily affected by configurational entropies. In particular, the entropy of the first excited state is higher than that of the ground state, which causes for increasing temperature an outright crossing of the two free energies and populations, a crossing which does not occur in the vibrational and rotational approximation. K. Since TIP4P bulk ice has a melting temperature of 232 K at 1 bar, 20 the higher temperatures correspond to a fast diffusing liquid-like regime, the lower ones to a slow diffusing solid-like one.

The (H 2 O) 9 cluster: RE-META calculations
Our level of description is entirely classical, and no quantum effects are included, either nuclear or electonic. Because therefore our treatment ignores all quantum freezing and zero point motion effects, the simulation results lose their strict validity at very low and zero temperatures, a limit where they should not be extrapolated. The molecular dynamics simulation was performed in the NVT ensemble and a cell volume of 18850 Å 3 , a volume where the probability of observing an isolated water vapour molecule is extremely small. A Nóse-Hoover chain thermostat was used to control the temperature, the integration time step was set to 0.2 fs, and configurations were stored for analysis every ps. Each replica was evolved for a total simulation time of 50 ns and independently biased by a metadynamics potential acting on a single collective variable S counting approximately the number of hydrogen bonds inside the cluster: where the sum runs over proton-oxygen pair bonds belonging to different molecules and r 0 = 2 Å.
The metadynamics potential has the form 17 where s (t) = S (r (t)), δ s = 0.15 and w = 0.004 meV. As customary in metadynamics, the collective variable defined in ?? is chosen on purely physical grounds, and is neither unique nor systematically definable in a general case. However, while this arbitrariness could be seen as endangering the reliability of metadynamics, the replica exchange strategy does, as shown in Ref., 14 effectively cure the problem. In our calculation the exchange moves between configurations taken at different temperatures were attempted every 2 ps, accepting the moves with a probability min(1, exp (−∆)) with where i, j is the index of the replica, V is the sum of the potential energy and the bias, r are the coordinates and β is the inverse temperature, which is different in each replica. 14 We found that after a typical equilibration time τ F ∼ 4 ns the history-dependent potential started growing evenly in CV space -meaning independent of S -in all replicas, indicating that -V G (S,t) will eventually, for large t, provide a good estimate of the free energy as a function of S. 21 After the equilibration time τ F we could therefore stop changing the history-dependent potential V G and start collecting the statistics of occurrence of properly defined relevant states, with the scope of computing their respective probability, and from that their free energies and entropies in thermal equilibrium.
In order to analyze the results of the procedure it is necessary to assign during the MD simulation each molecular coordinate configuration, to the closest reference structure -its relevant state.
If all the possible local minima of the potential energy surface were known, one could classify the configurations according to their distance from the reference structures, measured, for instance, by the appropriate root mean square deviation of coordinates. This however would require a preliminary search of all the minima. More importantly, in order to assign a configuration to the correct relevant state, one should consider all permutations between identical atoms, a procedure hugely cumbersome and computationally expensive. We apply here a more conservative procedure that uses ideas from the basin hopping technique. 6 Starting from all the coordinate configurations generated by MD, we carried out a static potential energy minimization by conjugate gradients evolving atomic coordinates until the force on each atom was less than 0.0001 eV/Å. If a reference structure with the same potential energy is already available, the instantaneous MD configuration is assigned to that structure. Otherwise, a new reference structure is introduced in the pool. Here two quenched structures are assumed to be the same if their energy difference is smaller than 0.0001 eV. This procedure is still computationally intensive, but allows an unambiguous assignment of all MD frames and is insensitive to permutation of the atoms. We now come to our main point. Besides providing a tool for the quick exploration of cluster structures, RE-META crucially yields an accurate estimate of the probability to observe each relevant structure S as a function of temperature. Following Ref, 23 we first compute the average bias potential where τ F is the equilibration time of V G , 4 ns in this work and τ sim is the total simulation time.
Denoting by r i all the configurations generated by RE-META at temperature T , we estimate the equilibrium probability to observe relevant structure α as where the notation i ∈ α means that the sum runs over the configurations that are assigned to structure α as explained above. In order to assess the reliability of this approach, we also computed the probability p α in normal replica exchange where the probability is simply proportional to the number of times the structure α is observed. As shown in [figure] [3][]3-a,b, the probabilities estimated in the two manners are consistent. The small remaining discrepancy between RE and RE-META is due to statistical uncertainties. We remark that the errors of methods based on statistical sampling are inevitably large at low T, as the relevant transitions are observed relatively rarely even if a powerful enhanced sampling approach is used.  Figure 3: (a) occupancy probability for the lowest energy relevant state ("0") versus temperature, with RE-META, normal RE, rigid-rotor-harmonic approximation ("Rot-Vib") and rigid-rotorharmonic approximation with an additional the ln 2 factor in the entropy of state 1 -state 8 ("Standard"). (b) occupancy probability for the second lowest energy relevant state ("1"). The statistical error bars in RE and RE-META represent the standard deviations estimated by separately comparing averages in three successive portions of the simulations. Note the substantial agreement between RE and RE-META and the large deviation of the rot-vib approximation, showing its clear inadequacy. Note also that the occupancies of states computed with RE and RE-META are fully consistent with those computed using the standard approach at low temperature, and how the probability of state 1 overtakes that of state 0 around 150 K, as a result of its larger configuration entropy.

Comparison between RE-META and rigid-rotor-harmonic vibration approximations
We calculate and compare in [figure] [3][]3-a,b the exact occupancy probabilities of the two lowest relevant structures with the approximate probabilities, now estimated using the rigid-rotorharmonic vibration approximation to the cluster vibrational spectrum of each state: where the rigid-rotor-harmonic free energies F α are calculated as where E α is the energy and where I α i are the principal moments of inertia, σ α is the rotational symmetry number and ω α m is the vibrational frequency of mode m obtained by diagonalizing the Hessian matrix of the potential in structure α at T =0. We finally calculated the entropy difference between states 0 and 1 where ∆F and ∆E are the free energy and the internal energy differences. The state 0 has a mirror symmetry, while state 1 does not. The broken symmetry implies that state 1 is two-fold degenerate. Thus, the relative configurational entropy of state 1 is roughly ln 2, which, at low temperature, approximately agrees with the relatively large occupancy of state 1.
Adding ln 2 to the rot-vib entropy for states (but not for state 0) leads to a positive "Standard"  Figure 4: Entropy difference between state 1 and state 0 as a function of T. Note the initially larger configurational entropy of state 1, further growing with temperature, due to anharmonicity. Both effects are missed by the rot-vib approximation, whereas the addition of ln(2) (the "Standard" approximation) provides a reliable result only well below 100 K.

Discussion and conclusions
The cluster description obtained with RE-META can naturally be used to predict measurable quantities. As an example, we computed the average dipole of the water nonamer as a function of T. In Beyond the deliberately simple case study (H 2 O) 9 presented here, RE-META method is applicable to larger water clusters, and of course to different systems, once the collective variables are properly selected. The workload and complexity grows considerably with size -not surprisingly in view of the exponential growth of phase space -but the effectiveness is not impaired. As an example we considered the larger water cluster (H 2 O) 14 , hundreds of times more complex than the nonamer with a 10 ns simulation. [figure] [6][]6 shows the number of independent structures explored as a function of simulation temperatures for RE and RE-META, underlining the strong  Figure 5: Thermal averaged dipole moment of a water nonamer, computed directly from the RE-META and RE trajectories and using the probabilities ( ??) estimated within the "Rot-Vib" and "Standard" approximation.
superiority of the latter, in particular its ability to unearth relevant structures that would otherwise be simply ignored. In summary, we have shown how the full configurational complexity of a small cluster can be directly, realistically and economically accessed by a suitable sampling, combining replica exchange with metadynamics techniques (RE-META). This proves extremely effective in searching the cluster's relevant states, and in measuring their thermal population as a function of T. In the water nonamer, these populations are shown to differ largely from those routinely predicted by rot-vib approximation. The crucial factor which influences the populations is the configurational entropy of the relevant states, an entropy which is generally inaccessible but which in RE-META is naturally and directly measured and understood. The method holds a good promise of applicability in a great variety of relevant systems, provided a set of good, physically motivated collective variables can be identified.