Molecular Dynamics Calculations on a Minicomputer

11 Molecular Dynamics Calculations on a Minicomputer

Downloaded by UNIV OF MASSACHUSETTS AMHERST on May 19, 2018 | https://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0057.ch011

PAUL A. FLINN Physics and Metallurgy Departments, Carnegie-Mellon University, Pittsburgh, PA 15213

Since its introduction by Rahman in 1964 (1), the technique of computer simulation of motion in liquids, generally known as "molecular dynamics", has played a vital role in increasing our understanding of the real nature of the liquid state. Applications of the technique to various liquids have been reviewed by McDonald and Singer (2), Rahman (3), and Fisher and Watts (4). The results of the calculations have been in excellent agreement with a variety of experimental measurements of the properties of liquids: the equation of state, the radial distribution function, inelastic scattering of neutrons, and diffusion. The molecular dynamics results also provide valuable tests of the adequacy of various approximate analytic theories of liquids. A major limitation of the technique has been economic: the calculations have required large amounts of time on large, expensive, computers. Fortunately, it is possible to carry out useful molecular dynamics calculations at greatly reduced cost on a minicomputer (or microcomputer); much more widespread use, including instructional use, of the technique, should now be possible. The calculation is, in principle, quite simple; it consists of numerical integration of the simultaneous nonlinear differential equations of motions for a number of particles constituting a small sample of the liquid. In the original work the Newtonian form of the equations of motion was used: d r. 1

m

.5-

f(r..)

dt where m is_^the p a r t i c l e mass, i s the position of the i t h p a r t i c l e , r-y i s t h e d i s t a n c e between t h e centers o f p a r t i c l e s i and j , and f i s t h e f o r c e a c t i n g between p a r t i c l e i and j . F o r the work d e s c r i b e d here, i t was more convenient t o use t h e Hamilton!an equations:

137 Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

138

MINICOMPUTERS

-4 dp. dt " . ^ j

U

r

SCALE

COMPUTATIONS

-• dr ρ dt " m

;

i j

A N D LARGE


where p. i s t h e momentum o f t h e i t h p a r t i c l e . The p o t e n t i a l used f o r a given l i q u i d i s g e n e r a l l y o f a form suggested "by t h e o r e t i c a l arguments, hut with parameters obtained from experimental data on the m a t e r i a l . To i l l u s t r a t e t h e method we use t h e case of argon, w i t h a Lennard-Jones p o t e n t i a l : V(r) =

[(σ/r)

1 2

6

- (σ/r) ]

—21 e = I.65 X 10"

and parameter values work o f Rahman.

5? J , σ = 3.U À taken from t h e

B a s i c P r i n c i p l e s o f Minicomputer Use. The c h a r a c t e r i s t i c f e a t u r e s o f most minicomputers are small word s i z e (l6 b i t s ) , l i m i t e d memory, and reasonable speed f o r integer arithemetic. F l o a t i n g point operations u s u a l l y r e q u i r e subroutines and are quite slow. The usefulness o f minicomputers f o r molecular dynamics c a l c u l a t i o n s r e s u l t s from the f a c t that t h e range o f values o f t h e v a r i a b l e s needed i s s u f f i c i e n t l y l i m i t e d t h a t i n t e g e r a r i t h m e t i c can be used, and l 6 b i t p r e c i s i o n i s adequate. The only p o t e n t i a l d i f f i c u l t i e s a r i s e i n connection w i t h t h e interatomic f o r c e f u n c t i o n , which may be o f complicated form, and has an unbounded magnitude. These problems can be f a i r l y e a s i l y circumvented: t h e interatomic f o r c e f u n c t i o n i s evaluated and t a b u l a t e d a t the beginning o f the c a l c u l a t i o n ; i n the body o f t h e c a l c u l a t i o n determination of t h e f o r c e i s simply a look-up operation. The wide range o f the magnitude of t h e f o r c e does not represent any r e a l problem, since t h e f o r c e becomes i n c o n v e n i e n t l y l a r g e only at d i s t a n c e s considerably shorter than those which a c t u a l l y e x i s t when the l i q u i d i s a t or near e q u i l i b r i u m . We can, t h e r e f o r e , truncate t h e magnitude of t h e f o r c e t o a constant value f o r d i s t a n c e s l e s s than some r . We a l s o , as i s customary, l i m i t t h e range of i n t e r a c t i o n by s e t t i n g the f o r c e equal t o zero f o r d i s t a n c e s beyond t h e c u t o f f d i s t a n c e . Our f o r c e law then has t h e form: s

