Ab initio metadynamics determination of temperature-dependent free-energy landscape in ultrasmall silver clusters

Ab initio metadynamics enables extracting free-energy landscapes having the accuracy of first principles electronic structure methods. We introduce an interface between the PLUMED code that computes free-energy landscapes and enhanced-sampling algorithms and the ASE module, which includes several ab initio electronic structure codes. The interface is validated with a Lennard-Jones cluster free-energy landscape calculation by averaging multiple short metadynamics trajectories. We use this interface and analysis to estimate the free-energy landscape of Ag5 and Ag6 clusters at 10, 100 and 300 K with the radius of gyration and coordination number as collective variables, finding at most tens of meV in error. Relative free-energy differences between the planar and non-planar isomers of both clusters decrease with temperature, in agreement with previously proposed stabilization of non-planar isomers. Interestingly, we find that Ag6 is the smallest silver cluster where entropic effects at room temperature boost the non planar isomer probability to a competing state. The new ASE-PLUMED interface enables simulating nanosystem electronic properties at more realistic temperature-dependent conditions.


I. INTRODUCTION
Metadynamics (MTD) is a free-energy estimation method that enables exploring the conformational space of a system at a given temperature. It relies on the theoretical relation between the free energy of the system and a bias potential that drives the system to cross barriers and explore new conformations 1,2 . In principle, such algorithms can be coupled to any energy-force level description of the system. However, most applications prefer classical to quantum methods.
Metadynamics with quantum methods has been used to simulate chemical and biochemical reactions in gas phase, solid and in solution using Car Parrinello 3 , Born-Oppenheimer Molecular dynamics 4 , QM/MM metadynamics 5 . Some applications, for example, are allyl cyanide to pirrole isomerization 6 , formation of silver-chloro complexes 7 and water splitting and H 2 evolution by Ru(II)-Pincer complexes 8 . To overcome the limitation of short trajectories characteristic of quantum methods, minimum activation barriers have been reported stopping the metadynamics trajectory once the first transition is achieved and averaging over a few resulting barriers 8,9 or continuing a single trajectory and stopping the dynamics after one recrossing has been achieved 10,11 .
The lack of a good estimation of the resulting errors (due to such short trajectories) is hindering a more extended use of this important free-energy estimation method.
Noble metal nanoclusters have attracted much attention due to their molecular-like properties and high luminescence with potential applications in catalysis, biosensing and bioimaging 12 . Silver nanoclusters both bare and ligand-stabilized have a particular ability to form diverse structural motifs and a rich variety of isomers 13 . Experimental and simulated absorption spectrum of ultrasmall bare silver clusters indicates the coexistence of several isomers even at low-temperature starting at N a =6 and a transition from planar to three-dimensional for its lower energy isomer at N a =7 atoms [14][15][16] . This transition to non-planar structures is then much faster than its gold equivalent, which is placed at N a =11 up to T=100K 17,18 .
How would isomerization of silver clusters depend on temperature? Could the 2D-3D transition depend on temperature and other experimental conditions? Such questions require an estimation of the free-energy landscape of small silver clusters, which is computationally challenging to compute with traditional unbiased first-principles descriptors. Some studies on gold clusters have started to address these questions with the combination of quantum methods and enhanced sampling methods. Metadynamics applied to Au 12 clusters 19 shows that at room temperature, there is an equiprobable mixture of isomers. Recent work 18 on gold clusters predicts that at T=300K with N a =8 atoms there are non-planar isomers with non-negligible probabilities competing with planar isomers.
Our goal is to determine accurately the temperature dependent free-energy landscape of small neutral clusters with the use of ab-initio metadynamics in a general purpose interface. These results are useful to address questions of isomerization and the influence of stabilizers like solvents or organic matter. Accurate free-energy landscapes can also be used as benchmark for classical force-field developments 20 .
This paper is organized as follows. In the first section, we present a new interface and its validation with a classical Lennard-Jones cluster free-energy estimation. In the second section, we show an application of the ab-initio metadynamics method for the determination of the Ag 5 and Ag 6 free-energy difference between its lowest planar and non-planar isomer, and the corresponding minimum activation barrier. We also include in this section a discussion on ab-initio metadynamics parameters determination. We finalize with the perspectives and conclusions of the work.

