12878
J. Phys. Chem. 1996, 100, 12878-12887
Ab Initio Molecular Dynamics Simulations Mark E. Tuckerman, P. Jeffrey Ungar, Tycho von Rosenvinge, and Michael L. Klein* Department of Chemistry, UniVersity of PennsylVania, Philadelphia, PennsylVania 19104-6323 ReceiVed: February 16, 1996X
Over the past decade, new simulation methodologies, such as the Car-Parrinello ab initio molecular dynamics technique, have become increasingly important as tools to study and characterize condensed phase molecular systems. We emphasize the versatility of these new approaches to simulation by reviewing selected applications to molecular crystals, liquids, and clusters, which highlight a range of interesting phenomena. The molecular crystals white phosphorus, nitromethane, and hydrogen chloride dihydrate exhibit molecular reorientation phenomena, methyl torsional motion, and proton-hopping events, respectively. We indicate how, in the latter examples, it is now possible to include quantum effects in the simulation of the proton motion. Ionic solvation and proton transport in water are used to illustrate the current status of simulations of liquid systems. The final topic in our survey deals with the possibility of including the quantum nature of nuclear motions into the simulation methodology of clusters.
I. Introduction Classical molecular dynamics (MD), which has its origins in the 1950s and 1960s, is now widely used to investigate condensed phase systems, ranging from micellar solutions to biomolecules, polymers, and inorganic materials. The standard MD technique, which is based on the use of empirical interatomic potential functions parametrized to experimental data, has also proved its worth in the study of ionic solutions, electrochemistry, tribology, and cluster science. One of the most exciting developments of the past decade is the advent of ab initio molecular dynamics, which, rather than requiring an empirical interaction potential as input, utilizes interatomic forces computed directly from the electronic structure. This new approach has now become established as an important tool in the investigation of chemically complex environments. In the Car-Parrinello (CP) approach1,2 to ab initio MD, the electronic structure is described using the Kohn-Sham formulation3 of the density functional theory,4 and the Kohn-Sham orbitals are expanded in a plane wave basis. The expansion coefficients are treated as a set of fictitious dynamical variables that are propagated adiabatically with respect to the nuclei, so that, at each time step, they describe the instantaneous ground state Born-Oppenheimer surface. In this way, the need to solve the Kohn-Sham equations explicitly is avoided. The CP approach has proved useful in a wide variety of physical and chemical applications. The present article illustrates the versatility of the CP approach to MD simulations by considering a variety of condensed phase systems. The first examples presented concern the simulation of molecular crystals. We begin with the low-symmetry, β-P4 phase crystal of white phosphorus. The common R-phase transforms to the β-phase below 196 K, at atmospheric pressure. The structure of this phase of phosphorus been determined5 using X-ray diffraction and verified6 using neutron scattering. The unit cell is triclinic and contains six somewhat distorted tetrahedral P4 molecules with inversion symmetry only (space group P1h). The CPMD scheme allows the study of both the structure of this low-order phase and the reorientational dynamics of the molecules, which are subject to weak intermolecular forces. Since several aspects of the crystal dynamics are X
Abstract published in AdVance ACS Abstracts, July 1, 1996.
S0022-3654(96)00480-7 CCC: $12.00
accessible from the neutron experiments, direct comparisons with the simulation results7 on β-P4 are possible. As an example of a crystal system with a torsional degree of freedom, we consider the nitromethane crystal. Nitromethane molecules contain a methyl group weakly bonded to an NO2 group. In the gas phase, the barrier against internal methyl rotation about this bond is very small (∼6 cal/mol), which raises interesting questions concerning the nature of this motion in the crystalline state. Since intermolecular interactions determine the dynamics of the crystal, there is the intriguing possibility of exploring the coupling between the rotational modes of each molecule and cooperative phenomena. Another aspect of molecular crystals that can be explored using the CP scheme is the displacement of protons through hydrogen bonds.8 The study of such processes allows the investigation of the role of structural relaxation involving other nuclei and the assessment of the contribution of tunneling. As an example, we consider the hydrogen chloride dihydrate crystal.9 Because of its low mass, the proton often exhibits quantum mechanical behavior, and so the usual approximation of classical nuclear motion is often too crude when applied to proton transfer processes in molecular crystals. In order to incorporate quantum effects into the CP scheme, we must go beyond this approximation. As yet, there is no satisfactory method for obtaining true quantum dynamics at finite temperature, but it is possible to obtain quantum mechanical averages of equilibrium properties using the path integral approach.10-14 In the path integral formulation, the quantum canonical partition function, Tr(e-βH), is expressed as a sum over paths in imaginary time that are periodic with a period equal to β ) 1/kBT. In practice, the path integral is represented as a multidimensional integral over discretized paths, a form that is equivalent to a closed, classical polymer chain with nearestneighbor harmonic coupling.15 Such an integral can be evaluated using MD techniques and therefore can be incorporated into the CP scheme. However, because of the stiffness of the harmonic coupling, a straightforward MD approach will suffer from the fact that the trajectories will not be ergodic and hence not sample the entire configuration space.16 However, a MD algorithm for path integrals exists11 that is almost as efficient as more traditional Monte Carlo schemes. The algorithm involves a change of variables that diagonalizes the harmonic © 1996 American Chemical Society
Molecular Dynamics Simulations interaction and a coupling of each degree of freedom to a Nose´Hoover chain thermostat17 to ensure ergodic motion of the new variables. An important application of the CPMD simulation technique presented here concerns the liquid state, the solvation of ions, and the transport of protons in ionic solutions. This work is particularly compelling because it gives insight into the mechanism behind the anomalously high mobility of protons in solution. In order to gain physical insight into this process, it is essential to have access to the microscopic motion of individual atoms during proton transfer events. Such information is often not obtainable experimentally but is given directly in a simulation. Since empirical potentials, which describe the interaction of an hydronium ion with a solvent like water, are difficult to construct, the ab initio approach is an attractive alternative. Here, we will review and discuss recent work on the solvation of H+ and OH- in liquid water.18 The ability not only to obtain solvation structures of these ions but also to shed light on the physical mechanism and rate-limiting steps of proton transport in these systems is among the significant results of this work. In the above-mentioned applications of CPMD, quantum effects of the proton are minimized by using deuterium instead of hydrogen. In order to assess the importance of quantum effects on solvation and proton transport, a fully quantum mechanical study of H+ and OH- in water should be carried out. However, since these calculations are computationally intensive, researchers have begun by considering small charged water clusters, where long runs and good statistics can be collected with a more reasonable expenditure of CPU of time.19 Finally, we mention possible future applications of ab initio MD techniques in biological systems. At present, ab initio MD methodology has been applied to small biological molecules20 and has been used in an ambitious attempt to elucidate the mechanism of proton transport in an ion channel containing a strand of hydrogen-bonded water molecules.21 The results of these studies are sufficiently encouraging to give optimism that as better algorithms are developed and enhanced computing power becomes available, even more complex biological systems can be investigated with the CPMD approach. Although we do not provide a detailed discussion in this article, researchers have also applied MD to systems involving other complications. For example, Silvestrelli et al.86 have studied molten Kx(KCl)1-x using a molecular dynamics technique based on finite temperature DFT in order to incorporate hot electrons. The CPMD method has also been extended to allow the simulation of systems at constant pressure. Recent applications discussed in the literature include the formation of hydrogenated amorphous carbon from the compression of polyacetylene87 and the prediction of a new high-pressure phase of ice.58 Finally, it is occasionally possible to obtain MD trajectories on an excited state potential energy surface, such as in the work of Bernasconi et al.88 on the solid state polymerization of acetylene under pressure. In this case, a triplet exciton self-trapped on a single, cis-bent molecule is incorporated using the local spin density approximation. The paper is organized as follows: In section IIA, the CP Lagrangian and equations of motion are reviewed. The discussion is generalized to allow the use of ultrasoft pseudopotentials22 for treatment of core electrons and involves the use of the recently introduced constraint nonorthogonal orbital (CNO) method.23 In section IIB, the path integral CP method is discussed, and the MD algorithm11 is described and compared to more straightforward approaches. Applications of these methods to the molecular crystals mentioned above are presented
J. Phys. Chem., Vol. 100, No. 31, 1996 12879 in subsequent sections. A review of recent work on solvation and charge transport is then given, and future applications are discussed. II. Simulation Methodologies A. Ab Initio Molecular Dynamics. The equations of motion for CPMD are derived from an extended Lagrangian in which electronic orbitals {φi} execute fictitious dynamics that allows them to follow the motion of nuclei having positions {RI} while remaining very close to the instantaneous ground state BornOppenheimer surface. This is accomplished by keeping the fictitious kinetic energy of the orbitals small compared to that of the nuclei and assigning to them a small fictitious mass parameter µ. The interatomic forces are computed via the Hellman-Feynman theorem from the instantaneous electronic structure at each nuclear configuration. The extended Lagrangian takes the form
L ) µ∑〈φ˙ i|φ˙ i〉 + 1/2∑MIR4 I2 - E[{φi},{RI}] + i
I
Λij(〈φi|Sˆ ({RI})|φj〉 - δij) ∑ i,j
(1)
where {MI} are the masses of the nuclei, and E[{φi},{RI}] is the Kohn-Sham energy density functional given by
E[{φi},{RI}] ) -
1
1
〈φi|∇2|φi〉 + ∫dr dr′ ∑ 2 i 2
n(r) n(r′) |r - r′|
+
Exc[n] + ∫dr Vloc(r) n(r) + ∑〈φi|Vˆ NL|φi〉 + U({RI}) (2) i
In eq 2, the external potential is expressed in terms of an atomic pseudopotential that contains a local component Vloc(r) and a nonlocal operator Vˆ NL. The use of pseudopotentials eliminates22,24 the need to treat core electrons explicitly. The term “ion” is used to refer to the unit consisting of the nucleus with core electrons attached. n(r) is the electronic density, Exc[n] is the exchange and correlation functional, and U({RI}) is the ionion interaction energy. The last term in eq 1 is not technically part of the Lagrangian, but we have placed it here for convenience of notation. A problem arises because Lagrange’s equations of motion assume that independent variations can be made in the orbitals {φi}, but they are actually subject to a set of constraints
〈φi|Sˆ ({RI})|φj〉 ) δij
(3)
These are incorporated in the last term of eq 1 using the method of Lagrange multipliers.25 The multipliers Λij are chosen to ensure the fulfillment of the constraints while allowing the equations of motion to be obtained from Hamilton’s variational principle. Note that eq 3 takes the form of a generalized orthonormality condition.22 The Hermitian overlap operator Sˆ is given by Sˆ ({RI}) ) ˆI + Pˆ ({RI}). The imposition of this constraint condition allows the use of ultrasoft pseudopotentials with eq 1, which also require augmentation of the electronic charge density n(r).22 For norm-conserving pseudopotentials, Pˆ is zero and Sˆ reduces to an identity operator, which gives the usual orthonormality constraint. For ultrasoft pseudopotentials, Pˆ is composed of projectors onto functions that characterize the pseudopotential. Details of the implementation of ultrasoft pseudopotentials in CPMD are given in the literature.26 Because the constraints couple both electronic and ionic degrees of freedom (through the dependence of Sˆ on RI), the usual SHAKE and RATTLE27,28 procedures used to enforce
12880 J. Phys. Chem., Vol. 100, No. 31, 1996
Tuckerman et al.
them must be iterated through the ionic positions. This complication can be avoided by introducing a transformation of the electronic orbitals from the set {φi} to a new set {ψi} given by |φi〉 ) ∑j|ψj〉Tji. The transformation matrix T ) O-1/2, where Oij ) 〈ψi|Sˆ ({RI})|ψj〉 is the matrix element of the operator Sˆ with respect to the new orbitals.23,29,30 Both O and therefore T are Hermitian. Thus, given the definition of the orbital transformation, the constraints in eq 3 are satisfied by construction. The orbitals {ψi} are used instead of {φi} to rewrite the CP Lagrangian. Since eq 3 is satisfied implicitly, no explicit constraint on the new orbitals is required. However, for use in a MD scheme, an arbitrary set of constaint conditions, σR[{ψi}] ) 0, is imposed on the orbitals {ψi} in order to prevent them from becoming linearly dependent. Incorporating the new orbitals and ad hoc constraints into the CP Lagrangian gives
L ) µ∑〈ψ˙ i|ψ˙ i〉 + 1/2∑MIR4 I2 - ∑Mji〈ψi|H ˆ |ψj〉 + i
i,j
I
∑R ΛRσR[{ψi}]
(4)
∂〈ψi| (5)
where (∆I)ij ) 〈φi|∇ISˆ |φj〉. FI is the force that would be derived from eq 1, and Hij is the matrix element of the Kohn-Sham Hamiltonian with respect to the orthogonal orbitals Hij ) ˆ |φj〉. The transformed electronic force |φ˜ i〉 is given by 〈φi|H
ˆ |φj〉 - ∑Sˆ |φk〉〈φk|H ˆ |φj〉)Tji |φ˜ i〉 ) - ∑(H j
E0({RI(τ)})]} (7) where E0({RI}) is the ground state energy surface (see ref 13 for details). Equation 7 can be written as the limit of a discretized “classical phase space” integral:10
{ [∑ (
Pf∞
I)1
P
∂σR
2 I ) FI + Tr(HδI) MIR
I)1
N
An where H ˆ is the Kohn-Sham Hamiltonian and M ) alternative formulation of the energy density functional in terms of nonorthogonal orbitals has been presented by Mauri and Galli.31 While no explicit matrix inversion is required in their scheme, departures from the exact energy density functional can occur when the orbitals are not at the precise minimum, making their approach less suitable for MD. The equations of motion derived from eq 4 are
R
N
Z ) ∫D[R1]...D[RN] exp{- ∫0 dτ [∑1/2MIR4 I2(τ) + β
Z ) lim [∏N ∫dRI(1)...dRI(P)∫dPI(1)...dPI(P)] ×
O-1.
˜ i〉 + ∑ΛR µ|ψ ¨ i〉 ) |φ
can be seen by writing down the imaginary time path integral expression for the quantum canonical partition function in the Born-Oppenheimer limit.32 With the approximation that the temperature kBT ) 1/β is small compared to the gap in the electronic spectrum, the path integral expression for the partition function takes the form
(6)
k
This is just the usual electronic force rotated by the matrix Tji and with Λij replaced by 〈φi|H ˆ |φj〉. While the constrained nonorthogonal orbital (CNO) method requires that some condition be imposed, a simpler condition than eq 3 can be chosen. In particular, the simple orthonormality constraint, 〈ψi|ψj〉 ) δij, and the norm constraint, 〈ψi|ψi〉 ) 1, conditions have been investigated and shown to be satisfactory choices.23 For norm constraints, it is necessary to monitor the off-diagonal elements of the overlap matrix 〈ψi|ψj〉. If they are observed to grow beyond a specified tolerance, the orbitals must be rotated using the unitary transformation that diagonalizes this matrix. This operation will set the off-diagonal elements once again to zero but will not effect the dynamics, since the CP Lagrangian is invariant under unitary transformation in the space of occupied states. A tolerance of 10-3 has been found to be sufficiently small to maintain a stable scheme.23 The norm constraint has the advantage that it constitutes only N conditions as opposed to the N2 conditions required by the full orthonormality constraint. B. Ab Initio Path Integral Molecular Dynamics. Quantum corrections to the equilibrium properties of the ionic subsystem can be computed by applying the CP method to path integral (PI) molecular dynamics.12-14 The validity of this approach
exp - β∑
N
s)1 I)1
(PI(s))2
)
1 + MIωP2(RI(s) - RI(s+1))2 + 2M′I 2 1
]}
E0({RI}(s))
P
(8)
where {RI}(s) denotes the set of all nuclear positions at imaginary time slice s, ωP is the chain frequency xP/β, {M′I} is a set of fictitious bead masses, and N is an overall normalization constant. The path periodicity condition RI(P+1) ) RI(1) is assumed. Note that the Gaussian integral over the fictitious momenta simply yields a temperature-dependent prefactor. From eq 8, it can be seen that the ground state energy surface must be evaluated at each imaginary time slice, a fact requiring the introduction of a set of Kohn-Sham orbitals for each time slice.12 Equation 8 can be evaluated for finite (but large) P using molecular dynamics techniques. An extended Lagrangian corresponding to eq 8 can be constructed12-14 that leads to the following set of equations of motion:
∂σR ˜ i(s)〉 + ∑Λ(s) µ|ψ ¨ i(s)〉 ) |φ R R ∂〈ψi(s)| 2 I(s) ) -MIωP2(2RI(s) - RI(s+1) - RI(s-1)) MIR
˜ (s) 1 ∂E P ∂R(s)
(9)
I
where E ˜ is the energy density functional expressed in terms of the nonorthogonal orbitals: E ˜ (s) ) ∑i,j Mji(s)〈ψ(s) ˆ (s)|ψ(s) i |H j 〉. Equation 9 can, in principle, generate the quantum canonical distribution. However, due to the presence of stiff harmonic forces arising from the quantum kinetic energy, the sampling efficiency will be very poor.16 An efficient PIMD algorithm exists11 that employs a change of variables to the staging modes used in path integral Monte Carlo calculations.33-36 The use of staging coordinates has the effect of diagonalizing the quantum kinetic energy. In addition, each degree of freedom is coupled to a Nose´-Hoover chain thermostat17 to ensure that the motion of each mode is ergodic. This method has recently been incorporated into the CP scheme.14 The staging transformation is defined by
uI(1) ) RI(1),
*
uI(s) ) RI(s) - RI(s) (s ) 2, ..., P) (10)
Molecular Dynamics Simulations
J. Phys. Chem., Vol. 100, No. 31, 1996 12881
*
where RI(s) ) [(s - 1)RI(s+1) + RI(1)]/s. The inverse staging transformation can be expressed as a recursion:
RI(1) ) uI(1),
RI(s) ) uI(s) +
s - 1 (s+1) 1 (1) + RI RI s s (s ) 2, ..., P) (11)
The s ) P term is used to start the recursion using the path periodicity condition RI(P+1) ) RI(1). A more general form of the staging transformation that is useful for very large P is given in ref 11. The staging transformation is substituted into eq 8 as a change of variables, giving the following expression for the partition function: N
Z ) lim [∏N ′ ∫duI(1)...duI(P) ∫dPI(1)...dPI(P)] × Pf∞
{ [(
)
I)1
P
exp -β∑
N
∑ s)1 I)1
(PI(s))2
1 + MI(s)ωP2 (uI(s))2 + 2M′(s) 2 I
1
the electronic subsystem. The thermostat forces are given by (s) (s) 2 GI,l,1 ) M′I(s)(u˘ I,l ) -
]}
E0({RI(uI(1),...,uI(P))}(s))
P
[
(s) (s) ) QR(ξ˙ I,l,M-2 )2 GI,l,M-1
(12)
where the masses MI(s) replace the physical masses MI and are referred to as the “staging masses”. They are given by MI(1) ) MI and MI(s) ) sMI/(s - 1) for s ) 2, ..., P. The fictitious masses M′I(s) are assumed to be a constant multiple of the staging masses, i.e., M′I(s) ) cMI(s). The notation {RI (uI(1),...,uI(P))}(s) is intended to indicate that, at each time slice s, the elements of the set {RI}(s) of Cartesian bead positions for the different atoms are functions of all the staging variables uI(1),...,uI(P). In order to produce the distribution in eq 12, the following set of equations of motion is employed:
∂σR - µη˘ (s) ˙ i(s)〉 ˜ i(s)〉 + ∑Λ(s) µ|ψ ¨ i(s)〉 ) |φ R 1 |ψ (s) R ∂〈ψi | (s) (s) M′I(s)u¨ I,l ) -MI(s)ωP2 uI,l -
Figure 1. Cumulative average of the primitive energy estimator (au) in eq 16 for a path integral simulation of a single H2O molecule using various approaches: eqs 9 (method I), ref 13 (method II), and the methods discussed in the text (method III).
(s) (s) GI,l,ν ) QR(ξ˙ I,l,ν-1 )2 -
)
(s) GI,l,M
[
-
(s) Qκη¨ κ(s) ) γκ(s) - Qκη˘ κ(s)η˘ κ+1
κ ) 1, ..., Melec - 1
(s) (s) (s) (s) ) γM - QMη˘ M-1 η˘ M QMη¨ M
1 βe
(s) 2 γκ(s) ) Qκ-1(η˘ κ-1 ) -
(13)
In eqs 13, the index l ) x, y, z runs over the Cartesian components of vectors uI(s). Both the electronic and the ionic subsystems have been coupled to independent sets of chain thermostats of lengths Melec and Mion, respectively. For the ions, each Cartesian component of each mode variable is coupled to its own Nose´-Hoover chain thermostat with the thermostat mass parameter QR ) 1/(βωP2). For the electrons, a single thermostat chain is used with masses Q1 ) 2Ee/ωe2 and Qκ ) 2Ee/(Neωe2), where Ee is the desired electron kinetic energy, ωe is a characteristic frequency of the fictitious electron dynamics, and Ne is the number of fictitious dynamical degrees of freedom in
] [
]
1 1 (s) 2 + QM(η˘ M ) βe βe
κ ) 2, ..., Melec - 2, Melec (14)
Using the recursive inverse transformation, eq 11, the forces on the staging modes can be computed straightforwardly in (s) terms of the primitive forces - ∂E ˜ (s)/∂RI,l
P ∂u(1)
(s) (s) QRξ˙ I,l,M-1 ξ˙ I,l,M
]
ν ) 2, ..., Mion - 2, Mion
(s) (s) γM-1 ) QM-2(η˘ M-2 )2 -
[
I,l
)
1
P
∑
∂E ˜ (s)
Pτ)1∂R(τ) I,l
]
∂E ˜ (s) ˜ (s) 1 (s - 2) ∂E ˜ (s) 1 ∂E ) + P ∂u(s) P (s - 1) ∂u(s-1) ∂R(s) I,l I,l I,l
I,l
(s) QRξ¨ I,l,M
] [
i
˜ (s) 1 ∂E
ν ) 1, ..., Mion - 1
1 1 (s) 2 + QR(ξ˙ I,l,M ) β β
γ(s) ˙ i(s)|ψ˙ i(s)〉 - Ee 1 ) ν∑〈ψ
˜ 1 ∂E (s) (s) - M′I(s)ξ˙ I,l,1 u˘ I,l P ∂u(s)
(s) (s) (s) (s) QRξ¨ I,l,ν ) GI,l,ν - QRξ˙ I,l,ν ξ˙ I,l,ν+1
1 β
1 β
(15)
Although eqs 13 and 15 have been written for the staging modes, they could have equally well been written in terms of the normal modes of the chain.37 A comparison between staging and normal modes was carried out in ref 14, and the two approaches were found to exhibit equal sampling efficiency. These methods were also compared to the primitive scheme given in eqs 9 and to the method proposed in ref 13. The results of this comparison are shown in Figure 1, where the convergence of the primitive energy estimator
Eprim )
P
N
1 - ∑ ∑ MIωP2 (RI(s) - RI(s+1))2 + 2β s)1 I)1 2 1 P E[{ψi}(s),{RI}(s)] (16) ∑ Ps)1
3NP
for the four methods is plotted. The system studied consisted of a single water molecule in a box of length 5.3 Å with periodic boundary conditions at a temperature of 100 K with P ) 32.
12882 J. Phys. Chem., Vol. 100, No. 31, 1996
Tuckerman et al.
Further details of this comparison can be found in ref 14. In Figure 1 the cumulative average of the estimator in eq 16 over a runs of length 10 000 steps is shown. It can be seen from the figure that the method described above with either staging or normal modes converged significantly faster than the other two approaches. Since PICP calculations are extremely demanding computationally due to the requirement of P sets of KohnSham orbitals, the increased sampling efficiency afforded by the method described above can result in a substantial savings in CPU time. III. White Phosphorus at Low Temperatures Both the R (room temperature) and β (low temperature) phases of white phosphorus consist of tetrahedral P4 molecules, but until the past 10 years there has been little structural data on either phase. The room temperature form is thought to be body-centered cubic (lattice parameter a ) 18.51 ( 0.03 Å) with 56 or more molecules per cell (space group I4h3m-Td3).38 The low-temperature phase is now known to have a triclinic unit cell (at 50 K, a ) 11.45 Å, b ) 5.503 Å, c ) 11.26 Å, R ) 71.84°, β ) 90.37°, and γ ) 71.56°) containing six P4 molecules.5,6 The neutron scattering experiments tend to confirm earlier speculation39 that R-P4 is a rotator phase while β-P4 is an orientationally ordered crystal. We performed ab initio CPMD simulations of β-P4 at 50 and 145 K using the coordinates obtained from the reduction of the neutron data as a starting point. A system consisting of two unit cells along the b axis was used to even out the overall simulation cell dimensions and to avoid unduly constraining the phonon modes along this direction. The standard approach40,41 was used for the local density approximation (LDA) to the exchange-correlation functional. Since the LDA has been found to fail badly in the case of weak interactions,42,43 we have added gradient corrections (GC) to the functional, which have been shown to give a much better description of weakly bonded systems.43,44 They are taken from Becke45 for the exchange and Perdew46 for the correlation. The phosphorus core electrons were treated using a generalized norm-conserving pseudopotential.24,47 The electronic orbitals were expanded in plane waves using a cutoff of 12 Ry. The plane wave coefficients were given a mass µ ) 850 au, and the time step was set to δt ) 7 au. Constant volume trajectories spanning about 10 ps were generated for the two temperatures. Before proceeding with the MD calculations, the crystal structure was allowed to relax fully. This resulted in very little change except that the (∼3%) distortions from ideal tetrahedral molecules, which were found experimentally, mostly disappeared. Subsequent analysis of the MD trajectories showed that these distortions reappeared at finite temperature, suggesting that they are a dynamical effect and not an artifact of the reduction of the neutron scattering data. Figure 2 shows the time-averaged pair distribution function for the phosphorus atoms. The inset shows satellite peaks near the base of the main peak, which corresponds to the typical P-P bond length of 2.23 Å. The main peak contains a little more than two of the three nearest neighbors. Therefore, the remaining atom is found, on average, at a distance about -2% to +4% different from the other two. The frequencies present in the system are obtained by calculating the power spectral density of the atomic velocities, which may be obtained either from the autocorrelation functions 〈v(t)·v(0)〉 or equivalently from summing the spectral densities of all the velocity components. The resulting spectrum shows sharp peaks at 337, 434, and 575 cm-1 from the internal molecular vibrations. These are red-shifted by about 30 cm-1 from published frequencies for the isolated P4 molecule.48
Figure 2. Pair distribution functions for phosphorus atoms taken from simulations of white phosphorus at 50 and 145 K (offset for clarity). The inset, which shows the intramolecular peak near 2.2 Å in enhanced detail, reveals the distortion of the phosphorus tetrahedral.
Figure 3. Power spectral densities of the molecular center-of-mass velocities (a) and (b) angular velocities taken from simulations of white phosphorus at 50 and 145 K.
The modes due to molecular translation and rotation may be separately analyzed using the power spectral densities of the molecular center-of-mass translational and angular velocities. (The angular velocity of a vibrating molecule can be calculated from its instantaneous angular momentum and inertia tensor about its center of mass.) Figure 3 shows the resulting spectra for (a) translation and (b) rotation. Both span a range from low frequency to 100 cm-1, which implies a correlation time for the velocities of ∼0.33 ps. There is little difference between the angular velocity spectra for the two temperatures. If a transition to a rotator phase had occurred, we would expect a significant increase in the spectral density at low frequencies. The rotational motion is therefore mostly libration with a center frequency in the vicinity of 40 cm-1, which is in fair agreement with inelastic neutron scattering measurements below the β f R transition temperature.6 These show peaks centered near 27 cm-1 at 50 K and 20 cm-1 at 180 K, extending to 80-100 cm-1. Casual inspection of snapshots of the system at times 4-8 ps apart shows that the molecules do succeed in rotating completely around, but it is not clear whether or not the molecules are orientationally disordered. Further insight into the reorientational dynamics may be gained through the use of the tetrahedral rotor functions Mn(Ω), n ) 1-7, which are suitable normalized combinations of the seven l ) 3 spherical
Molecular Dynamics Simulations
J. Phys. Chem., Vol. 100, No. 31, 1996 12883
Figure 4. Mean reorientational correlation functions from CPMD simulation of white phosphorus at 50 and 145 K. The curves shown are averages over the normalized autocorrelation functions of the seven 7 tetrahedral rotor functions; i.e., R(t) ) ∑n)1 Rn(t).
harmonics. In Cartesian coordinates, they are
M1 ) (27/16)1/2∑xyz
(17)
M2 ) (9/320)1/2∑(5x2 - 3xr2)
(18)
M5 ) (27/64)1/2∑x(y2 - z2)
(19)
where r ) (x,y,z) is a unit vector from the center of a molecule to one of its atoms. The rotor functions not listed, n ) 3, 4 and n ) 6, 7, are obtained by cyclic permutation of the (x,y,z) variables in the M2 and M5 functions, respectively. The sum is to be taken over all for atoms for a given molecule, and the 7 Mn2 ) 1.49 normalization is fixed as ∑n)1 The rotor functions can serve as rotational order parameters for the system. For example, if the phase is orientationally disordered but only Td-type orientations are present (with respect to the selected x, y, and z axes), then 〈M1〉 ) 0 but 〈M12〉 ) 1. Our simulations show that β-P4 has persistent orientational order at both temperatures studied, since 〈Mn〉 * 0. It happens that the expectation values for the n ) 6 rotor function are particularly large, which indicates that one D2d orientation is favored (corresponding to one bond aligned with the z axis and the opposite bond aligned with the x axis). More information that can quantify the reorientation dynamics is available from the autocorrelation of the rotor functions
Rn(t) ) 〈Mn(t) Mn(0)〉/〈Mn2〉
(20)
For illustration, Figure 4 shows averages over all seven autocorrelation functions for both temperatures. The long-term memory of the phase is indicated by the nonzero final value of the autocorrelation at long times ≈〈Mn〉2/〈Mn2〉. (The horizontal lines represent the final values for the two temperatures estimated from the values of 〈Mn〉 and 〈Mn2〉.) On the other hand, the rapid decay to the vicinity of the final value represents the short-term memory loss. From the graph, the time scale for this decay is in the 0.4 ps range, which is consistent with the estimated angular velocity correlation time discussed earlier. IV. Dynamics of the Nitromethane Crystal The nitromethane crystal (CH3NO2) was studied as an example of a system in which rotational motion of a methyl group plays an important role. In the gas phase, the barrier against rotation of the unit about the C-N axis is very small (about 6 cal/mol).50 This raises the possibility that in the solid state (melting point 301.7 K) the barrier might remain small
Figure 5. Atomic velocity power spectrum for solid nitromethane at T ) 285 K.
enough that rotational reorientation could be observed. Furthermore, in earlier neutron scattering experiments with an energy resolution of 0.053 meV, no definite resolution of the tunnel splitting of rotational states was possible,51 suggesting that tunneling phenomena have a negligible effect on the rotational motion at temperatures near the melting point. Hence, the approximation of classical nuclei is sufficient to study the dynamics of the nitromethane crystal, and the techniques described in section IIA can be applied. Ab initio MD simulations were performed using two unit cells (eight molecules) of the crystal (orthorhombic, space group P21212152). The lattice parameters were taken to be a ) 5.404 Å, b ) 6.513 Å, and c ) 8.996 Å, which allows for a small thermal expansion from the T ) 0 values. Core electrons were treated using Vanderbilt’s ultrasoft pseudopotential scheme,22 and a plane wave cutoff of 35 Ry was used. Exchange was treated using the gradient corrected functional of Becke,45 while correlation was treated within the LDA. The electron mass parameter µ ) 600 au, and the time step δt ) 5 au. The equations of motion were integrated using the reversible CP method based on Gauss’ equations as described in ref 56. In addition, combining the scaled mass scheme56 with nonorthogonal orbitals gives roughly a factor of 3 savings over standard approaches to CP dynamics with Vanderbilt pseudopotentials.84 A trajectory of total length 4 ps was generated. The temperature was set at T ) 285 K, just below the melting point, to allow sufficient rotational reorientation to occur over the course of the simulation. In Figure 5 the frequency spectrum is shown, which is computed from either the velocity component power spectral densities or the Fourier transforms of the velocity autocorrelation functions. This shows the frequencies present in the system. In order to obtain intensities, it is necessary to compute the dipole moment spectra or autocorrelation functions.53 Since the pure modes measured by IR spectroscopy in the crystal (T ) 78 K)54 and liquid (T ) 354 K)55 are not substantially different, some comparison with the results presented here is possible. The spectrum shown in Figure 5 exhibits the same qualitative features of the IR spectra.54,55 The use of ultrasoft pseudopotentials often leads to a small red shift in the observed frequencies. In Figure 5, there is a consistent 2%-6% red shift in comparison to measured spectra. We can attempt to examine the contributions from the motions of the methyl and nitrate groups separately by computing the velocity power spectra of only the hydrogens and only the oxygens, respectively. The resulting hydrogen vibration spectrum retains all the qualitative features of the full spectrum, which means that, at the temperature considered, the rotation of the methyl group is coupled to the other modes of the
12884 J. Phys. Chem., Vol. 100, No. 31, 1996 molecule. Similarly, the oxygen vibration spectrum exhibits many of the features of the full spectrum, except for the highfrequency CH stretching modes. From these results it can be concluded that a one-dimensional model, based only on the rotation angle θ, would likely be insufficient to describe properly the rotational motion of the methyl group at this temperature. V. Crystalline Hydrogen Chloride Dihydrate Hydrogen chloride dihydrate is a monoclinic crystal with space group P21/c. The lattice parameters are9 a ) 3.991 Å, b ) 12.055 Å, c ) 6.698 Å, and β ) 100.58°. The structure contains puckered layers of Cl- ions connected by hydrogen bonds from water molecules. These water molecules are held together in pairs by a short hydrogen bond (2.41 Å), forming H5O2+ ions. Such a crystal could present an interesting opportunity to study quantum contributions to the motion of the proton in the hydrogen bond joining the two water molecules in each H5O2+Cl- unit. Since the H5O2+ ions appear in an ordered, crystalline environment, it is expected that quantum effects may be important. Thus, the methods presented in section IIB were applied with the hope of comparing the quantum and classical CP calculations. The system simulated consisted of a single unit cell of the crystal (containing four H5O2+Cl- units) at a temperature of 200 K. As in the previous example, Vanderbilt pseudopotentials were used to treat core electrons. A plane wave cutoff of 30 Ry and the Becke corrections to the exchange with Lee-YangParr functional (BLYP)89 were used along with an electron mass parameter of µ ) 600 au and a time step of 5 au. Path integral CP calculations were carried out using the staging MD approach with P ) 8. While this value of P is probably too small for the chosen temperature, it is sufficiently large to make a qualitative comparison between the quantum and classical simulation results. The fictitious masses of the ions were chosen to be twice their physical values (i.e., M′I ) 2MI, M′I(s) ) 2MI(s)). Nose´-Hoover chain thermostat lengths of Melec ) 4 and Mion ) 4 and mass parameters determined by 4ωP and ωp were coupled to the electrons and ions, respectively. The thermostated equations of motion, eqs 13, were integrated using the techniques described in ref 14 with a seventh-order Suzuki/ Yoshida integrator for the thermostats. The equations were evolved for 24 000 and 8000 steps for the classical and quantum systems, respectively. We checked the convergence of the simulation by plotting the cumulative average of both the primitive (cf. eq 16) and virial57 energy estimators. Both estimators show convergence within the first 1000 steps for both the quantum and classical simulations. Although we wished to study the properties of the proton in the O-O hydrogen bond, this goal has not yet been achieved. The simulation proceeded smoothly at first, but after approximately 1 ps in the classical simulation the configuration of the system changed. In the true crystalline structure, each H5O2+ molecule hydrogen bonds to four chlorine atoms. In the new shifted configuration, two of the H5O2+ molecules have rotated slightly in such a way that one of each of their hydrogen atoms (which used to be bonded to a chlorine atom) has shifted and now bonds to the other oxygen atom. Thus, an infinite chain of oxygen atoms hydrogen bonded together is formed. Since one of these oxygens is bonded to four hydrogens and the other to three hydrogens, the structure is hydronium that is bonded twice to water (to a water molecule in the same unit cell and the equivalent one in the next cell along the a axis). Although the other two H5O2+ molecules maintain their form, their bonding pattern to the chlorine atoms is different. The
Tuckerman et al. quantum simulation, which was started from this new configuration, showed no signs of returning to the original structure. We minimized both of these configurations and found that while they were both stable the shifted one had an energy that was approximately 0.015 hartree lower. This energy difference is fairly small and probably within the accuracy of our use of DFT; furthermore, it is possible that while the correct structure has a lower energy in nature, the shifted structure may have an energy only slightly higher. The imposition of a fixed simulation cell based on the experimental parameters may also force the model system to assume the new structure. It will be necessary to let the cell parameters vary to find the minimumenergy configuration without this constraint. It may also be necessary to construct a better pseudopotential for chlorine that takes into account more than just the outermost valence electron of the chlorine atom, based on experience gained in the simulation of other systems involving Cl. Although we have found the BLYP functional to be the best available for this work, it may also be inadequate for a system as complex as this one. These results help to stress that great care must be taken when modeling hydrogen-bonded systems. VI. Hydronium and Hydroxyl Ions in Liquid Water Spectroscopic experiments in solution60-62 and diffraction studies of hydrated crystals63 were able to show that charge defects introduced by excess or missing protons in water occur naturally in the form of the hydronium (H3O+) and hydroxyl (OH-) ions, referred to as the “water ions”. Recent ab initio simulations have been able to verify these results.18,64 These simulations were performed using 32 water molecules in a box of length 9.8 Å with an added or missing proton. Vanderbilt ultrasoft pseudopotentials were used to treat the core electrons. Becke’s gradient-corrected functional was used for the exchange, and the correlation was treated within the LDA. The electronic orbitals were expanded in plane waves up to a cutoff of 25 Ry. When solvated, the H3O+ ion is observed to form hydrogen bonds to three neighboring water molecules to make an H9O4+ solvation complex. The existence of this complex had been predicted by Eigen and De Maeyer in the late 1950s.65 These simulations have also predicted that the OH- ion forms a fourfold coordinated solvation complex in which the ligands lie in a single plane that is roughly perpendicular to the OHbond axis. These two “species” will be referred to as primary solvation complexes. Proton migration is observed to occur via the Grotthuss mechanism,66 which is characterized by a series of proton “hopping” events from one oxygen site to another. These hopping events are actually successive transfers of protons through the hydrogen bonds along a contiguous path in the hydrogen bond network. During a Grotthuss event, intermediate complexes are observed to occur. In the H+ case, this complex is the H5O2+ complex in which the excess proton is located midway along the hydrogen bond joining two water molecules. This complex has received considerable attention in the literature (for a review, see ref 67). In the OH- case, the intermediate complex is a tetrahedral H7O4- complex. Figures 6 and 7 show the radial distribution functions gO*O and gO*H, where O* denotes the oxygen atom in H3O+ or OH-. Figure 6 shows the radial distribution functions for each solvation complex (H9O4+ and H5O2+) in comparison to the distribution functions for pure water.44 Similarly, Figure 7 shows the radial distribution functions for the H9O5- and H7O4complexes. In both cases the solvation shells are shorter than the ordinary water-water H bond. However, in the OH- case, these H bonds are somewhat longer, and hence weaker, than those in the H3O+ solvation shell. Rough measurements of the
Molecular Dynamics Simulations
Figure 6. Oxygen-oxygen (top) and oxygen-hydrogen (bottom) radial distribution functions for H+ in water. The distribution functions are measured with respect to the hydronium oxygen (O*). In each figure, the solid line and short dashed line respectively correspond to the subset of configurations in which the H5O2+ and H9O4+ complexes are present. The long dashed curves are the O-O and O-H radial distribution functions for pure water, taken from the simulations of Laasonen et al.84
lifetimes of these complexes in solution as well as studies of them in vacuum suggest that they are stable and hence can be referred to as secondary complexes. A Grotthuss hop is accompanied by a reorganization of the local hydrogen bond network to accommodate the formation of a solvation complex on a new site. Thus, a solvation complex is observed to migrate through the hydrogen bond network, giving rise to a structural diffusion model of proton migration.68-71 The H+ and OH- systems are qualitatively different. In the H+ case, the H9O4+ complex is “active” in the sense that it is the diffusing structure, moving from site to site via formation of the H5O2+ intermediate. In the OH- case, the H9O5complex is “inactive” owing to the presence of a distinct solvation shell around it. The active H7O4- complex forms upon breaking of one of the hydrogen bonds in the plane defined by the OH- bond axis. It is this active complex that is observed to diffuse through the network. Further analysis of the simulated trajectories have allowed the rate-limiting step in proton transfer for each system to be identified. In the case of H+, a change in the coordination number of a water molecule in the first solvation shell of the hydronium ion (i.e., second solvation shell dynamics) from 4 to 3 allows the hydronium structure to migrate to the lowcoordinated site. In the case of OH-, the analysis reveals that purely geometrical factors determine which of the two solvation complexes permits structural diffusion. In the H9O5- complex, the angle between the ligands and the OH- bond axis is roughly 90°. Thus, proton transfer from water molecules in the first solvation shell to the OH- site is hindered since a new water molecule with an angle of 90° would form as a result of the transfer. By contrast, the H7O4- complex has an approximately
J. Phys. Chem., Vol. 100, No. 31, 1996 12885
Figure 7. Oxygen-oxygen (top) and oxygen-hydrogen (bottom) radial distribution functions for OH- in water. The distribution functions are measured with respect to the hydroxyl oxygen (O*). In each figure, the solid line and short dashed line respectively correspond to the subset of configurations in which the H7O4- and H9O5complexes are present. The long dashed curves are the O-O and O-H radial distribution functions for pure water, taken from the simulations of Laasonen et al.84
105° angle between the OH- bond axis and the ligands and can easily accept the transferring proton, thus forming a water molecule. These geometrical differences are responsible for the distinction in activity of the two complexes. There is, in addition, a secondary rate-limiting step in the OH- case. In the active state, the H7O4- complex is hydrogen-bonded to three neighboring water molecules, each of which is solvated by four neighbors. As in the H+ case, a fluctuation in the coordination number of one of these molecules from 4 to 3 must precede proton transfer so that the active complex can be re-formed at the new site, and once again second solvation shell dynamics plays a role in the proton transfer process. On the other hand, since the dominant rate-limiting step is the formation of the active complex from the inactive one (which also depends on the breaking of a hydrogen bond), second solvation shell dynamics can only be said to have a secondary effect here. VII. Charged Water Clusters As a prelude to understanding the possible role of quantum effects in ion solvation charged water clusters, two prototypical species, H5O2+ and H3O2-, have been investigated.19 The simulation studies employed the pseudopotential scheme of Martins and Troullier72 with a cutoff of 70 Ry and the so-called BLYP functional to treat exchange and correlation.45,73 Each complex was placed at the center of a cubic box with length 8 Å to which cluster boundary conditions were applied.85 Nuclear quantum effects have been incorporated using the path integral CP techniques described above.19 The same clusters have been studied by traditional quantum chemical techniques.74-76 Figures 8 and 9 show contours of the two-dimensional distribution function of ROO (the oxygen-oxygen distance) and
12886 J. Phys. Chem., Vol. 100, No. 31, 1996
Tuckerman et al.
Figure 8. Contours of the two-dimensional distribution function of the oxygen-oxygen distance, ROO, and the asymmetric stretch coordinate, d12 ) rO1H* - rO2H*, for the isolated H5O2+ complex at 300 K. H* denotes the shared proton. The top curve shows the results from a CP simulation using the classical nuclei approximation. The bottom curve shows the results from a path integral CP simulation.
Figure 9. Contours of the two-dimensional distribution function of the oxygen-oxygen distance, ROO, and the asymmetric stretch coordinate, d12 ) rO1H* - rO2H*, for the isolated H3O2- complex at 300 K. H* denotes the shared proton. The top curve shows the results from a CP simulation using the classical nuclei approximation. The bottom curve shows the results from a path integral CP simulation.
the asymmetric stretch coordinate (d12 ) rO1H* - rO2H*, where H* is the shared proton). These figures show the results for both the classical and path integral simulations of each complex. In the path integral simulations, P ) 16 imaginary time discretizations were used. The results of the simulations suggest that, at a temperature of 300 K, quantum effects in the H5O2+ cluster do not significantly alter the qualitative classical picture, in which the proton is most likely to be found midway between the two oxygen atoms. Comparing quantum and classical simulations at the same temperature, one can see that the quantum distribution function describing the spatial location of the proton is somewhat wider than the classical one but that the general features of these distributions are the same. In the H3O2cluster, however, quantum effects appear to alter even the qualitative picture. The classical distribution has a definite bimodal character in which the proton is located on either one oxygen site or the other. However, in the quantum picture, the proton is most likely to be located midway between the oxygen atoms, as in the case of positive charge. These results are consistent with the picture obtained from the ion solvation simulations. The hydrogen bonds in the first solvation shell of the hydronium ion are short, and therefore the potential surface felt by the proton is not likely to have a large barrier. In the hydroxyl simulations, however, the hydrogen bonds in the solvation shell of both the active and inactive complexes are long (∼2.5-2.7 Å), probably leading to the presence of a significant barrier. The results of the cluster calculations suggest that quantum effects soften this barrier
enough that the proton has a high probability of being located midway between the oxygens. VIII. Summary and Future Outlook This article has focused on a few applications of ab initio MD simulations to condensed phase problems linked to physical chemistry. Codes for classical equilibrium MD simulations of such systems have been available commercially for some time. Their use is spreading and, while they are gaining in sophistication, the field of nonequilibrium molecular dynamics (NEMD)83 still presents unique challenges.78-82 The methodology of purely classical equilibrium MD also continues to evolve.77 The topics we have discussed show that it is now possible to carry out MD simulations of condensed phase systems without recourse to empirical parameters to describe the interatomic interactions. Even the quantum mechanical behavior of the nuclear motion is now accessible via path integral CPMD methods. Although the system sizes that can be ordinarily handled at the present time (typically, a few hundred atoms) are modest compared to classical MD simulations with empirical potentials, ab initio MD methods will no doubt grow in importance. For example, in the future to investigate the function of an enzyme, it will likely become a routine procedure to use empirical potentials to handle the protein and its host environment, while using ab initio techniques for the active site of interest. A future challenge will be to develop more rigorous approaches and simulation methodologies for such mixed quantum-classical schemes.
Molecular Dynamics Simulations Acknowledgment. It is a pleasure to thank our colleagues around the world for sharing their new developments and most recent results prior to publication. This research described herein received generous support from the National Science Foundation through Grants CHE 92-24536, STI 94-13373, and DMR 91-20668. M.E.T. thanks the NSF for a Postdoctoral Research Associateship in Computational Science and Engineering (ASC95-04076) and P.J.U. the Natural Sciences and Engineering Research Council of Canada for a Postdoctoral Fellowship. References and Notes (1) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (2) Galli, G.; Parrinello, M. In Computer Simulations in Materials Science; Meyer, M., Pontikis, V., Eds.; Kluwer: Dordrecht, 1991. (3) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (4) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B86. (5) Simon, A.; Bormann, H.; Craubner, H. Phosphorus Sulfur Relat. Elem. 1987, 30, 507. (6) Kamitakahara, W. A.; Gompf, F.; Neumann, D. A. Mater. Res. Soc. Symp. Proc., in press. (7) Ungar, P. J.; Laasonen, K.; Klein, M. L. Can. J. Phys. 1995, 73, 710. (8) Trommsdorff, H. P. In Disorder in Materials; Powell, B. M., Ed.; 1995; Vol. 10, No. 2. (9) Lundgren, J.-O.; Olovsson, I. Acta Crystallogr. 1967, 23, 966. (10) Parrinello, M.; Rahman, A. J. Chem. Phys. 1984, 80, 860. (11) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. J. Chem. Phys. 1992, 97, 2635. (12) Marx, D.; Parrinello, M. Z. Phys. B 1994, 95, 143. (13) Marx, D.; Parrinello, M. J. Chem. Phys., submitted. (14) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. J. Chem. Phys., submitted. (15) Chandler, D.; Wolynes, P. G. J. Chem. Phys. 1981, 74, 4078. (16) Hall, R. W.; Berne, B. J. J. Chem. Phys. 1984, 81, 3641. (17) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (18) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150. (19) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. Manuscript in preparation. (20) Carloni, P.; Blo¨chl, P.; Parrinello, M. J. Phys. Chem. 1995, 99, 1338. (21) Sagnella, D.; Laasonen, K.; Klein, M. L. Biophys. J., in press. (22) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892. (23) Hutter, J.; Tuckerman, M.; Parrinello, M. J. Chem. Phys. 1995, 102, 859. (24) Bachelet, G. B.; Hamann, D. R.; Schlu¨ter, M. Phys. ReV. B 1982, 26, 4199. (25) See, for example: Fetter, A. L., Walecka, J. D. Theoretical Mechanics of Particles and Continua; McGraw-Hill: New York, 1980; Chapter 3. (26) Laasonen, K.; Pasquarello, A.; Car, R.; Lee, C.; Vanderbilt, D. Phys. ReV. B 1993, 47, 10142. (27) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (28) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (29) Stich, I.; Car, R.; Parrinello, M.; Baroni, S. Phys. ReV. B 1989, 39, 4997. (30) Galli, G.; Parrinello, M. Phys. ReV. Lett. 1992, 24, 3547. (31) Mauri, F.; Galli, G. Phys. ReV. B 1994, 50, 4316. (32) Cao, J.; Berne, B. J. J. Chem. Phys. 1993, 99, 2902. (33) Pollock, E. L.; Ceperley, D. M. Phys. ReV. B 1984, 30, 2555. (34) Sprik, M.; Klein, M.; Chandler, D. J. Chem. Phys. 1985, 83, 3942. (35) Sprik, M.; Klein, M.; Chandler, D. Phys. ReV. B 1985, 31, 4234. (36) Coker, D. F.; Thirumali, D.; Berne, B. J. J. Chem. Phys. 1987, 86, 5689. (37) Martyna, G. J.; Cao, J. J. Chem. Phys., submitted. (38) Corbridge, D. E. C.; Lowe, E. J. Nature 1952, 170, 629. (39) Corbridge, D. E. C. Phosphorus, an Outline of Its Chemistry, Biochemistry, and Technology, 4th ed.; Elsevier: Amsterdam, 1990; pp 61-66. (40) Ceperley, D. M.; Alder, B. J. Phys. ReV. Lett. 1980, 45, 566. (41) Perdew, J. P. Phys. ReV. B 1981, 23, 5048. (42) Ortiz, G.; Ballone, P. Phys. ReV. B 1991, 43, 6376. (43) Sim, F.; St. Amant, A.; Papai, I.; Salahub, D. R. J. Am. Chem. Soc. 1992, 114, 4391. (44) Laasonen, K.; Parrinello, M.; Car, R.; Lee, C.; Vanderbilt, D. Chem. Phys. Lett. 1993, 207, 208.
J. Phys. Chem., Vol. 100, No. 31, 1996 12887 (45) Becke, A. D. Phys. ReV. A 1988, 38, 3098; J. Chem. Phys. 1992, 96, 2155. (46) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (47) Hamann, D. R. Phys. ReV. B 1989, 40, 2980. (48) Bosworth, Y. M.; Clark, R. J. H.; Rippon, D. M. J. Mol. Spectrosc. 1973, 46, 240. (49) Klein, M. L.; McDonald, I. R.; Ozaki, Y. J. Chem. Phys. 1983, 79, 5579; Phys. ReV. Lett. 1982, 48, 1197. (50) Tannenbaum, E.; Myers, R. J.; Gwinn, W. D. J. Chem. Phys. 1956, 25, 42. (51) Trevino, S. F.; Rymes, W. H. J. Chem. Phys. 1980, 73, 3001. (52) Trevino, S. F.; Prince, E.; Hubbard, C. R. J. Chem. Phys. 1980, 73, 2996. (53) Wilson reference. (54) McKean, D. C.; Watt, R. A. J. Mol. Spectrosc. 1976, 61, 184. (55) Hill, J. R.; Moore, D. S.; Schmidt, S. C.; Storm, C. B. J. Phys. Chem. 1991, 95, 3037. (56) Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1994, 101, 1302, 1316. (57) Herman, M. F.; Bruskin, E. J.; Berne, B. J. J. Chem. Phys. 1982, 76, 5150. (58) Benoit, M.; Bernasconi, M.; Focher, P.; Parrinello, M. Phys. ReV. Lett. 1996, 76, 2934. (59) Komatsuzaki, T.; Ohmine, I. Chem. Phys. 1994, 180, 239. (60) Gigue`re, P. A.; Turrell, S. Can. J. Chem. 1976, 54, 3477. (61) Zundel, G. In The Hydrogen Bond; Schuster, P., Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vol. 2, Chapter 15. (62) Bratos, S.; Ratajczak, H.; Viot, P. In Hydrogen Bonded Liquids; Dore, J. C., Texeira, J., Eds.; Kluwer: Dordrecht, 1991. (63) Lundgren, J.-O.; Olovson, I. in The Hydrogen Bond; Schuster, P., Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vol. 2, Chapter 10. (64) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Phys. Chem. 1995, 99, 5749. (65) Eigen, M.; De Maeyer, L. Proc. R. Soc. London, Ser. A 1958, 247, 505. (66) Daneel, H. Z. Elektrochem. 1905, 11, 249. (67) Scheiner, S. Annu. ReV. Phys. Chem. 1994, 45, 23. (68) Conway, B. E. In Modern Aspects of Electrochemistry; Bockris, J. O’M., Conway, B. E., Eds.; Butterworths: London, 1964; Vol. 3, p 43. (69) Erdey-Gruz, T.; Lengyel, S. In Modern Aspects of Electrochemistry; Bockris, J. O’M., Conway, B. E., Eds.; Butterworths: London, 1964; Vol. 12, p 3. (70) Stillinger, F. H. In Theoretical Chemistry, AdVances and PerspectiVes; Eyring, H., Henderson, D., Eds.; Academic Press: New York, 1978; Vol. 3, p 177. (71) Kreuer, K. D.; Dippel, Th.; Hainovsky, N. G.; Maier, J. Ber. BunsenGes. Phys. Chem. 1992, 96, 1736. (72) Troullier, N.; Martins, J. Phys. ReV. B 1991, 43, 1993. (73) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University: Oxford, 1989. (74) Grimm, A. R.; Backsay, G. B.; Haymet, A. D. J. Mol. Phys. 1995, 86, 369. (75) Lee, Y. T. 1995. (76) Wei, D.; Salahub, D. R. J. Chem. Phys. 1994, 101, 7633. (77) Martyna, G. J.; Tuckerman, M. L.; Tobias, D. J.; Klein, M. L. Mol. Phys., submitted. (78) Edberg, R.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1986, 84, 6933. (79) Edberg, R.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1987, 86, 4555. (80) Daivis, P. J.; Evans, D. J. J. Chem. Phys. 1994, 100, 541. (81) Mundy, C. J.; Siepmann, J. I.; Klein, M. L. J. Chem. Phys. 1995, 102, 3376. (82) Mundy, C. J.; Siepmann, J. I.; Klein, M. L. J. Chem. Phys. 1995, 103, 10192. (83) Evans, D. J.; Morriss, G. P. Statistical Mechanics of Nonequilibrium Liquids; Academic Press: New York, 1990, and references therein. (84) Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1993, 99, 9080. (85) Marx, D.; Fois, E.; Parrinello, M. Int. J. Quantum Chem. 1996, 57, 655. (86) Silvestrelli, P. L.; Alavi, A.; Parrinello, M.; Frenkel, D. Europhys. Lett. 1996, 33, 551. (87) Bernasconi, M.; Parrinello, M.; Chiarotti, G. L.; Focher, P.; Tosatti, E. Phys. ReV. Lett. 1996, 76, 2081. (88) Bernasconi, M.; Chiarotti, G. L.; Focher, P.; Parrinello, M.; Tosatti, E. Phys. ReV. Lett., submitted. (89) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785.
JP960480+