r < r r

g

^ r * ^

f(r) f

(

r

f(r ) g

-13

) c r

r > r c

f(r) = ' v

c r 2

0

and t h e lookup t a b l e need cover only t h e range r ^ r ^ r . F o r t h i s c a l c u l a t i o n , r was taken as 2.82 A, and r_ as 5.07 A. r

c

Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

11.

FLiNN

Molecular Dynamics Calculations

139


System Hardware. The minicomputer used f o r t h i s work was a Texas Instruments 96OA, borrowed from i t s normal use as a Mossbauer spectrometer and data processor (5>). The system i n c l u d e s 8192 words (l6 b i t ) of semiconductor memory, an i n t e r f a c e t o a t e l e t y p e with paper tape punch and reader, and a CRT d i s p l a y . The o p t i o n a l extended i n s t r u c t i o n set of the 96OA includes the f o l l o w i n g hardware operations: m u l t i p l y two l6 b i t words, 32 b i t (double word) product; d i v i d e double word by s i n g l e word, s i n g l e word quotient and s i n g l e word remainder; double word add; double word subtract; double word l e f t and r i g h t s h i f t operations. The d i s p l a y i s provided by a Tektronix 603 storage d i s p l a y u n i t , d r i v e n by two D a t e l DAC k^lOB 10 b i t analog t o d i g i t a l converters, i n t e r f a c e d through the communications r e g i s t e r u n i t (CRU) of the computer. Alphameric d i s p l a y i s provided by software; no character generating hardware i s used. The o r i g i n a l cost of the computer was about $7000 (1972). The 96OA i s no longer made, but equipment w i t h s i m i l a r performance can now be obtained at a much lower cost. For i n t e g e r a r i t h m e t i c , the 96OA i s only moderately slower than t y p i c a l l a r g e computers, such as the Univac 1108. The times i n microseconds f o r some t y p i c a l i n s t r u c t i o n s are:

96 OA

1108

Add

3.583

1.50

Subtract

3.583

1.50

Multiply

8.583

3.125

10Λ17

Load

3.333

3.875 1.50

Store

3.583

1.50

Divide

System Software. The operating system used f o r t h i s work was one o r i g i n a l l y w r i t t e n f o r the Mossbauer spectrometer a p p l i c a t i o n , and described i n more d e t a i l elsewhere (_5). I t c o n s i s t s of a monitor, i/o r o u t i n e s , and a f l o a t i n g p o i n t a r i t h m e t i c package. The monitor provides f o r the l o a d i n g of programs from paper tape, i n i t i a t i o n of execution, recovery from e r r o r t r a p s , dump, patch, and debug facilities. The i/o r o u t i n e s provide f o r t e l e t y p e input and output of decimal, hexadecimal, and alphanumeric data, and CRT d i s p l a y of alphanumeric data by software character generation. The f l o a t i n g point package i n c l u d e s a d d i t i o n , s u b t r a c t i o n , m u l t i p l i c a t i o n , d i v i s i o n , i n t e g e r t o f l o a t i n g p o i n t , and f l o a t i n g p o i n t t o i n t e g e r conversion. A f i x e d point square root r o u t i n e was w r i t t e n and included f o r t h i s c a l c u l a t i o n . The system


MINICOMPUTERS

140

AND

LARGE SCALE

COMPUTATIONS

software occupies 2128 words, with an a d d i t i o n a l 82 words f o r t h e square root r o u t i n e . A l l programming was done i n assembly language, and converted t o object code by a macro i n s t r u c t i o n processor and a cross assembler run on an IBM 360/67.