II. ASE-PLUMED INTERFACE AND VALIDATION
Atomic Simulation Environment (ASE) is an open-source code written in python, with an object-oriented structure that allows setting, manipulating and running atomistic simulations. A central Atoms object obtains the energy and forces from a calculator object that in turn can be stand-alone or in interface to many of the widely used quantum and classical atomistic simulation codes. Therefore, with the aim to perform biased simulations, we created an interface to the opensource PLUMED library 21-23 by developing a calculator called Plumed that adds the energies and forces of bias to the forces and energies of other ASE calculators, which can be selected by the user. This interface between ASE and PLUMED opens the possibility to implement enhancedsampling methods with first-principles accuracy. We note that with PLUMED, it is possible to carry out several enhanced sampling methods and post-processing analysis of trajectories.
We focus here on the implementation of well-tempered metadynamics (WT-MTD) 2 . Metadynamics is a method that adds an artificial history-dependent bias potential over a small set of collective variables (CVs) 24 , s. Typically, the bias is accumulated as the sum of Gaussians centered along the CV trajectory. This pushes the system to explore different configurations and enhanced the sampling. In particular, in WT-MTD the bias potential at time t is where W is the initial height of the bias, τ is the time between deposited Gaussians, β is the inverse of k B T, γ is a bias factor and σ i is the width of the Gaussians for the i-th CV in s. Note that the first exponential decreases the height of the deposited Gaussians where previous bias energy has been added. This reduction of the Gaussian height reduces the error and avoids exploration towards high free energy states that are thermodynamically irrelevant. The rate at which the magnitude of the new added bias decreases is regulated by the bias factor γ: the lower the bias factor, the faster the decrease. The last exponential is a product of Gaussians in the direction of CV i with width σ i centered at the CV value at time t . In this way, the system's dynamics is enhanced, enabling it to explore different conformations.
For a sufficiently long exploration of the conformational space, it is possible to extract the free-energy landscape over the CVs (F(s)) using the information from the bias potential 2 This is the main advantage of MTD-derived methods.
As an accuracy test for our new Plumed calculator, we used a tutorial from the PLUMED documentation 22 as a benchmark system. This consists of WT-MTD/Langevin simulations for a simple system formed by seven atoms with Lennard-Jones (LJ) interactions in a planar space.
The LJ cluster has several stable isomers, which can be distinguished with the CVs of second and third central moments of the distribution of the coordination numbers (labeled by SCM and TCM respectively). The nth central moment µ n of the N a -atoms cluster is defined as where X i is the coordination number of the i-th atom: with r i j the distance between atoms i and j, and d a reference parameter. We used LJ dimensionless reduced units. The parameters of the simulation are d = 1.5, k B T = 0.1, friction coefficient fixed equal to 1, initial bias height of 0.05, Gaussian's width of 0.1 (for both CVs), and a bias factor of 5.
For this system, we compared the free-energy obtained by PLUMED as a stand-alone code and the free-energy estimated when using our new Plumed calculator that adds a bias force to a LJ-force calculator in ASE. For both cases, we ran 121 independent trajectories of WT-MTD, starting from the same configuration and random initial velocities, of duration 10 6 steps. In After testing the interface, we move into the problem of studying the temperature-dependence of the free energy landscape of small silver clusters.

III. FREE-ENERGY LANDSCAPE OF SMALL SILVER CLUSTERS
In the following sections, we present the DFT and WT-MTD setup for Ag 5 and Ag 6 cluster simulations, the free-energy surfaces and results extracted at different temperatures.

