A new algorithm for Reverse Monte Carlo simulations Fernando Lus B. da Silva, Bo Svensson, Torbjörn Åkesson, and Bo Jönsson Citation: The Journal of Chemical Physics 109, 2624 (1998); doi: 10.1063/1.476861 View online: http://dx.doi.org/10.1063/1.476861 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/109/7?ver=pdfcov Published by the AIP Publishing This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 186.217.234.225 On: Tue, 14 Jan 2014 15:11:34 http://scitation.aip.org/content/aip/journal/jcp?ver=pdfcov http://oasc12039.247realmedia.com/RealMedia/ads/click_lx.ads/www.aip.org/pt/adcenter/pdfcover_test/L-37/586982248/x01/AIP-PT/JCP_CoverPg_101613/aipToCAlerts_Large.png/5532386d4f314a53757a6b4144615953?x http://scitation.aip.org/search?value1=Fernando+Lu�s+B.+da+Silva&option1=author http://scitation.aip.org/search?value1=Bo+Svensson&option1=author http://scitation.aip.org/search?value1=Torbj�rn+�kesson&option1=author http://scitation.aip.org/search?value1=Bo+J�nsson&option1=author http://scitation.aip.org/content/aip/journal/jcp?ver=pdfcov http://dx.doi.org/10.1063/1.476861 http://scitation.aip.org/content/aip/journal/jcp/109/7?ver=pdfcov http://scitation.aip.org/content/aip?ver=pdfcov JOURNAL OF CHEMICAL PHYSICS VOLUME 109, NUMBER 7 15 AUGUST 1998 This a A new algorithm for Reverse Monte Carlo simulations Fernando Luı́s B. da Silva,a) Bo Svensson, Torbjörn Åkesson, and Bo Jönsson Physical Chemistry II, Chemical Centre, P.O. Box 124, Lund University, S-221 00 Lund, Sweden ~Received 3 February 1998; accepted 12 May 1998! We present a new algorithm for Reverse Monte Carlo~RMC! simulations of liquids. During the simulations, we calculate energy, excess chemical potentials, bond-angle distributions and three-body correlations. This allows us to test the quality and physical meaning of RMC-generated results and its limitations. It also indicates the possibility to explore orientational correlations from simple scattering experiments. The new technique has been applied to bulk hard-sphere and Lennard-Jones systems and compared to standard Metropolis Monte Carlo results. ©1998 American Institute of Physics.@S0021-9606~98!50731-3# an n he he e y th ha es ce in t th is y o ol a es h ns fo b th te ar M b di s n it- om ul- al it ions ral s and of ions en- on- eri- in ies. put r, it pos- ces, h- mo- ite of n- or- s h e a nd n- om- Fi- ns s n- I. INTRODUCTION The Reverse Monte Carlo~RMC! technique1–6 has emerged as a tool for determining the structure of liquids glasses.2,3 The experimentally determined radial distributio functions~rdfs!, or structure factors, are used as input in t simulation. The particles of the simulated system are t displaced using a Monte Carlo~MC! procedure to generat configurations that correspond to the input~experimental! rdf. During the simulation other structural information ma be collected. In recent years, a number of applications of RMC technique to atomic and molecular systems appeared.1–8 RMC is appealing in that various other sourc of information may be taken advantage of. For instan knowledge in terms of bond lengths, bond angles, coord tion numbers and distance constraints may be used in simulation to complement the input data. However, while RMC simulation generates many-body correlations there question of how relevant they are to the experimental s tem. The aim of this work is to present a detailed analysis the three-body correlation functions generated by Metrop Monte Carlo~MMC! and RMC. In addition, we describe new RMC algorithm with improved convergence properti We will focus on bulk fluids of spherical particles wit interactions given by pairwise additive central forces. Eva9 noted that for such a system, the pair potential, and there all higher ordered distribution functions, are determined the rdf. To confirm that the RMC technique generates correct three-dimensional structure, a simple numerical was proposed. First, a conventional MMC simulation is c ried out and the generated rdf is used as input in an R simulation. The three-body correlation functions obtained the two simulations are then compared. Previous stu along this route have focused on the comparison of the called bond angle and coordination distributio functions.2,4,7,8As we will show, however, these are not su able for a critical comparison. Instead, the deviations fr the Kirkwood superposition approximation10 at the three par- ticle level appear more useful for a careful examination. a!Also affiliated with Departamento de Quı´mica, Fac. de Cieˆncias, Univer- sidade Estadual Paulista, Unesp 17033-360 Bauru, SP, Brazil. 2620021-9606/98/109(7)/2624/6/$15.00 rticle is copyrighted as indicated in the article. Reuse of AIP content is sub 186.217.234.225 On: Tue d n e s , a- he e a s- f is . re y e st - C y es o A number of more or less serious convergence diffic ties have been reported for the RMC technique. In gener seems that better convergence is obtained with simulat using large system sizes, typically on the order of seve thousand particles.2,6 In difficult cases, annealing procedure have been proposed as a way to improve convergence avoid getting trapped into local minima. Another feature the algorithm used is that large errors are obtained in reg of short distances. Thus, at small particle separations, imp etrable hard-core conditions are introduced to exclude c figurations with too close encounters.6,8 Some of this may be of little consequence in practical applications, as the exp mentally determined distribution functions inevitably conta various sources of errors which may lead to inconsistenc Thus, one has to allow for some difference between the in and generated structure. In theoretical tests, howeve should be important to obtain as good an agreement as sible. In this work, we will present a modified algorithm which gives excellent convergence even for short distan and it allows one to avoid the introduction ofad hochard- core constraints. Comparisons between the different simulation tec niques are almost entirely based on structural data. For lecular systems in particular, such tests may be qu involved.11 As a simple complement we suggest the use thermodynamic functions as well, e.g., configurational e ergy and excess chemical potentials. The latter is an imp tant bulk property which implicitly contains contribution from the higher order distribution functions. Armed wit these tests, one can investigate if the RMC results hav physical meaning. As model systems we will use bulk hard-sphere a Lennard-Jones fluids. At first, the new algorithm and its co vergence properties are presented. This is followed by c parisons of the energy and excess chemical potentials. nally, an analysis of various three-body correlation functio is made. II. A NEW RMC ALGORITHM The new RMC algorithm that we propose here follow the original idea.1 The principle is to generate a set of co 4 © 1998 American Institute of Physics ject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: , 14 Jan 2014 15:11:34 ion he - ng - e a- te e ein ie ns is am is s - n o be tio on al - , m nd ion s la- ore ari- ed uc- al- tial. ed he to C. ne va- of us d. r vi- d in d ms le at in on- he sis an ter ri- ith that ter. e d. ed a 2625J. Chem. Phys., Vol. 109, No. 7, 15 August 1998 da Silva et al. This a figurations which corresponds to an input radial distribut function or structure factor. Instead of trying to minimize t free energy like MMC,12,13 we minimize the difference be tween the calculated and experimental rdf. The followi notation is used:ge(r ) is the experimental~input! rdf, gs(r ) is the simulated one, andga(r ) the final answer, the corre sponding histograms used to create these rdfs will app with the same subscript—hists(r ) and hista(r ). We can sum- marize our procedure as follows: ~1! As in ‘‘traditional’’ simulations, we start with an initial configuration of particles in a cubic box. This configur tion could be either generated from random coordina or starting from a lattice. It is of importance that th density is the same as in the experimental system b studied. Usual periodic boundary conditions are appl in all directions. To determine the atomic separatio the minimum image convention is assumed. At th stage, we check the distances and build the histogr hists(r ), which contains the number of atoms at a d tance betweenr and r 1Dr from a central atom. The number of observations,nobs, is equal toN21, whereN is the number of atoms. We duplicate hists(r ) in two other histograms, histsnew (r ) and histsold (r ). ~2! We attempt single particle random displacements a standard MMC.12,13 Using the new configuration, we re calculate the distances and add the new observatio histsnew (r ). The original ~old! distances are added t histsold (r ). Two trial rdfs are calculated as, gsnew ~r !5 histsnew ~r ! nobs4pr 2Drr ~1! and gsold ~r !5 histsold ~r ! nobs4pr 2Drr , ~2! wherer is the number density of the system. It should noticed that both histsnew (r ) and histsold (r ) are based on accumulative distance observations since the simula beginning, and contain the same number of observati ~3! These are compared to the experimentalge(r ), xnew 2 5( i 51 nl @gsnew ~r i !2ge~r i !# 2 ~3! and xold 2 5( i 51 nl @gsold ~r i !2ge~r i !# 2, ~4! where nl is the number of points in the experiment ge(r i). The new configuration is accepted ifxnew 2