Volume 55. number 3 Cff EIMICAL PHYSICS LETTERS I May 1978 ON A NOVEL MONTE CARLO SCHEME FDR SLMULATING WATER AND AQUEOUS SOLUTiONS* C. PANGALI**, M. RAO and B.J. BERNE CbiumBk Wiiversity, New York, New York fUU27, US.4 Received 3 February I978 The usual hfetropolis ltfonte Carlo afgorithm, when applied to highiy associated liquids iike ST-2 water is _Ehown to be very siow in establishing equilibrium, and often Ieads to bottlenecks in configuration space. A new Jfgorithm is presented in which the Monte Carfo moves are biased in the direction of the forces and torques acting on tbe individual molecule. Comparison with the hfetropoiis scheme shows that this new method is mush more rapidly convergent. En this note we introduce a novel rapid& conver- gent Monte Carlo technique and apply it to the study of water where the common Metropolis procedure is very slowly convergent and often leads to bottlenecks in configuration space. In the Monte Carlo method [ 1,2] , various configu- rations of the system are sampled according to the I3oltzmann distributions exp HW~ t , . ..~~.I @k --- qip (II where Rj is a vector specifying bo*& the positions and o~entations of molecule j. For rigid linear molecufes Ri is a five component vector whereas for rigid non- iinear molecuies Rj is a Six component vector. dRj is the appropriate volume element. To generate a set of contigurations that are dis- tributed according to eq. (I), Metropolis et al. [ 1,2J, devised a scheme according to which a new configu- ration R’ = (R’ t , .-Rh) is generated from an aid con- f~~r~~o~ R = (42 1, -RN) by sampling a prescribed transition probabihty, T(R’lR), and this new con- figuration is accepted with probability * Research supported by National Science Foundation un- der grant #NSF-CHE-76-02414 and NIH-ROi NS 12714- 02. ** In partial fulf&nent of the Ph.D. in Chemistry at Cohen- bia University. and rejected with probab~~ty 4 = I -p_ In tbls way a random walk over con~guratiou space is generated. Averages over the distribution (1) are obtained by weighting old configurations with Q and new con- figurations by p f3] _ Up to this point the method is quite general. Usu- ally the configuration is changed by single particie moves. First particle 1 is moved and this move is either accepted or rejected according to eq. (3) then this is repeated for psrtiele 2,3 . ..K The full cycle of moves is then repeated. Thus aI1 that need be speci- fied is the transition probability for a single particfe move. For atomic fluids, particlei is dispiaced from Rj to .Ri where Ri is randomly sampfed between Ri - $ARo andRj-t-$A.&. The maximum value of the move is $ARu = (fhuu, qAyu. $A+) and the larger this parameter the Isrger the region of configuration space sampled in 3 fixed number of iterations, but AR0 cannot be made arbitrariiy large because most of the moves will then be rejected. The transition probability for the single particle move is then QR;$) = (Ax0 Ayu AZ&‘, i3Rj ED; =O, 6Rj $ D* (3) Here HQ denotes the particular choice of Sx, Sy, and Sz and D denotes the domain defined by the above inequalities, Similarly for moiecuIes, the translational and rotational moves are sampled randomly from 4x3Volume 55, number 3 CHEMICAL PHYSICS LETTERS 1 May 1978 some domain D. For uniform sampling the ratio of the transition probab~~ti~s in eq, (2) is unity. This procedure, correct tl~~.~gh it is, does not con- verge rapidly when applied to strongly angle depen- dent potentials like the water potential where it leads to bottlenecks in configuration space. In fact, previ- ous Monte Carlo and molecular dynamics studies on water do not agree 143 _ A little thought shows that the moves in molecular dynamics are usually biased in the direction of the intermolecular forces, arzd torques, whereas the moves sampled according to eq_ (4) by the usual Metropolis prescription gives equal weight to moves par&ei to the forces as to those antiparallel to those forces. It occurred to us that if we built this bias into ?‘(JQ&), rapid con- vergence might result*. We propose to modify the standard Metropolis scheme by biasing the moves to be along the Forces and torques_ If T(RJ-IRi) a e -BVGR’) - IS substituted in- to eq. (2) it is readily seen that the acceptance prob- ability with be unity_ Of course this procedure is not feasible since it presupposes knowledge of the very distribution required_ Expansion of V(R’) around R, gives for the one particle move ZV$V$.) = c exp C_P[O,i V(R)] *“R/3, VSRi ED; = 0, otherwise, (4) where c is the normalization constant and where V, _ is the gradient operator with respect to the position6 and angles. For mouatom~c fluids, - VR _ Jf(R I, . . ..R~) is the force fj on particle j_ T then simdlifies to c exp (Pr;;:-SRi>. This clearly shows that displacements are more often made ir, the direction of the force than otherwise_ Detailed balance however is still sat- isfied. For polyatomic molecules 7’ = c exp /3(+55 + Nj.8 a/> where fj, NP 6~ and 6wi are respectively the total force, torque, center of mass displacement * In several previous papers, various attempts have been made to improve the efficiency of the Monte Cado technique. In the study of many-body quantum system ES] a force biasing technique was used_ Because T was chosen to be linear in F-6r, and because N particle moves were made, the acceptance ratio decreased by a factor of ten, and the method was not pursued. The transition probability has al- so been modified in various attempts to use importance sampling (see ref. f2) ) and to simulate the gas-liqrid inter- face (see ref. [6 1). 414 and angular displacement of molecule j. The precise form of N~*~~~ depends, of course, on the particular angles used to defme the molecular orientations. in a more extensive paper we show how to appIy eq. (3) to 2 wide variety of representations including simuia- tions involving reshaping using holonomic constraints 171, and simulations