Dynamics and thermodynamics of decay in charged clusters
Mark A. Miller, David A. Bonhommeau, Christian P. Moerland, Sarah J. Gray, Marie-Pierre Gaigeot
aa r X i v : . [ phy s i c s . a t m - c l u s ] A p r April 5, 2018 Molecular Physics preprint
To appear in
Molecular Physics
Vol. 00, No. 00, Month 2015, 1–10
RESEARCH ARTICLEDynamics and thermodynamics of decay in charged clusters
Mark A. Miller a ∗ , David A. Bonhommeau b , Christian P. Moerland c , Sarah J. Gray a andMarie-Pierre Gaigeot da Department of Chemistry, Durham University, South Road, Durham, UK b GSMA CNRS UMR7331, Universit´e de Reims Champagne-Ardenne, UFR SciencesExactes et Naturelles, Moulin de la Housse, BP 1039, F-51687 Reims Cedex 2, France c Department of Applied Physics, Technische Universiteit Eindhoven,Postbus 513, 5600 MB Eindhoven, The Netherlands d Universit´e d’Evry val d’Essonne, LAMBE CNRS UMR8587, Boulevard F Mitterrand,Bˆatiment Maupertuis, F-91025 Evry, France ( Received 00 Month 200x; final version received 00 Month 200x )We propose a method for quantifying charge-driven instabilities in clusters, based on equi-librium simulations under confinement at constant external pressure. This approach makesno assumptions about the mode of decay and allows different clusters to be compared on anequal footing. A comprehensive survey of stability in model clusters of 309 Lennard-Jonesparticles augmented with Coulomb interactions is presented. We proceed to examine dynamicsignatures of instability, finding that rate constants for ejection of charged particles increasesmoothly as a function of total charge with no sudden changes. For clusters where many par-ticles carry charge, ejection of individual charges competes with a fission process that leadsto more symmetric division of the cluster into large fragments. The rate constants for fissiondepend much more sensitively on total charge than those for ejection of individual particles.
Keywords: charged clusters; instability; simulation; Lennard-Jones; Rayleigh limit
1. Introduction
When a droplet or cluster carries a sufficiently high net electrostatic charge, thecharge drives decay into two or more fragments. The first analysis of this phe-nomenon, due to Rayleigh [1] is now well over a century old, but the resultingprediction for the critical charge at which decay occurs retains its central impor-tance in modern research in the area. Rayleigh’s analysis is based on a uniformlycharged, structureless sphere of diameter D that may undergo shape fluctuations,leading simultaneously to an increase in surface energy and a decrease in electro-static energy. Deformation leads to a barrierless decrease in the sum of these en-ergies beyond the limiting charge Q R = π p ǫ γD , where γ is the surface tensionof the droplet interface and ǫ is the permittivity of free space. Rayleigh reachedthis classic result within the first page of his elegantly short article [1], concealinga somewhat involved derivation [2].Since Rayleigh’s time, his work has gained new significance due to the impor-tance of charged droplets in contemporary fields such as atmospheric chemistry [3] ∗ Corresponding author. Email: [email protected]
ISSN: 0026 8976 print/ISSN 1362 3028 onlinec (cid:13) pril 5, 2018 Molecular Physics preprint and electrospray ionisation mass spectrometry [4]. Nevertheless, it has been recog-nised that unstable global deformations are not the only mechanism by which acharged droplet may decay. For example, the ion evaporation model [5] provides acontinuum-based understanding of how small ions may be ejected from a dropleteven below the Rayleigh limit. When macroions are dissolved in a droplet, othermechanisms come into play. Evaporation of the droplet may lead to some of thecharge being deposited on the solute [6]. Unravelling of polymeric solutes in watercan lead to more complex processes [7, 8], including the extrusion of a solvent-free segment of the chain or the formation of subdroplets along extended chainconformations [9].In principle, a cluster at finite temperature in an unbounded environment isthermodynamically at most metastable, even if it is charge-neutral, since unlim-ited translational entropy can be gained by dividing the cluster into fragments.Hence, the concept of stability implicitly involves a sense of time scale, such asthe typical time of flight of an electrosprayed droplet in a mass spectrometer. In aprevious paper, some of us devised an equilibrium approach to identifying charge-driven instability in small clusters [10] based on confining the cluster to a sphericalcontainer and monitoring the fraction of equilibrium configurations in which thecluster had emitted at least some of its charged particles. In the present article,we propose a different equilibrium-based approach that provides some significantadvantages in terms of the generality of its applicability and the physical basis ofits interpretation. We then proceed to a quantitative characterisation of the dy-namics of charge-driven instabilities in a simple model, focusing on the competitionbetween ejection and fission processes.
2. Model
As in our previous study [10], we use a Lennard-Jones (LJ) cluster of N = 309particles to explore methods of characterisation. For simplicity, all particles interactwith the same LJ well-depth u and length parameter σ . In addition, particles maycarry a charge, introducing a Coulomb term to the total potential energy: E = 4 u N X i 43, which lies just above the broadened melting transition,for most of our calculations. In this liquid-like regime, the highly ordered globalpotential minimum structure is not seen. As might be expected, the addition ofcharge lowers the cluster’s melting temperature [13].The LJ potential with charges localised on individual particles is clearly a highlyidealised model. It does not incorporate the complexity associated with hydrogenbonding in water clusters or mechanisms for transferring excess charge to the clus-ter surface by orientation of dipole moments [14]. We note that it is possible torefine the model by the inclusion of polarisability, with induced dipoles solved self-consistently [13], though this is a costly addition to a simple model. Nevertheless,idealised models of this type do provide scope for probing different mechanisms ofdecay [15, 16], and can even provide a useful reference for quantitative analysis ofexperimental data on the evaporation of real charged droplets [17]. 3. Droplets at equilibrium As pointed out in the introduction, an unbounded cluster decays over time even ifit is uncharged. Clusters bound by van der Waals-like interactions such as the LJpotential gradually evaporate, typically one particle at a time. Hence, the equilib-rium thermodynamics of clusters must always be characterised within a container.A spherical enclosure is usually used, and an appropriate radius must then be cho-sen by some criterion [18, 19]. In our previous article [10], we used a containerradius of R c = 6 . σ for the 309-atom LJ cluster, with the centre of the containertracking the centre-of-mass of the cluster [18]. This choice of R c allows a clear sep-aration between any particles that have been decisively ejected and the remainingsubcluster. However, the equilibrium between evaporated and condensed particlesis influenced by the choice of radius [19] and so the arbitrary choice of R c is anunattractive feature of the approach. It is also difficult to compare clusters of dif-ferent size N , since it is not entirely clear how R c should change with N for a trulyunbiased comparison.In the present work, we take a different approach that offers a number of ad-vantages. We place the cluster in a fluctuating container whose volume changes inresponse to the cluster inside it and to a fixed external pressure p . This pressuremay be set at a physically meaningful value by reference to, for example, the back-ground medium or carrier gas in a particular application. The same value of p maystraightforwardly be used for different N if a comparison between cluster sizes isrequired.The isobaric, isothermal approach that we propose is slightly different from typi-cal constant- N pT Monte Carlo simulations, which are normally encountered whenmodelling bulk systems as represented by a periodic simulation cell. In a periodicsystem, stochastic changes in the cell volume are accompanied by a uniform scalingof all particle positions in order to respect the boundary conditions in the resizedcell, and particle displacement steps are made in coordinates that have been scaledto lie in a unit cube. This protocol leads to a Monte Carlo acceptance criterion fortrial volume changes of [20] P acc = min " , exp (cid:18) − ( E ′ − E ) k B T (cid:19) exp (cid:18) − p ( V ′ − V ) k B T (cid:19) (cid:18) V ′ V (cid:19) N , pril 5, 2018 Molecular Physics preprint 25 30 35 40 45 50total charge, Q *681012 m ea n c on t a i n e r r a d i u s , 〈 R c * 〉 n = 10 n = 30 n = 60 Figure 1. Detection of decay using the fluctuating-container method at reduced pressure p ∗ = 0 . Q ∗ distributed amongst three different numbers n of charge-carriers. Instabilityappears as a sharp increase in the reduced mean container radius. where V is the volume of the cell and primes indicate new values after the trialmove. The factor containing E is due to the energy change associated with scal-ing the particle coordinates, while the factor ( V ′ /V ) N comes from the change ofvariables to scaled coordinates in a unit cell. However, our cluster is not a peri-odic system, and it is more efficient to allow the container to fluctuate withoutunnaturally expanding or contracting the cluster by a scaling of the particle co-ordinates. If the configuration is left unchanged during trial volume changes thenthe acceptance for trial volume changes is the simpler P acc = min " , exp (cid:18) − p ( V ′ − V ) k B T (cid:19) N Y i =1 Θ( R ′ c − | r ′ i | ) , where r ′ i is the position of particle i relative to the centre of mass (after the trialmove). The product of Heaviside step functions Θ enforces the requirement that allparticles lie within the container. Hence, any trial decrease in volume that wouldresult in the exclusion of an outlying particle must be rejected. In effect, sucha proposed state would have infinite energy because the container is hard. Trialvolume changes constitute 1% of Monte Carlo steps in our simulations.In addition to the modified volume-change steps, 20% of Monte Carlo steps inour simulation are devoted to attempted exchanges between neutral and chargedparticles, leaving their positions fixed. These swaps are accepted according to thestandard Metropolis criterion [20, 21] and make exploration of the equilibriumconfigurations more efficient [10].Figure 1 shows the mean reduced container radius h R ∗ c i as a function of totalcharge Q ∗ at a reduced pressure of p ∗ = pσ /u = 0 . p ∗ ≈ . × Monte Carlo steps per particle. A sharpincrease in h R ∗ c i occurs at a particular charge, indicating that the fragmentation ofthe cluster is sufficient to force the container outwards against the external pressure.From a thermodynamic point of view, the sudden change is akin to the finite-systemanalogue of a first-order liquid–gas phase transition, but here driven by chargerather than by temperature. The value of Q ∗ at which the jump occurs increaseswith the number of particles n over which the charge is divided, reinforcing ourearlier result that the maximum charge a cluster can sustain is generally largerwhen the charge is divided over more particles [10].4 pril 5, 2018 Molecular Physics preprint n m a x i m u m c h a r g e , Q * m a x p * = 0.002 p * = 0.005 p * = 0.01 Figure 2. Maximum charge as detected by the fluctuating-container method at three reduced pressures p ∗ . We take the point at which the numerical derivative dR ∗ c /dQ ∗ is largest to be themaximum charge Q ∗ max that the cluster can sustain. This transition is more sharplydefined than the point at which dissociated configurations become a majority in acontainer of fixed radius [10]. Q ∗ max based on the fluctuating-container method isplotted as a function of n in Fig. 2, showing qualitatively similar behaviour at threedifferent pressures. At large n , Q ∗ max shows a small but systematic decrease. As willbe explored in the next section, this feature corresponds to the onset of competitionbetween ejection of individual charged particles and fission of the cluster into largerfragments. The fact that the fluctuating container detects both these modes ofdecay using the same procedure and without prior knowledge of the change inmechanism is one of the method’s appealing features. Even in the more complexcases of droplets of polar molecules with or without a polymeric solute, the decaypathways all involve spreading out of the charged particles by a decisive deviationof the droplet away from a spherical or near-spherical shape [7, 8]. Such changeswould always lead to a sudden increase in the radius of the container at somecritical charge. Hence, the fluctuating-container method should be quite generallyapplicable to the problem of detecting and defining charge-driven instability.At very small n , Fig. (2) shows that there is a shallow minimum in Q ∗ max . Thisfeature seems to be related to the rapid increase in energy with respect to thenumber of charges at small n . In the hypothetical case where the charged parti-cles lay neatly at the surface of a perfectly spherical cluster, the cluster would bedescribed by the so-called Thomson problem [22], which poses the question of theoptimal arrangement and corresponding energy of n unit charges on the surface ofa unit sphere. In the continuum limit, a sphere with a uniformly charged surfacehas an electrostatic energy that scales as the square of the surface charge density.The global potential energy minima of the Thomson problem do approach a con-stant ratio of energy to n at large n [23]. However, for small n , the energy risesconsiderably more steeply than n [24]. 4. Dynamics of decay We now turn to dynamic signatures of charge-driven instability. The clusters mustbe prepared in a well defined intact state before being allowed to decay. Our proto-col is to equilibrate the cluster using constant-temperature Monte Carlo simulationsin a tight-fitting spherical container of fixed radius 5 . σ , which is the maximum5 pril 5, 2018 Molecular Physics preprint t *01234567 l n M Q * = 40506070 Figure 3. First-order decay plots for first ejection of a charged particle averaged over an ensemble of 1000clusters carrying n = 100 charges at four values of the initial total charge Q ∗ . M is the number of intactclusters remaining at time t ∗ . radius of the icosahedral global minimum structure [11] plus a margin of σ . Thecontainer is then removed and random Maxwell–Boltzmann-distributed velocitiesare assigned to the atoms. From this point, the cluster is allowed to evolve viaconstant-energy molecular dynamics (MD) with a standard Verlet integrator [20]and a time-step of δt ∗ = 0 . t ∗ = t p u/mσ and m is the mass of one particle.To detect evaporation and ejection events, we monitor the distance of all particlesfrom the centre of mass. When a particle reaches a large distance (22 σ ) from theremaining subcluster, we can be confident that it will not return. The actual timeof the evaporation or ejection is then backdated to the MD step at which thevelocity of the particle first started on its outward trajectory, i.e. , the earliest time t ∗ e such that v i ( t ∗ ) · r i ( t ∗ ) > t ∗ ≥ t ∗ e , where v i is the velocity of particle i . This procedure allows temporary excursions of particles to be ignored while stillproviding an accurate timestamp for ejection events.We have measured the time of first ejection of a charged particle in an ensembleof 1000 independent trajectories. From the list of ejection times, we may constructa first-order kinetics plot of the logarithm of the number M of clusters that havenot yet decayed as a function of time. Results for a selection of ensembles with n = 100 charged particles and different initial total charges are shown in Fig. 3,confirming that ejection of the first charge is a first-order process. From the slopes,we obtain reduced rate constants k ∗ , which are plotted (circles) as a function ofinitial charge in Fig. 4.The total charge spanned in Fig. 4 starts well below and finishes considerablyabove the range of Q ∗ max measured in our equilibrium simulations (Fig. 2) for n = 100. The rate constant increases smoothly with Q ∗ , never showing a decisivejump that could be taken as a dynamically-defined threshold for instability.We may also extract an effective activation energy E ∗ a for ejection from Arrhe-nius plots of ln k ∗ against inverse temperature. Remembering the proximity of ourchosen temperature T ∗ = 0 . 43 to the melting transition of the neutral cluster [12]on the one hand, and the rapidly increasing tendency for thermal evaporation athigher temperatures [19] on the other, we have varied 1 /T ∗ over only a narrowrange. Arrhenius behaviour is nevertheless clearly evident, as shown in the inset ofFig. 5, allowing the activation energy to be plotted as a function of total charge inthe main panel of the figure. In keeping with the smoothly changing rate constantsof Fig. 4, the activation energy shows no sudden changes with charge and appears6 pril 5, 2018 Molecular Physics preprint 30 40 50 60 70 80 90initial charge on cluster, Q * −3 −2 −1 r e du ce d r a t e c on s t a n t , k * Figure 4. Rate constants for the ejection of the first particle for a cluster carrying n = 100 chargedparticles (blue circles) and n = 309 charged particles (red squares), and rate constants for fission of the n = 309 cluster (green plusses). 40 45 50 55 60 65 70total charge, Q *23456 r e du ce d ac ti v a ti on e n e r gy , E a * T * −4−3−2−1 l n k * Figure 5. Activation energy E ∗ a for particle ejection for a cluster carrying n = 100 charges. Bars show thestandard error in the linear regression coefficient. The inset shows example Arrhenius plots whose slopesare − E ∗ a : Q ∗ = 40 (blue circles), Q ∗ = 50 (red squares), Q ∗ = 60 (green diamonds). to depend rather linearly on charge in this case.In highly charged clusters, the ejection of a charged particle can lower the overallpotential energy of the remaining subcluster and, depending on the kinetic energycarried away, can in turn lead to slight increase in the effective temperature of thesubcluster. However, we have not observed any “run-away” cascades where the rateof decay of an ensemble increased with successive ejection events.The same observation of smoothly changing first-ejection rates at n = 100 appliesto the case of n = 309 (squares in Fig. 4), where the charge is spread over all theparticles in the cluster. As might be predicted from the higher value of Q ∗ max atlarge n (Fig. 2), the rate constant is lower than for the n = 100 case at a giventotal charge, despite the fact that there are more charged particles that could beejected at n = 309.The crucial difference between n = 100 and n = 309 is that ejection of individualcharged particles is not the only mechanism of decay for the latter. At large n ,clusters may also undergo a more symmetric fission process, in which two largefragments are produced. Snapshots of ejection and fission processes are shown inFig. 6. Fission may be detected by performing a cluster analysis based on a pairwiseparticle connectivity criterion of 1 . σ and defining fission to have taken place at7 pril 5, 2018 Molecular Physics preprintFigure 6. Left: ejection of a particle from a cluster carrying n = 10 charges (red spheres), total charge Q ∗ = 24. Right: fission of a cluster of n = 309 charges, Q ∗ = 67, into two sub-clusters. the first time at which two separate clusters, each of at least N/ n = 309 and the absence offission at smaller n can be understood from an idealised analysis by Consta andMalevanets of a spherical droplet splitting into two daughter droplets with con-servation of total volume and of total charge [2]. Symmetrical decay (fission) isfavoured when electrostatic energy dominates over surface energy, since the de-crease in electrostatic energy of separating two highly charged subclusters com-pensates for the large increase in surface area caused by creating two spheres fromone with the same total volume. This is the regime in our clusters at large n . Incontrast, when n is small, a fractionally significant decrease in electrostatic energycan be achieved by the ejection of a single charged particle without incurring alarge increase in surface energy.Our procedure for determining the decay rate constants for ejection and fissionrelies on being able to separate trajectories cleanly into portions corresponding toreactants (intact clusters) and products (separated particles or fragments). Themethods of backdating outbound ejected particles and monitoring the connectivityof subclusters in fission amount to selecting reactant–product dividing surfacesin phase space. Our choices of surface are not unique but have the advantagesof being straightforward to evaluate and, crucially, not being prone to recrossingby trajectories. Hence, crossings can be interpreted decisively as decay events. Indroplets of polar molecules, where decay may occur via departure of a small cluster[25], it should be possible to adapt the fission criterion to detect division of thedroplet into uneven fragments simply by adjusting the threshold in the number of8 pril 5, 2018 Molecular Physics preprint particles at which a detached fragment is defined.To gain more detailed insight into the fluctuations leading to decay is a moredifficult goal that requires construction of a reaction coordinate capable of faith-fully describing the ensemble of decay pathways. The task has been approached byConsta and coworkers [25, 26] for the case of water droplets containing simple ions.Such clusters develop a bottleneck as the ions start to separate prior to decay of thedroplet, even when decay is uneven with respect to fragment size. A special trans-fer reaction coordinate was introduced to identify the resulting dumbbell-shapedconfigurations, as distinct from more general prolate distortions. This approachallowed differences to be observed in the diffusive nature of trajectories near thetop of the free energy barrier for clusters containing different ions [26]. As in anyactivated process, the height and shape of the free energy barrier for a given pro-cess and system depend on the choice of reaction coordinate, but the overall rateconstant should not, provided that a converged transmission coefficient can beobtained. 5. Concluding remarks Our fluctuating-container simulations effectively cast charge-driven instability as atransition under well defined equilibrium thermodynamic conditions. This approachhas the advantages of not requiring any assumptions about the mechanism of theinstability, of permitting comparisons between different clusters if required, andof being quite sharply defined even for small clusters (where temperature-driventransitions are normally broadened). The method is also oblivious to any transientevaporation of neutral fragments that are of no interest in the context. The onlyparameter that must be fixed is the imposed pressure, but this quantity does havea direct physical interpretation.It is not straightforward to compare our results for the stability of the 309-atom LJ cluster with the prediction of Rayleigh’s classic formula, not least becausethe surface tension is not known. Clusters in this size range have a dramaticallydepressed melting temperature in comparison with the bulk, and the liquid-likeregime that we have studied around T ∗ = 0 . 43 is well below the triple point ofbulk Lennard-Jones, which lies at about T ∗ = 0 . 69 [27, 28]. Below this point, nobulk equilibrium value of γ exists.The first mode of deformation to become unstable in Rayleigh’s model is thesecond-rank spherical harmonic, corresponding to oblate–prolate deviations fromthe sphere. Although Rayleigh’s work makes no formal prediction about the mecha-nism or products that result from the instability, we note that a prolate deformationmuch more closely resembles the motion that precedes the fission mechanism thatwe observe in the LJ cluster at large number n of charges. Nevertheless, even wheresingle-particle ejection dominates the decay and fission is absent, the clusters aretypically distorted in a prolate direction.For simplicity, we have given all particles in our cluster the same LJ energy andlength parameters. However, if charged particles interacted with stronger well depth u or had a larger diameter σ , one could envisage emission of “solvated” chargesrather than bare particles. Clearly, a number of other refinements could be includedto make the model a more realistic depiction of a particular system. However, theidealised model has allowed a systematic survey of charge limits, the decay of largeensembles, rate constants, and activation barriers to be performed. The methodsdeployed here should be useful in future investigations of other charged clustersand droplets. 9 pril 5, 2018 Molecular Physics preprint Acknowledgements SJG is grateful to Durham University’s Institute for Advanced Research Computingfor a summer seedcorn grant. DAB is indebted to CNRS (Centre National de laRecherche Scientifique) for the award of a Chaire d’excellence. Travel associatedwith this collaboration was supported by the Alliance Programme of the BritishCouncil and the French Minist`ere des Affaires Etrang`eres. The authors thank DrFlorent Calvo for helpful discussions. References [1] Lord Rayleigh, Philos. Mag. , 184 (1882).[2] S. Consta and A. Malevanets, Mol. Simul. , 73 (2015).[3] A. Hirsikko, T. Nieminen, S. Gangn´e, K. Lehtipalo, H.E. Manninen, M. Ehn, U. H˜orrak,V.M. Kerminen, L. Laakso, P.H. McMurry, A. Mirme, T. Pet¨aj¨a, H. Tammet, V. Vakkari,M. Vana and M. Kulmala, Atmos. Chem. Phys. , 767 (2011).[4] J.B. Fenn, Angew. Chem. Int. Ed. , 3871 (2003).[5] J.V. Iribarne and B.A. Thomson, J. Chem. Phys. , 2287 (1976).[6] M. Dole, L.L. Mack, R.L. Hines, R.C. Mobley, L.D. Ferguson and M.B. Alice, J. Chem. Phys. , 2240 (1968).[7] S. Consta, J. Phys. Chem. B , 5263 (2010).[8] E. Ahadi and L. Konermann, J. Phys. Chem. B , 104 (2012).[9] S. Consta and A. Malevanets, J. Chem. Phys. , 044314 (2013).[10] M.A. Miller, D.A. Bonhommeau, C.J. Heard, Y. Shin, R. Spezia and M.P. Gaigeot, J. Phys.Cond. Mat. , 284130 (2012).[11] D. Romero, C. Barr´on and S. G´omez, Comp. Phys. Comm. , 87 (1999).[12] E.G. Noya and J.P.K. Doye, J. Chem. Phys. , 104503 (2006).[13] D.A. Bonhommeau and M.P. Gaigeot, Comp. Phys. Comm. , 873 (2013).[14] R. V´acha, S.W. Rick, P. Jungwirth, A.G.F. de Beer, H.B. de Aguiar, J.S. Samson and S.Roke, J. Am. Chem. Soc. , 10204 (2011).[15] I. Last, Y. Levy and J. Jortner, Proc. Natl. Acad. Sci. USA , 9107 (2002).[16] Y. Levy, I. Last and J. Jortner, Mol. Phys. , 1227 (2006).[17] R. Ho lyst, M. Litniewski, D. Jakubczyk, M. Zientara and M. Wo´zniak, Soft Matter , 7766(2013).[18] J.K. Lee, J.A. Barker and F.F. Abraham, J. Chem. Phys. , 3166 (1973).[19] C.J. Tsai and K.D. Jordan, J. Chem. Phys. , 6957 (1993).[20] D. Frenkel and B. Smit, Understanding Molecular Simulation , 2nd ed. (Academic Press, SanDiego, 2002).[21] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. , 1087 (1953).[22] J.J. Thomson, Philos. Mag. Ser. VI , 237 (1904).[23] D.J. Wales and S. Ulker, Phys. Rev. B , 212101 (2006).[24] H.A. Munera, Nature , 597 (1986).[25] S. Consta, J. Mol. Struct. (Theochem.) , 131 (2002).[26] S. Consta, K.R. Mainer and W. Novak, J. Chem. Phys. , 10125 (2003).[27] E.A. Mastny and J.J. de Pablo, J. Chem. Phys. , 104504 (2007).[28] J.P. Hansen and L. Verlet, Phys. Rev. , 151 (1969)., 151 (1969).