A. DFT parameters
Small neutral silver clusters have planar low-energy isomers, and as the number of atoms increases, the 3D isomers get closer in energy to the lowest energy isomer until Ag 7 , where the trend changes and the lowest energy configuration is a 3D structure 15,16 .
In the case of Ag 5 cluster, experimental studies with Raman and optical photoabsorption spectroscopy agree that the lowest energy isomer of this system has a planar trapezoidal shape 25 . Isomers of Ag 5 have been studied by means of computational methods such as Hartree-Fock, coupledcluster CCSD(T) and Density Functional Theory (DFT) 15,16,26 . All computational methods predict a 3D bipyramidal isomer which is about 0.4-0.5 eV higher in energy than the lowest energy state. Various methods predict other planar isomers with energies that are more strongly methoddependent. There is a planar isomer (edge-capped square) predicted by PBE, N12 and TPSS exchange-corrrelation functionals to be between the trapezoid isomer and the three-dimensional bipyramidal isomer, but not found with CCSD(T) method. Moreover, a planar isomer denoted bow-tie is found by CCSD(T) method 16 at the same energy of the three-dimensional isomer but other DFT functionals, like PBE and N12, place bow-tie isomer energy 0.2-0.3 eV below the threedimensional bipyramidal isomer energy 15 . In Fig. 2 (top), the trapezoidal lowest isomer is isomer 1 and the 3D bipyramidal is represented as isomer 3. The edge-capped square isomer is isomer 2 in Fig. 2 (top). Bow-tie isomer is isomer 4 in Fig. 2 (top).
Likewise, Ag 6 cluster has been studied experimentally and computationally, using DFT 15 and CCSD(T) method 16 . All computational methods suggest a triangular-planar isomer as configuration of minimum energy, followed by a 3D pyramidal isomer with a difference of energy of 0.1-0. In the column TPSS of Table I, we report optimization of isomers with TPSS exchangecorrelation functional and finite difference basis. As it was reported earlier for gold clusters 28 , this functional gives a good accuracy and reproducing the order and energies of CCSD(T) calculations. We suggest its use in future simulations, although it was out of reach for the computational resources used in this work. Here, we used PBE exchange-correlation functional with the LCAO basis which gives the right 2D-3D ordering with very good efficiency.
In order to take in account variations due to entropic effects and temperature, we apply WT-MTD to these two systems using the CVs described below.

B. Collective variables for WT-MTD
CVs are functions of the coordinates of the atoms that help to extract biophysical properties, separate relevant metastable states and reduce the dimensionality. CVs project the multidimensional system onto a small set of relevant and, in most cases, interpretative degrees of freedom. CVs are used in MTD and other enhanced-sampling methods (such as umbrella sampling 29 ) to add bias to the system. To guarantee a convergence, the CVs have to be chosen carefully, such that they differentiate the isomers and the transition states.
We studied several choices of CVs for the Ag 5 and Ag 6 systems. In particular, we focused on the coordination number C and radius of gyration R, which were previously used for studying the conformations of Au 12 cluster 19 . The average coordination number is where X i is defined in Eq. 4 and the reference distance d was set to 2.8Å to include all first neighbors distances in the silver isomers. This CV measures the number of bonds in the system.
The radius of gyration is where r i is the position of atom i, r CM is the center of mass of the cluster and N a is the number of atoms of the cluster. This CV gives information about how disperse the system is with respect to the center of mass. C and R enable extracting information about the shape of the cluster and permit differentiating the free-energy minima found by DFT optimization, which are expected to be metastable states in the free-energy landscape.
To evaluate more precisely the suitability of this set of CVs, we performed unbiased of MD in ASE for 10000 steps with a 5 fs time-step using the Born-Oppenheimer Approximation. The electronic distribution was obtained with LCAO-pvalence basis in a cell of 16 Å, using the GPAW calculator 30 . The temperature was controlled with a Berendsen thermostat at 10 K with a τ t of 50 fs (that will be the same setup used in WT-MTD for a wider range of temperatures). Starting from the trapezoid state and bipyramidal state, we observe that the form of the basins in the space of these collective variables were tilted ellipsoids in the unbiased MD ( Fig. S1 and Fig. S2 in Supplementary). Therefore, there are regions of the space that are thermodynamically irrelevant.
To avoid enhancing the exploration toward these regions, we created a new set of CVs (CV1 and CV2) that are a rotation of C and R, over which we could easily apply a wall. The rotated CVs are defined as Using this CV setup for WT-MTD, we added walls using repulsive semi-harmonic potentials that act when CV1 is lower than 5 with harmonic constant 10 eV and when CV2 is greater than 3 with harmonic constant 50 eV for Ag 5 (dashed lines in Fig. 4). In the case of Ag 6 , the walls repel values of CV1 lower than 8 with harmonic constant 10 eV and values of CV2 greater than 3.3 with harmonic constant 50 eV. We note that all isomers of Ag 5 and Ag 6 appear discriminated in the space of the CV1 and CV2. Moreover, this combination of the coordination number and radius of gyration in the CV2 variable is a good collective variable that allows to represent the FE of Ag 6 along a 1D profile (as will be described below).

