2 Monte Carlo Simulation of Water
Downloaded by COLUMBIA UNIV on February 15, 2015 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch002
C. S. PANGALI, M. RAO, and B. J. BERNE Columbia University, New York, NY 10027
The usual Metropolis Monte Carlo(1,2) procedure when applied to the study of water at low temperatures is found to lead to bottlenecks in configuration space. In the Monte Carlo (M.C.) method different configurations of the system are sampled according to the Boltzmann distribution: where Rj is a six component vector in the case of water. Metropolis et al. devised a scheme according to which a new configuration R' = (R'1,...R'N) is generated from an old configuration R = (R1,...RN) by sampling a transition probability T(R'|R), and R' is accepted with probability
and rejected with probability q = 1 - P. Normally the configuration is changed by selecting a particle j at random, or cyclinically, and displacing it from Rj to R'j where dR.j = Rj - Rj is sampled uniformly between —— and + ——. The transition probability for an atomic fluid is therefore given by
where D defines a domain about each point R. in which the moves are to be sampled. Such a method when suitably modified for polyatomic fluids does not converge rapidly, particularly at low temperatures. In fact the M.C. method gives higher potential energies and lower values for the specific heat than the 0-8412-0463-2/78/47-086-029$05.00/0 © 1978 American Chemical Society In Computer Modeling of Matter; Lykos, P.; ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
30
COMPUTER MODELING OF MATTER
molecular dynamics method f o r the same state(3_). I f we —8 v (R ) T(R!|R_.) « , i t i s seen from Eq. (2) t h a t p = 1.
take
1
e
v(R')
i s not known, we
expand i t around C
T (R'. IR.)
e
-BV^ ~R.
R, and
, OR. e D otnerwise
v(R)«6R. ~ ~j
3
=
Since
so
where C i s a n o r m a l i z a t i o n constant and V
(4)
i s the g r a d i e n t
R
J
Downloaded by COLUMBIA UNIV on February 15, 2015 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch002
operator with r e s p e c t to p o s i t i o n s and angles. I t i s easy to see t h a t t h i s procedure s a t i s f i e s the c o n d i t i o n o f d e t a i l balance: e"
B v ( R )
T(R'|R) P(R'|R) =
e
"
B
v
(
R
,
)
T(R|R«) P(R|R')
(5)
We adopted the Barker-Watts procedure(4) f o r performing the u s u a l M e t r o p o l i s M.C. walk. The change i n o r i e n t a t i o n i s accomplished by s e l e c t i n g the x, y or z a x i s a t random and performing a r o t a t i o n through an angle 69 about the chosen a x i s . The angle A0 A9 69 i s u n i f o r m l y d i s t r i b u t e d on (—— , ~^~~) • Our new f o r c e - b i a s 0
0
method samples t r a n s l a t i o n a l displacements 6x., 6y. and 6z. according to the components of the f o r c e F. acting^on the p a r t i c l e j , while 69. w i l l be sampled s u b j e c t to the component of the torque along t h e a x i s s e l e c t e d f o r the r o t a t i o n . 3
3
Results and D i s c u s s i o n S t a r t i n g from a l a t t i c e c o n f i g u r a t i o n with 216 H O molecules i n t e r a c t i n g with an ST-2 potential(5) a t T = 283 K, the standard M e t r o p o l i s method generated a walk i n which the p o t e n t i a l energy p e r ^ p a r t i c l e f l u c t u a t e d around -130.5e ± 0.2e (e = 5.31 x 10 e r g ) . A l t o g e t h e r n e a r l y 500,000 moves were attempted with t h i s method. S t a r t i n g with one of the " e q u i l i b r a t e d " s t a t e s from t h i s method, our f o r c e b i a s scheme generated a q u i t e d i f f e r e n t walk. The r e s u l t s are summarized i n Table I. Not o n l y does the f o r c e b i a s method g i v e a s t a t e with a lower mean p o t e n t i a l energy, but the s t a t e i s r e p r o d u c i b l e . S t a r t i n g with a con f i g u r a t i o n generated by M.D. a t 283°K with V = -135.5e, the standard scheme generates a walk i n which = -135.0c ± 0.2e whereas the forge b i a s method y i e l d s = -134.0e ± 0.4e. At T = 391 K the two M.C. methods and the molecular dynamics method a l l give the same e q u i l i b r i u m s t a t e with s t a t i s t i c a l l y i n d i s t i n g u i s h a b l e p o t e n t i a l energies. In the molecular dynamics method the p a i r f o r c e i s s e t equal to zero when the p a i r d i s t a n c e exceeds a c e r t a i n d i s t a n c e r . Thus the t r a j e c t o r i e s correspond to a s h i f t e d p o t e n t i a l Q
In Computer Modeling of Matter; Lykos, P.; ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
2.
Monte Carlo Simulation of Water
PANGALI ET AL.
31
(J)(r) - ^ ^ Q ) where
