Chemical Potential Calculations In Dense Liquids Using Metadynamics
The calculation of chemical potential has traditionally been a challenge in atomistic simulations. One of the most used approaches is Widom’s insertion method in which the chemical potential is calculated by periodically attempting to insert an extra particle in the system. In dense systems this method fails since the insertion probability is very low. In this paper we show that in a homogeneous fluid the insertion probability can be increased using metadynamics. We test our method on a supercooled high density binary Lennard-Jones fluid. We find that we can obtain efficiently converged results even when Widom’s method fails.
Eur. Phys. J. Spec. Top. DOI: 10.1140/epjst/e2016-60094-x
Chemical potential is an important thermodynamic quantity that regulates many phenomena, such as phase transitions, reaction equilibria and electro-chemistry JobEJP2006 (); BaierleinAJP2001 (). The ability to calculate it via a reliable and fast technique is therefore of great value. Several numerical methods have been proposed, ranging from Widom test-particle insertion WidomJCP1963 () to Grand-Canonical Monte-Carlo approaches AdamsMP1974 (); AdamsMP1975 (). Extensive reviews on the efficiency and reliability of these and other techniques can be found in the literature KofkeMP1997 (); ShirtsJCP2005 (); DalyCPC2012 (). However, despite these efforts, chemical potential calculations still remain a difficult challenge for many real systems, such as dense fluids and solids.
At constant volume and temperature the chemical potential is equal to the Helmoltz Free-Energy difference associated to the addition of an extra molecule. The most common approach to calculate this quantity, namely the Widom insertion method, consists in sampling the insertion energy of a test-particle. In this method the largest contribution comes from those configurations in which the insertion of an extra molecule is energetically favorable. This occurs if a cavity sufficiently large to accommodate the extra particle is present. In a low density system many such cavities are present, and Widom method is efficient. On the contrary, in dense systems the formation of a particle-sized cavity is unlikely, and the calculation becomes inefficient. Thus elaborate methods, based on Voronoi tessellation, have been suggested as a way to look for the appropriate cavities in a dense system BieshaarMS1995 (). Further suggestions have been among others the introduction of energy-biased sampling of the insertion space DelgadoJCP2005 (), the use of external potentials to generate controlled inhomogeneities within the system PowlesJCP1994 (); MooreJCP2011 () or the application of an adaptive resolution simulation scheme AgarwalJCP2014 ().
In the present paper we shall focus on the calculation of the chemical potential for an homogeneous liquid, and extend the range of densities at which the particle insertion method can be applied by using Well-Tempered (WT) metadynamics BarducciPRL2008 (). Here the role of metadynamics is to enhance fluctuations in a controlled manner, promoting those rare configurations in which particle insertion is favored. This leads to an efficient and reliable calculation method, even in those cases in which Widom’s method fails.
The paper is organized as follows: in Sec. II metadynamics main features are recalled, while its application to chemical potential calculation is described in Sec. III. The implementation of our technique for the selected test case, namely the binary Lennard-Jones (LJ) fluid, is addressed in Sec. IV, while the results of our calculations are reported in Sec. V. Finally the conclusions and perspectives of present work are reported in Sec. VI.
Before discussing the chemical potential calculation, we recall some of the features of WT metadynamics. For a more detailed description we refer to the by now vast literature on the subject (see e.g. ValssonRev2015, ; BarducciWIR2011, ).
Metadynamics introduces a bias potential to enhance the sampling of rare but important fluctuations. This bias potential is taken to be a function of the so called Collective Variables (CV), a set of order parameters that depend on the atomic coordinates . The CVs are built to describe the degrees of freedom of the system whose sampling we want to enhance. Let us restrict for simplicity to the case of a single collective variable . The equilibrium distribution of is given by:
where is the Boltzmann distribution defined by potential energy function . is the inverse temperature and is the partition function of the system. Introducing the free-energy surface , where is a constant that plays an inessential role, Eq. (1) can be rewritten as
If a bias potential is introduced, a different distribution is sampled:
where is the Boltzmann distribution of the biased ensemble. is the corresponding partition function.
In WT metadynamics the bias is periodically updated by adding a small repulsive potential, that usually takes the form of a Gaussian of height and width :
is centered on , that is the value of the CV at the moment of the update. During the simulation the height of the deposited Gaussians is modified in a way that can be described by the following iterative equation:
where indicates the iteration and is the so called biasing factor, a key parameter in WT metadynamics. After the -th iteration the system evolves under the combined action of and of .
As a consequence, the biased distribution is related to the unbiased one by the relation:
This result is central to the present paper. In fact Eq. (7) shows that narrow peaks in can be made broader in by increasing . Thus rare but important fluctuations, that are hidden in the tails of are more likely to occur in the biased ensemble.
Furthermore, the theory of metadynamics allows to calculate the unbiased ensemble average of any position dependent variable via the reweighting TiwaryJPCB2015 (); BonomiJCC2009 () procedure described by:
where the suffix indicates that the average is performed under the action of the time varying metadynamics bias potential . The position independent constant is given by TiwaryJPCB2015 ():
Iii Choice of the Collective Variable
We now define the proper CV that will be used to calculate the chemical potential via metadynamics. Let us begin by considering a system of particles with coordinate vector . As pointed out by Widom WidomJCP1963 (), the excess chemical potential, that is the chemical potential from which the free particle contribution has been subtracted, can be written as:
where is the system volume and
where is the potential energy of the system with an extra particle in . The average in Eq. (10) is performed over the Boltzmann distribution defined by the particle interaction potential . Widom suggested to estimate by evolving the system and periodically computing at a randomly generated set of positions . Also in other methods, first the coordinate vector is sampled and later insertions are attempted BieshaarMS1995 (); BennettJCP1976 ().
Here we take an alternative route and first we fix a number of insertion points and later we let adjust around them so that the sampling of is enhanced. In order to see how this is done we restrict to the case of an homogeneous fluid, in which the average is independent. It follows that one can choose a set of ’s and rewrite the chemical potential as:
This form of suggests the introduction of the following CV:
In principle even the choice , namely that of using a single insertion point , leads to the correct result, however it requires much longer sampling times. In this paper the insertion points will be arranged on a regular grid, but other choices are possible.
A consequence of this CV definition is that the chemical potential can be rewritten as the expectation value of a function of :
or, in terms of the free energy associated to , as:
This choice of CV has a practical drawback. In fact in most models the atoms interact at short distances via a highly repulsive and rapidly varying potential. Thus, if any of the insertion points comes close to the atomic positions, varies very rapidly and can even diverge. This can be remedied by modifying the insertion potential at short distances using a regularized potential in the definition of a new and more smoothly behaved CV:
The idea is that should be large but finite when is close to one of the atomic positions, so that the derivatives are of manageable size. One can then use to bias the metadynamics simulation and calculate using the reweighting procedure in Eq. (8) TiwaryJPCB2015 (). However, in the test case that will be discussed below, we define so that the change of CV has no practical effect on the calculation of , and the reweighting is not necessary (see also the Supporting Information (SI)).
Iv Test Case
In order to test the method we shall consider a binary Lennard-Jones (LJ) fluid. In particular we consider the system studied by Kob and Andersen KobPRL1994 (); KobPRE1995 () for its glass forming properties, that it is known to remain fluid at low temperature and high density. Thus it provides a challenging test for our method. The system is composed by two atomic species, and , of equal atomic mass. The interaction potential is:
where is the atom-atom distance and . In such a system one defines an excess chemical potential for each species, and generalizes Eq. (12):
where and are the coordinates of the atoms and indicates the ensemble of the two combined sets of coordinates. is the number of specie particles.
Each species will have its own , therefore Eq. (13) generalizes to:
Analogously, Eq. (14) becomes:
and Eq. (15) becomes:
As mentioned before, in order to avoid divergences and the appearance of a large bias, we shall use the regularized version of the CV, defining the following modified potential:
where is the distance below which and differ.
V Results and Discussion
Kob and Andersen model KobPRE1995 () is a mixture of particles and particles. In units of and the other potential parameters are , , and . In our simulations the potentials are cut and shifted at , while the distance parameter of the regularized potential is set at . The simulated system is composed of atoms in a cubic box of length and Periodic Boundary Conditions (PBC) are applied. The temperature is kept at by using the stochastic velocity rescaling thermostat BussiJCP2007 (). At these thermodynamic conditions the system is in the liquid phase.
The MD simulations were carried out using the GROMACS HessJCTC2008 () software, equipped with a private version of the PLUMED2 plug-in TribelloCPC2014 (). The integration time-step was chosen to be (which corresponds to for the LJ parameters of Argon: , and Å). The Widom’s calculations were performed using the appropriate GROMACS subroutines.
One advantage of using this particular binary system is that the majority of atoms () have a larger radius, leaving sufficient interstices to place the smaller atoms. This makes the calculation of relatively easy and accessible to Widom’s method, providing a testing ground for our technique. On the other hand the calculation of is challenging and its successful completion demonstrates the usefulness of our approach. We performed two separate calculations, for and . The Gaussian initial height was chosen to be of and the width is equal to . New Gaussians were deposited every and the parameter was .
We experimented with different values of number of insertion points (see Eq. (20)). They all lead to the same results within statistical error. However, the speed of convergence of was different, with the limiting case being particularly slow (see the SI). Here we present the results for a regular mesh of points. While was not carefully optimized, this choice seemed to yield a good compromise between speed of convergence and computational cost. We note that the calculations can be easily parallelized and it pays to distribute wisely the load among the concurrent processors, thus the determination of the optimal depends also on the computational platform.
As discussed earlier we expect the chemical potential calculation to be easier for the smaller atoms () and more challenging for the large component (). This is brought out in Fig. 1, where the resulting from metadynamics simulations is compared to the calculated in an unbiased MD run. In the case (Fig. 1a), even without applying the bias potential, the range of spontaneously occurring fluctuations, indicated in figure by the domain of the unbiased , is large enough to overlap with the negative region. This region yields the more substantial contribution to the numerator of Eq. (22), and thus to the value (see SI for further details). Thus can be reliably computed with unbiased sampling. Note however that even in this favorable case the spontaneous fluctuations do not cover entirely the tail of the distribution. As noted in Ref. LuJCP2001a () this leads to a systematic error. In contrast, no such problem is present in the metadynamics calculation, where the fluctuations of cover thoroughly all the important region.
In the case of the larger atoms (Fig. 1b) the situation is much more dramatic and the unbiased results show little or no overlap between the fluctuations of and the region that contributes the most to . In contrast, our metadynamics calculations explore in detail the important region of . This is reflected by the convergence of and calculations, shown in Fig. 2, where we compare metadynamics results (obtained with Eq. (8)) to standard Widom calculations. In the case (Fig. 2a) the agreement between metadynamics and standard Widom method is within the statistical error. In the case however (Fig. 2b) Widom’s method faces severe difficulties and convergence is jumpy and very slow, even with a very large number of insertion points sampled. On the contrary our method converges to an accuracy of using only insertion points. This underlines the necessity of metadynamics to visit the extreme tails of the distribution and calculate accurately.
In conclusion we have shown that this new approach can be very competitive with previous ones. Three features should be noted: the use of the homogeneity of the system, the change of paradigm from an atom position-centric to an insertion-centric point of view and the use of metadynamics to enhance fluctuations. The calculation is no more difficult than a standard WT metadynamics one, and can be performed using any of the codes with which the PLUMED2 plug-in can be integrated. The recently developed variationally enhanced sampling ValssonPRL2014 (), in its WT version ValssonJCTC2015 (), could be also profitably applied to perform these calculations and lead to further improvement in the efficiency of chemical potential evaluations. Although here our method has been developed for homogeneous systems, its extension to non-homogeneous situations is possible and it is currently under study. A relatively straightforward application of this method is the calculation of the solvation energy of molecules in liquids, another quantity of great practical interest. Moreover, if properly adapted, these ideas could be extended to free energy perturbation calculations ZwanzigJCP1954 (); Allenbook1989 (); KollmanCR1993 (); Frenkelbook2001 ().
Acknowledgements.The authors are pleased to acknowledge M. Nava, B. Smit and O. Valsson for the useful discussions. The computational resources were provided by the Swiss Center for Scientific Computing, and the Brutus Cluster at ETH Zurich. The authors acknowledge also the VARMET European Union Grant ERC-2014-ADG-670227 and the National Centres of Competence in Research “Materials Revolution: Computational Design and Discovery of Novel Materials” project for funding.
- (1) G. Job, F. Herrmann, Eur. J.Phys. 27(2), (2006) 353
- (2) R. Baierlein, Am. J. Phys. 69(4), (2001) 423
- (3) B. Widom, J. Chem. Phys. 39(11), (1963) 2808
- (4) D. J. Adams, Mol. Phys. 28(5), (1974) 1241
- (5) D. J. Adams, Mol. Phys. 29(1), (1975) 307
- (6) D. A. Kofke, P. T. Cummings, Mol. Phys. 92(6), (1997) 973
- (7) M. R. Shirts, V. S. Pande, J. Chem. Phys. 122(14), (2005) 144107
- (8) K. B. Daly, J. B. Benziger, P. G. Debenedetti, A. Z. Panagiotopoulos, Comput. Phys. Commun. 183(10), (2012) 2054
- (9) R. Bieshaar, A. Geiger, N. N. Medvedev, Mol. Simulat. 15(3), (1995) 189
- (10) R. Delgado-Buscalioni, G. De Fabritiis, P. V. Coveney, J. Chem. Phys. 123(5), (2005) 054105
- (11) J. G. Powles, B. Holtz, W. A. B. Evans, J. Chem. Phys. 101(9), (1994) 7804
- (12) S. G. Moore, D. R. Wheeler, J. Chem. Phys. 134(11), (2011) 114514
- (13) A. Agarwal, H. Wang, C. Schütte, L. Delle Site, J. Chem. Phys. 141(3), (2014) 034102
- (14) A. Barducci, G. Bussi, M. Parrinello, Phys. Rev. Lett. 100(2), (2008) 020603
- (15) O. Valsson, P. Tiwary, and M. Parrinello, Annu. Rev. Phys. Chem. 67(1), (2016)
- (16) A. Barducci, M. Bonomi, M. Parrinello, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1(5), (2011) 826
- (17) J. F. Dama, M. Parrinello, G. A. Voth, Phys. Rev. Lett. 112(24), (2014) 240602
- (18) P. Tiwary, M. Parrinello, J. Phys. Chem. B 119(3), (2015) 736
- (19) M. Bonomi, A. Barducci, M. Parrinello, J. Comput. Chem. 30(11), (2009) 1615
- (20) C. H. Bennett, J. Comput. Phys. 22(2), (1976) 245
- (21) W. Kob, H. C. Andersen, Phys. Rev. Lett. 73(10), (1994) 1376
- (22) W. Kob, H. C. Andersen, Phys. Rev. E 51(5), (1995) 4626
- (23) G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126(1), (2007) 014101
- (24) B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, J. Chem. Theory Comput. 4(3), (2008) 435
- (25) G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi, Comput. Phys. Commun. 185(2), (2014) 604
- (26) N. Lu, D. A. Kofke, J. Chem. Phys. 114(17), (2001) 7303
- (27) O. Valsson, M. Parrinello, Phys. Rev. Lett. 113(9), (2014) 090601
- (28) O. Valsson, M. Parrinello, J. Chem. Theory Comput. 11(5), (2015) 1996
- (29) R. W. Zwanzig, J. Chem. Phys. 22(8), (1954) 1420
- (30) M. P. Allen, D. J. Tildesley, Computer simulation of liquids (Oxford university press, 1989)
- (31) P. Kollman, Chem. Rev. 93(7), (1993) 2395
- (32) D. Frenkel, B. Smit, Understanding molecular simulation: from algorithms to applications (Academic press, 2001)