Molecular Dynamics Program. The working program c o n s i s t s o f s e v e r a l p a r t s : initial i z a t i o n , c o o l i n g , i n t e g r a t i o n , c a l c u l a t i o n o f s t a t i s t i c s , CRT d i s p l a y , and t e l e t y p e output. The program storage requirements, i n l 6 b i t words, are: i n i t i a l i z a t i o n , 210; c o o l i n g , 38; i n t e g r a t i o n , 192; s t a t i s t i c s , 100; d i s p l a y and output, kQ; t o t a l , 588. The data storage requirements are: f o r c e t a b l e , 102^; p o s i t i o n , momentum, i n i t i a l p o s i t i o n and i n i t i a l momentum, 38+ each; c o r r e l a t i o n f u n c t i o n s , 102^; t o t a l data 358^. The o v e r a l l space r e q u i r e d i s ^172 words. 1

Units and S c a l i n g . I t i s customary i n molecular dynamics c a l c u l a t i o n s t o use a system o f u n i t s based on t h e p r o p e r t i e s o f the system under study: f o l l o w i n g the choice o f Tsung and Maclin (6), we take t h e p a r t i c l e mass as t h e u n i t o f mass, t h e parameter σ as t h e u n i t o f length, and a conveniently short time ( l O " ^ sec.) as t h e u n i t o f time. With t h i s convention, t h e momentum i s numerically equal t o t h e velocity. In order t o avoid unnecessary l o s s o f p r e c i s i o n i n t h e i n t e g e r a r i t h m e t i c , i t i s necessary t o scale t h e v a r i a b l e s o f t h e problem i n t o proper i n t e g e r u n i t s . We consider f i r s t t h e length s c a l e : we have a system o f Ν p a r t i c l e s , with a volume V.per p a r t i c l e , f o r a t o t a l volume To use t h e f u l l p r e c i s i o n o f the computer we represent t h i s length by 2 ^ . We choose our u n i t o f time f o r one i n t e g r a t i o n step as 62.5 femtoseconds; t h i s i s scaled as l / l 6 o f an i n t e g e r u n i t , since m u l t i p l i c a t i o n by At i s accomplished by a s h i f t o f k b i n a r y places t o the r i g h t . One i n t e g e r unit o f time i s t h e r e f o r e 0.1 picosecond. T h i s choice o f distance and time s c a l e s f i x e s t h e v e l o c i t y (and momentum) s c a l e s . We d e f i n e the "temperature" o f the system i n terms o f the k i n e t i c energy: 1

Τ - 3k/m. 2

Startup

of Calculation.

The f i r s t step i n the c a l c u l a t i o n i s the generation o f t h e f o r c e lookup t a b l e f o r t h e range o f R needed (1023^ t o 18^26). Since t h e memory space a v a i l a b l e was quite l i m i t e d , steps o f 8 were used, so that only 102^ l o c a t i o n s were r e q u i r e d . I n t e r p o l a t i o n from the t a b l e was planned f o r intermediate values of R, but proved t o be unnecessary.



11.

FLiNN


141

Next, t h e i n i t i a l c o n f i g u r a t i o n o f t h e system i s constructed by a s s i g n i n g random values t o t h e p o s i t i o n and momentum coordinates of t h e p a r t i c l e s . Random numbers are generated by the m u l t i p l i c a t i v e congruence method: repeated m u l t i p l i c a t i o n by 3125 and r e t e n t i o n o f t h e l e a s t s i g n i f i c a n t h a l f of t h e double word product. T h i s i n i t i a l c o n f i g u r a t i o n has, o f course, an extremely high energy. I t i s necessary t o " c o o l " t h e system by g r a d u a l l y removing k i n e t i c energy. T h i s i s done by p e r i o d i c a l l y reducing the magnitude o f each component o f momentum o f each p a r t i c l e by some f r a c t i o n o f i t s value. The c o o l i n g must be done g r a d u a l l y t o avoid f r e e z i n g i n a nonequilibrium s t a t e ; t h e c o o l i n g r a t e ' s h o u l d not exceed t h e r a t e a t which t h e i n i t i a l l y extremely high p o t e n t i a l energy o f t h e system can be converted i n t o k i n e t i c energy, so t h a t approximate e q u i p a r t i t i o n i s maintained. Main C a l c u l a t i o n Loop. The c a l c u l a t i o n proper i s c a r r i e d out i n a nest of three loops. The outermost (T loop) i s a time loop; each execution corresponds t o one time i n t e g r a t i o n step. The intermediate loop (J loop) i s over a l l p a r t i c l e s ; one execution corresponds t o a c a l c u l a t i o n o f the net f o r c e on one p a r t i c l e , and an updating o f the p o s i t i o n and momentum o f t h a t p a r t i c l e . The innermost loop (I loop) i s a l s o over a l l p a r t i c l e s ; one execution corresponds t o a c a l c u l a t i o n o f t h e f o r c e on one p a r t i c l e due t o one other particle. We use t h e f o l l o w i n g n o t a t i o n t o d e s c r i b e t h e c a l c u l a t i o n : X l ( l ) , X 2 ( l ) , X 3 ( l ) : p o s i t i o n coordinates of the I t h p a r t i c l e . P l ( l ) , P 2 ( l ) , P 3 ( l ) : momentum coordinates o f t h e I ' t h p a r t i c l e . DX1, DX2, DX3: components o f t h e v e c t o r from p a r t i c l e J t o p a r t i c l e I; e.g., DX1 = X l ( l ) - X l ( j ) . R: t h e d i s t a n c e from p a r t i c l e J t o p a r t i c l e I . F: t h e f o r c e exerted by p a r t i c l e I on p a r t i c l e J . F l , F2, F3: t h e components o f the net f o r c e on p a r t i c l e J ; t h i s i s c a l c u l a t e d as a running sum over t h e I p a r t i c l e s i n t h e inner loop. The c a l c u l a t i o n proceeds as f o l l o w s : Zero t h e time r e g i s t e r and enter t h e Τ loop. I n i t i a l i z e t h e J r e g i s t e r and enter t h e J loop. C l e a r F l , F2, F3 t o zero. I n i t i a l i z e t h e I r e g i s t e r and enter t h e I loop. Test and s k i p i f I = J . C a l c u l a t e RR - DX1**2 + DX2**2 + DX3**2 as a double word sum. Test and s k i p i f RR > RRC. (Separation beyond c u t o f f range). C a l c u l a t e R = SQRT(RR) and s t o r e . Form R-RS and set equal t o zero i f negative. S h i f t r i g h t 3 places ( d i v i d e by 8) and use as index t o look up F. Form the components of F: (F*DXl/R), (F*DX2/R), (F*DX3/R), and add t o F l , F2, F3. T


MINICOMPUTERS


142

AND LARGE SCALE

COMPUTATIONS

Increment I and continue. On e x i t from I loop c a l c u l a t e momentum changes as F l , F2, F3 s h i f t e d r i g h t k places ( d i v i d e d "by 16, corresponding t o At = l / l 6 ) and update P l ( j ) , P 2 ( j ) , and P5(J). C a l c u l a t e p o s i t i o n changes as P l ( j ) , P 2 ( j ) , and P 3 ( j ) , s h i f t e d r i g h t k p l a c e s , and update X l ( j ) , X 2 ( j ) , and X5(J). D i s p l a y new p o s i t i o n . Increment J and continue. On e x i t from J loop, store panel switches and t e s t f o r e x i t t o monitor, c o o l i n g , or temperature c a l c u l a t i o n . Increment Τ r e g i s t e r and continue. D i s p l a y and Output o f R e s u l t s . The c a l c u l a t i o n s produce two s o r t s o f r e s u l t s : p a r t i c l e p o s i t i o n s as a f u n c t i o n o f time, and s t a t i s t i c a l f u n c t i o n s o f t h e system. The p a r t i c l e p o s i t i o n s as a f u n c t i o n o f time are d i s p l a y e d on t h e storage CRT. F o r reasons o f c l a r i t y , only those p a r t i c l e s i n t h e f i r s t octant ( a l l components o f p o s i t i o n p o s i t i v e ) are d i s p l a y e d . A f t e r the new p o s i t i o n o f t h e J ' t h p a r t i c l e i s c a l c u l a t e d , t h e t h r e e components o f p o s i t i o n a r e t e s t e d , and, i f a l l are p o s i t i v e , t h e χ and y components o f p o s i t i o n are t r a n s m i t t e d through t h e CRU t o t h e 10 b i t analog t o d i g i t a l converters which d r i v e t h e CRT d i s p l a y u n i t . We thus d i s p l a y t h e p r o j e c t i o n on t h e x-y plane o f t h e content o f t h e f i r s t octant. P l a c i n g t h e d i s p l a y u n i t i n storage mode r e s u l t s i n the development o f t r a c e s o f t h e paths o f the centers o f t h e p a r t i c l e s . Some t y p i c a l t r a c e s a f t e r v a r y i n g lengths o f time (θ.36 ps, 0.9 ps, 1.8 p s ) a r e shown i n Figures 1, 2 and 3. Such d i s p l a y s are quite valuable f o r v i s u a l i s i n g t h e nature o f a l i q u i d ( i t r a p i d l y becomes obvious t h a t a l i q u i d i s n e i t h e r g a s - l i k e nor s o l i d - l i k e ) , but, obviously some q u a n t i t a t i v e c h a r a c t e r i s t i c s a r e needed. Two widely used s t a t i s t i c s o f a l i q u i d a r e t h e mean square displacement f u n c t i o n , and t h e v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n . We take t h e mean square displacement f u n c t i o n , χ ( t ) , as t h e ensemble average: 2

x ( t ) = i

i

1

It i s , o f course, equal t o t h e time average f o r any p a r t i c l e : 2

x ( t ) = x

T

but t h e f i r s t form i s more convenient here. To evaluate i t , we choose a s t a r t i n g time a f t e r t h e system has reached e q u i l i b r i u m , as determined by t h e constancy o f t h e "temperature". We take t h i s time as t = 0, and store t h e values o f XI, X2, and X3 f o r a l l t h e



11. F L i N N


143

Figure 1. Projection on the x-y plane of the tracer of the motion of the center of simulated argon atoms in thefirstoctant of the system. Temperature is 90 K; elapsed time 0.36 picoseconds; computation time, 2 minutes.

^ •

^

j ^

* \

φ

•η) ν \.

w

Figure 2. Projection of tracer of motion of argon atoms as in Figure 1, but after total elapsed time of 0.9 picoseconds; computation time, 5 minutes

Figure 3. Projection of tracer of motion of argon atoms as in Figures 1 and 2, but after total elapsed time of 1.8 picoseconds; com putation time, 10 minutes


MINICOMPUTERS

144

A N D LARGE SCALE

COMPUTATIONS

p a r t i c l e s as X S 1 , X S 2 , and X S 3 . At each time step, a f t e r completion of the J loop, t h e current value o f x ( t ) i s evaluated by summing ( X l ( l ) - X S l ( l ) ) * * 2 + ( X 2 ( l ) - X S 2 ( l ) ) * * 2 + ( X 3 ( l ) - X S 3 ( l ) ) * * 2 over a l l p a r t i c l e s and d i v i d i n g by N. The r e s u l t i n g f u n c t i o n can be d i s p l a y e d a t any time on the CRT or punched out on paper tape at the conclusion o f a run f o r p l o t t i n g on a pen and ink p l o t t e r . A t y p i c a l p l o t i s shown i n Figure k The normalized v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n i s c a l c u l a t e d i n a s i m i l a r way. I t i s defined as: 2


m

2

0 ( t ) = /. At the t = 0 chosen as d e s c r i b e d above, we store the values o f PI, P 2 , P3 f o r a l l the p a r t i c l e s as PSI, P S 2 , P S 3 . A f t e r each time step we form P1*PS1 + P 2 * P S 2 + P3*PS3 f o r each p a r t i c l e , sum over a l l p a r t i c l e s , and normalize by d i v i s i o n by the sum of P S 1 * * 2 + P S 2 * * 2 + P S 3 * * 2 f o r a l l p a r t i c l e s . This f u n c t i o n a l s o can be d i s p l a y e d on the CRT or punched out f o r e x t e r n a l p l o t t i n g . A t y p i c a l p l o t o f t h i s f u n c t i o n i s shown i n F i g u r e 5 · Comparison with Standard C a l c u l a t i o n s . The r e s u l t s obtained i n t h i s i n v e s t i g a t i o n are c o n s i s t e n t with those obtained i n conventional l a r g e machine c a l c u l a t i o n s , but the cost per computation i s very much lower. The q u a l i t a t i v e f e a t u r e s seen i n Figures 1-5 are the same as those reported f o r the standard c a l c u l a t i o n s ; a q u a n t i t a t i v e t e s t i s provided by the d i f f u s i o n c o e f f i c i e n t , which i s a s e n s i t i v e t e s t o f the technique. Levésque and V e r l e t ( 6 ) have summarized the r e s u l t s o f t h e i r c a l c u l a t i o n s f o r argon with the e m p i r i c a l formula i n reduced units: D = 0.006^23 T/p

2

+ 0 . 0 2 2 2 - 0 . 0 2 8 0 p.

Converted t o SI u n i t s , t h i s becomes: D = 5.639 X 1 0 "

5

T/p

2

+ 8.270 X 1 0 "

9

- 6.2C7 X 1 0 '

1 2

p.

For t h e c o n d i t i o n s corresponding t o the data shown i n Figure k, Τ = 113 Κ, and Ρ = l kh6 X 1 0 " 9 m /s, t h e i r equation p r e d i c t s D = 2 . 3 ^ X 1 0 " n r / s . The data o f F i g u r e 2 correspond t o a D - 2 . 7 7 Χ Ι Ο " m /s. T h i s d i f f e r e n c e i s of the same order as t h e s c a t t e r of standard c a l c u l a t i o n s , and the discrepancy between c a l c u l a t i o n and experiment. The speed o f the c a l c u l a t i o n was quite reasonable: 282 time steps i n 1 0 minutes, or 1692 steps per hour. Each second of machine time corresponds t o 3 X 1 0 " l 5 seconds i n argon. F o r comparison, the reported r a t e achieved on a l a r g e machine, a CDC 6 6 0 0 , was I5OO steps per hour f o r a somewhat l a r g e r system (Q6k p a r t i c l e s ) ( 2 ) . With proper programming, the c a l c u l a t i o n 2

m

9

9

2



11. F L i N N


Time

145

(Picoseconds)

Figure 4. Mean square displacement of simuhted argon atoms at 113 Κ as a function of time

I 0.0

ι 0.2

ι

ι

ι

ι

0.4 0.5 0.8 1.0 Time (Picoseconds)

ι

ι

1.2

1.4

Figure 5. Normalized velocity autocorrelation function of simulated argon atoms at 113 Κ


MINICOMPUTERS


146

AND LARGE SCALE

COMPUTATIONS

time v a r i e s approximately as N. With t h i s allowance, i t appears t h a t c a l c u l a t i o n on a minicomputer i s slower "by roughly a f a c t o r of 6 than on a l a r g e machine. To estimate t h e r e l a t i v e cost o f computation, we take t h e i n i t i a l cost of t h e system and d i s t r i b u t e i t over t h r e e years (a conservative procedure, s i n c e our machine has been i n continuous use f o r f i v e years with no maintenance contract and n e g l i b l e s e r v i c i n g ) . T h i s corresponds t o a cost of $6Λθ a day or $0.27 per hour. I f we assume a l a r g e machine cost of about $100 per hour, t h e cost of equivalent c a l c u l a t i o n s i s lower on t h e small machine by a f a c t o r o f about 50. Future Prospects. The use o f c u r r e n t l y a v a i l a b l e hardware, i n s t e a d of t h e obsolete 96ΟΑ, would make p o s s i b l e both g r e a t e r savings and much more ambitious c a l c u l a t i o n s . In p a r t i c u l a r , a t h r e e dimensional a r r a y o f microprocessors (such as t h e T I 9900), each assigned a p o r t i o n o f the volume under study, could be used t o i n c r e a s e t h e speed o f c a l c u l a t i o n by more than an order of magnitude f o r a c o s t i n c r e a s e o f about a f a c t o r o f two.

Literature Cited. (1) Rahman, Α., Phys. Rev. (1964), 136, A405. (2) McDonald, I. R. and Singer, Κ., Quart. Rev. (1970), 24, 238. (3) Rahman, A. in "Interatomic Potentials and Simulation of Lattice Defects", ed. by Gehlen, P. C., Beeler, J . R. J r . , and Jaffee, R. I., p. 233, Plenum, N.Y., 1972. (4) Fisher, R. A. and Watts, R. O., Aust. J . Phys. (1972), 25, 529. (5) Flinn, P. Α., in "Mössbauer Effect Methodology", vol. 9, ed. by Gruverman, I. J., Seidel, C. W., and Dieterly, p. 245, Plenum, N.Y. 1974. (6) Levesque, D. and Verlet, L . , Phys. Rev. (1971), A2, 2514.


Molecular Dynamics Calculations on a Minicomputer

Recommend Documents