C. WT-MTD parameters
We used the unbiased MD trajectories to determine the optimal parameters for the WT-MTD simulation. By monitoring the CVs as a function of time, we can estimate the MTD Gaussian width, which should approximate the amplitude of the CV at each minimum (bars in Fig. 3). In other words, the Gaussian widths are on the same order as the variation of the CVs in the unbiased simulation. Therefore, we choose the values of σ CV 1 and σ CV 2 fixed to 0.3 and 0.03, respectively. For choosing the other WT-MTD parameters, we performed several simulations using classic metadynamics to obtain an idea of the barrier height between different states, and extracting an optimal setup. From this exploration step, we decided to fix the initial height to 0.3 eV for Ag 5 and 0.2 eV for Ag 6 . The bias factor was fixed in 500, 100 and 50 for the temperatures 10, 100 and 300 K, respectively, in such a way that enabled the system to jump from the deepest minimum, but the Gaussians decrease sufficiently fast to achieve convergence in the simulated steps.

D. Ag 5 FES from low to room temperature
Using the paramters and CVs described above, we performed WT-MTD on Ag 5 clusters for 11 independent replicas and 50000 steps, resulting in a total of 250 ps. We obtained a clear difference in the exploration of the configurations compared with the unbiased simulation. For WT-MTD, the CVs filled the metastable state smoothly and then many transitions between states were observed. As a consequence, the system explored a large range of values in comparison with the maximum range covered by the simulations without bias (Fig. 3). This demonstrates that the free-energy landscape is being filled by the bias potential and that the system is recrossing the relevant metastable states. Importantly, we remark that it is not enough to find one single transition because the free-energy reconstruction will be poor. Therefore, the simulations ran until the error (calculated using N trajectories) was in the order of tens of meV, that condition usually implies more than 4 transitions between minima. are shown as dots.
We estimated the free-energy surface in the space of CV1 and CV2 as the average of the FE for the 11 trajectories for three different temperatures: 10K, 100K and 300K (Fig. 4 (top)). For all temperatures, the free-energy landscape contain just two minima, although four minima are obtained from optimization. At these temperatures, states 2 and 4 (shown in Fig. 2), corresponding with the edge-capped square and bow-tie isomers, are just saddle points that belong to state 1 (the trapezoid isomer). Therefore, only states 1 and 3 (shown in Fig. 2) are representative configura- For low temperatures, we expect only a small variation between the computed zero temperature energy obtained via DFT optimization and the 10K free energy obtained with ab initio metadynamics because of the underlying assumptions of fixed Boltzmann statistics in metadynamics and Born-Oppenheimer electron-nucleus decoupling in DFT. We use therefore the 10K FES values as convergence check inside our trend study, which converge to the expected values. We note that quantum effects that are not valid under these assumptions will not be captured in our simulations.
The effect of increasing the temperature is given by a decrease in the minimum activation barrier and the free-energy difference between minima between the planar and non-planar isomers. These results are shown in Fig. 5, finding that the transition state decreases approximately 0.03 eV and the free energy difference decreases by 0.09 eV from 10K to 300K.
A more dramatic change is captured when calculating the relative population of the basins using the Boltzmann factor. We define a basin as the region where the free-energy surface is less than the value of the minimum activation barrier (Fig. 5). Then, we obtained the probability of each state by integrating the Boltzmann factor over the correspondent basin, i.e. the probability of state 1 is P 1 = 1 exp(−β F(s))ds where β = 1/(k B T ) and F(s) is the free-energy at s. Interestingly, the probability associated to all the non-planar isomers is negligible for all temperature ranges (even 300K), namely, the probability to find a planar configuration is 100% for Ag 5 .
Using statistical bootstrapping with 50 resamples, we explored how many independent simulations are required to extract an error by varying the number of samples in each resampling. This gives a notion of how the predicted result changes as a function of the number of simulated replicas. In Fig. 6, we show the mean value (dots) and the standard deviation (bars) of the difference in free-energy between the isomers 1 and 3 of Ag 5 . This result demonstrates the importance of running at least 4 replicas for obtaining a reliable free-energy difference estimate. We note that when using only one metadynamics simulation, the results can significantly change, even up to 0.15 eV, which is a large variation compared with the value of this observable. This demonstrates the importance of considering several replicas for relatively short simulations. We note that the exact convergence rate will depend also on the complexity of the particular system, the simulation length and WT-MTD setup. E. Ag 6 FES from low to room temperature We also studied the free-energy landscape of the Ag 6 cluster, running 4 independent trajectories with 136000 steps, resulting in a total of 680 ps. The cluster has three stable isomers according to the optimization analysis (Table I and Fig. 2). In Supplementary Fig. S3, we show the FES along both CVs at the simulated temperatures. We note that only the isomers 1 and 2 are stable states of the free-energy landscape. The incomplete hexagon isomer, the third Ag 6 isomer, appears as part of the basin of the isomer 1. From that figure, it is clear that states 1 and 2 are properly separated along CV2. Therefore, integrating CV1 enables a clear representation of a free-energy profile along CV2 (i.e. exp(−β F(CV 2)) = exp(−β F(CV 1,CV 2)) dCV 1.). The averaged profile over the 4 trajectories, and a shaded region representing the standard error are shown in Fig. 7.
For the three cases, the standard error remains lower than 0.04 eV.
From 10 K to room temperature, the free-energy difference decreases approximately by 0.2 eV, but the barrier with respect to the global minima decreases by only 0.02 eV. In terms of probability, however, the change is drastic. At room temperature the probability of non-planar isomer reaches 10% which shows that the system reaches a new equilibrium where planar and non-planar isomers are competing. This is a marked difference between Ag 5 and Ag 6 , and it is in accordance to an observed change in optical spectrum experiments 27 . for the isomers 1 and 2 in Fig. 2.

IV. CONCLUSION
We developed a new calculator for ASE called Plumed, which can be used for running simulations of enhanced-sampling methods by patching the open-source code ASE and the plugin PLUMED. This calculator was tested with a simple system of seven LJ atoms as benchmark. This ASE-PLUMED interface was used for studying Ag 5 and Ag 6 clusters at different temperatures.
We found crucial thermal effects over the Ag 6 system, which changes from a planar-dominated population at low temperatures to a state with mixture of planar and non planar isomers at room temperature. Because no changes in population are found in Ag 5 at the same temperature range, it follows then that Ag 6 is the smallest silver cluster with a 2D-3D isomer equilibrium at room temperature.
As an important result, we studied the convergence of the WT-MTD technique, showing that when considering trajectories with few recrossings using just one replica could generate unreliable results but averaging over several independent replicas decreases significantly the error. All-in-all, we foresee that the ASE-PLUMED interface as a general purpose tool for enhanced sampling simulations having the accuracy of first principles electronic structure methods will expand this method use in the material's science community.

V. DATA AVAILABILITY STATEMENTS
The data that supports the findings of this study are available within the article and its supplementary material. The code developed in this study is openly available in gitlab (https://gitlab.com/ase/ase).