A sample of water, consisting of 216 rigid molecules at mass density 1 gm/cm³, has been simulated by computer using the molecular dynamics technique. The system evolves in time by the laws of classical dynamics, subject to an effective pair potential that incorporates the principal structural effects of many-body interactions in real water. Both static structural properties and the kinetic behavior have been examined in considerable detail for a dynamics “run” at nominal temperature 34.3°C. In those few cases where direct comparisons with experiment can be made, agreement is moderately good; a simple energy resealing of the potential (using the factor 1.06) however improves the closeness of agreement considerably. A sequence of stereoscopic pictures of the system’s intermediate configurations reinforces conclusions inferred from the various “run” averages: (a) The liquid structure consists of a highly strained random hydrogen-bond network which bears little structural resemblance to known aqueous crystals; (b) the diffusion process proceeds continuously by cooperative interaction of neighbors, rather than through a sequence of discrete hops between positions of temporary residence. A preliminary assessment of temperature variations confirms the ability of this dynamical model to represent liquid water realistically.

I. Introduction

Although water occupies a preeminent position among liquids, this substance has not enjoyed the attention of a rapidly developing body of statistical mechanical theory devoted specifically to its own properties.¹ One obvious reason for this retardation is the internal structure of the water molecule, which at the very least requires considering orientational degrees of freedom. In addition, the potentials of interaction for water molecules have until very recently^2-5 been imperfectly known. Furthermore, it now appears that these interactions are nonadditive to a significant degree.^6-9 Finally, the maximum cohesive binding between pairs of molecules, in units of k_BT at the triple point, is roughly an order of magnitude greater than the same quantity for the theoretically popular liquified noble gases. This combination of complications renders impractical a large part of conventional liquid state theory for studying water. One must forego reliance on the integral equation approaches to static pair correlation functions on the one hand, while it is clear on the other hand that the fundamental theory of kinetic processes becomes even more complex than usual. Under these circumstances, the most promising approach at present seems to be the direct simulation of liquid water by electronic computer. Both the “Monte Carlo” method¹⁰ and the technique of "molecular dynamics"¹¹ are available for this purpose. The former offers the possibility of generating canonical ensembles of given temperature, but it is entirely restricted to a study of static structural properties. Molecular dynamics (which is nominally microcanonical) however can probe both static and kinetic behavior, so in principle it is the more powerful tool. We have therefore chosen to utilize this more powerful approach. This paper provides details of computational strategy, and initial results, in our molecular dynamics investigation of liquid water. The following Sec. II specifies the Hamiltonian used for our dynamical water model. The individual molecules are treated as rigid asymmetric rotors, i.e., their internal bond lengths and bond angles are invariant. Classical mechanics describes the temporal evolution of our model system. The coupled differential equations for translational and rotational motion are considered in Sec. III for the model water system. Special choices are introduced there for system size (216 molecules), boundary conditions (periodic unit cell corresponding to liquid at 1 gm/cm³), and time increment for numerical integration of the coupled dynamical equations (delta_t=4.355 X 10^-l6 sec). Discussed as well in Sec. III is a force truncation scheme. Section IV presents a body of results thus far accumulated which specifically bears on the static molecular structure for our water model. Separate radial correlation functions are reported for the three distinct types of nuclear pairs present (O-O, O-H, and H-H), and these are used to synthesize the hypothetical x-ray scattering pattern for the model liquid. Several aspects of the elaborate local orientational order are presented in Sec. IV by examining dipole direction correlation in successive concentric shells about a given molecule and by analyzing the oxygen-oxygen pairs in separate icosahedral sectors about a fixed molecule. The character of hydrogen bonding in the liquid has also been examined using the distribution function for pair interaction energy. Having thus described the main features of equilibrium molecular order, we pass on to kinetic properties of the water model in Sec. V. Several autocorrelation functions are presented which reveal distinctive characteristics of translational and rotational motion. These autocorrelation functions permit one in principle to calculate the self-diffusion constant, the dielectric relaxation spectrum, neutron inelastic scattering, and NMR spin-lattice relaxation. In order to supplement these conventional molecular dynamics quantities, we have also produced stereoscopic photographs (of a cathode-ray display) which visually present instantaneous configurations of the 216 molecules during the system’s temporal evolution. Unfortunately, a printed paper such as this one does not provide an effective direct way to communicate these elaborate stereoscopic pictures. However, we have attempted verbally to summarize their contribution to our own understanding of liquid water at the appropriate points in Secs. IV and V. In particular, these pictures allow one to perceive the global features of liquid water hydrogen-bond patterns, and to appreciate details of local cooperativity in molecular Brownian motion. The results in Secs. IV and V refer to a single computer “run” corresponding to water at a fixed temperature. Some early results for a substantially lower temperature are mentioned in Sec. VI. Although we reserve most of the details concerning temperature variations for a later publication, these few observations strengthen our conviction that the model used is a relatively faithful representation of real water. Several items are taken up for discussion in the final Sec. VII. We list there some extensions of the present project that appear to us to have relatively high scientific merit. Included among these are possible modifications of the Hamiltonian that could well be required at the next precision level of computer simulation for aqueous fluids. In particular, we stress a simple energy rescaling that may be applied to the present results, which seems to produce substantial improvement in agreement with experiment.