Ewald artifacts in liquid state molecular dynamics simulations Paul E. Smith and B. Montgomery Pettitt Department of Chemistry, University of Houston, Houston, Texas 77204-5641 (Received 23 May 1996; accepted 4 June 1996) An investigation into the effects of the anisotropic nature of the Ewald potential for the treatment of long range electrostatic interactions in liquid solutions has been performed. The rotational potential energy surface for two simple charge distributions, and a small protein, have been studied under conditions typically implemented in current biomolecular simulations. A transition between hindered and free rotation is observed which can be modeled quantitatively for simple charge distributions. For most systems in aqueous solution, the transition involves an energy change well below kBT. It is argued that, for solvents with a reasonably high relative permittivity, Ewald artifacts will be small and in many cases may be safely ignored. © 7996 American Institute of Physics. [S0021-9606(96)02034-X] I. INTRODUCTION The correct treatment of long range electrostatic interactions is of considerable importance for the reliability and reproducibility (precision) of modern computer simulations of polar/charged systems. Approximate methods,1 while often simple and computationally cheap, have been shown to produce unacceptable artifacts which bring into question the validity of many of the results so obtained.2-4 The Ewald technique is a well established method for a rigorous treatment of electrostatic interactions in periodic systems which contain a set of full or partial atomic charges.5 While the method has been successfully applied to a variety of systems in condensed phase physics, its use for simulating liquids6'7 and, in particular, biomolecular systems, has only recently emerged.8'9 Some of the reasons for the slow acceptance of the technique include the added expense and complexity involved in implementing the method, and the possibility of artificial enhanced periodicity effects for finite size systems.10-12 Progress has been made in reducing the computational effort required in implementing the Ewald technique.13-15 However, the question of artificial periodicity effects has received less attention.16 Artificial periodicity effects can be thought of by simply comparing the force lines generated by a simple charge distribution, a dipole for example, in isolation or at infinite dilution (normally the target conditions), with the field lines obtained for a infinite array of dipoles. The field lines differ at large distances and this has been a perceived problem with the technique. Here, we consider the significance of this artifact for a number of test cases. Our study is focused on two simple model charge distributions, a linear dipole and quadrupole, together with a study of the protein, bovine pancreatic trypsin inhibitor (BPTI). For all three systems we were interested in the effect of the anisotropic nature of the Ewald potential on the rotational properties of the system. A linear dipole, for instance, will tend to align preferentially along one of the axes of the periodic cell; a ferroelectric configuration of the system. This study attempts to determine under what conditions the alignment will dominate to a degree that free thermal rotation is no longer possible. The evidence suggests that for most typical liquid state biomolecular simulations Ewald artifacts will be small and can be safely ignored. In Section II, we develop equations which prove useful in modeling the anisotropic nature of the Ewald potential for simple charge distributions. We show the dependence of the rotational potential energy surface on charge, charge separation, relative permittivity, and box length. The results for a linear dipole and quadrupole, together with those for BPTI, are then presented. Finally, the implications for simulations using the Ewald potential in liquid state simulations of mac-romolecules are discussed in the conclusions. II. THEORY Here, we develop equations to model the orientational effects of simple linear dipoles and quadrupoles observed when using Ewald electrostatics. The Ewald electrostatic potential energy of a periodic system of TV charges ({q}) with a total charge of Q in a cubic box of length L and relative permittivity er, surrounded by a dielectric continuum of relative permittivity =° is given by5'17 1 ^ q,q V= 4Trere, N 2 ferfc(a|r,,|) • 3 (x2y2z2) As the V0 and V2 terms do not depend on the relative orientation of the dipole in the box, the leading term for the rotational energy is given by, Va=V. 1 277-r ?4! Va(6,), (8) where Va( 6, ) describes the variation in energy (spherical polar coordinates) on rotation of the dipole. We will show that further terms provide a negligible contribution to the relative rotational potential energy surface. We now switch our attention to the properties of this potential energy surface. The transition from hindered to free rotation may be thought of as analogous to that of a rotational melting (phase) transition. However, we know that heat capacity of such systems changes dramatically at the melting temperature. Hence, we have chosen to follow the heat capacity curves as a function of temperature, generated from the rotational potential energy surface, in an attempt to characterize the transition from hindered to free rotation. For a system in the canonical (NVT) ensemble the constant volume heat capacity is given by, CV(T)- 1 kuT- i(E-(E))2), and the derivative with respect to temperature by T4^f=^T((E-{E))3)-2T3Cv, (9) (10) where kB is the Boltzmann constant, E is the total energy of the system, and the angular brackets denote an ensemble average. We denote a maximum in Cv as defining a transition temperature Tr. Below Tr rotation is severely restricted due to orientation along the field lines generated by the lattice, while above Tr rotation is relatively free and lattice artifacts will be small. Tr may be obtained from the expres- + 15D2(x4y2 + x4z2 + x2y4 + x2z4 + y4z2 + y2z4)l where x = rxl\r\, y = ryl\r\, and z = rzl\r\ are unit vectors defining the orientation of the dipole and, E„n2a, Cv E„na, C2=S Enn2an\, £>i=S Enn6a, D2= Zj Enn4anl, D3= X E„n2ji2Rn2 (7) for which a ¥= /3 ¥= y. The parameters A, B, C, and D are solely dependent on ah and, by symmetry, are invariant on interchange of a, /3, or y. In addition, the terms within the square brackets in Eq. (6) are only dependent on the relative orientation of the dipole in the central box. The above expansion is similar to previous Ewald approximations.18-20 2kBTr ((E-(E)T/)Tr ~~{{e-{e)tý)t; (ii) where the subscript Tr denotes that ensemble averages are determined at Tr. One can rewrite the expressions for the heat capacity and Tr using our expression for the energy as a function of the rotation of the dipole. It will be useful to define a reduced energy e such that, 1 q2 1 /27r|r|\4 ATTere0 ttL 4! (12) Considering just the configurational contribution to the rotational heat capacity in terms of the reduced temperature T* = knTlB we have, C0(T*) 1 ^2((Va-(Va))2) (13) and J. Chem. Phys., Vol. 105, No. 10, 8 September 1996 P. E. Smith and B. M. Pettitt: Ewald artifacts in molecular dynamics simulations 4291 2T* ~~(V2a)T* (14) (a) if we chose (Va)T* as our zero of energy. For our interest the ensemble average of X is given by (X)- J2^^X(B,cf>)e - V„IT* sin B dB d

. The calculation is for charges +1/— 1 separated by 1.0 nm in a box of length 2.5 nm and a relative permittivity of 1.0. (a) Calculated using the full Ewald potential, (b) just the V4 term, and (c) the V4 + V6 terms. Contours are displayed at 0.05 (dot-dashed), 0.10 (dotted), 0.15 (dashed), and 0.20 (long-dashed). 1(c) display the maps obtained for the dipole using the V4 term and the V4 + V6 terms, respectively. The maps are very similar in the region of the minima with small differences in J. Chem. Phys., Vol. 105, No. 10, 8 September 1996 4292 P. E. Smith and B. M. Pettitt: Ewald artifacts in molecular dynamics simulations FIG. 2. Constant volume heat capacity curves for a linear dipole (solid) and quadrapole (dashed) in an Ewald field as a function of the reduced temperature. curvature when approaching the maximum. The relative energy of the maximum in Figs. 1(b) and 1(c) is 0.251 and 0.247, respectively. The results obtained with just the V4 term are sufficiently accurate for our purposes. The rotational energy surface for a linear quadrupole may be obtained by inverting the dipole potential surface and scaling by 3/4. Hence, linear quadrupoles tend to orientate with their axis pointing towards the corner of the cell. In Fig. 2 we show the variation of Cv with reduced temperature for the linear dipole and quadrupole systems. Both display a maximum and a slow decay to zero as T* approaches infinity. Values of T* for the dipole and quadrupole were 0.037 and 0.007, respectively, or 131 and 24 K for the conditions described in the methods section. Above T* there are no large potential energy barriers to rotation. However, this does not imply free rotation. As a more realistic measure of free rotation we have taken the value of T* for which the average rotational kinetic energy (kBT) is equal to twice the energy difference between the maximum and minimum. Using this approach, free rotation occurs at reduced temperatures of 0.502 (1775 K) and 0.377 (1332 K) for the dipole and quadrupole, respectively. The numbers in parentheses correspond to a relative permittivity of unity, which is an extreme case. Both Tr and the temperature for free rotation are inversely scaled by the relative permittivity. Therefore, under aqueous conditions both temperatures will be reduced by almost two orders of magnitude. Further justification for dominance of the V4 term for linear dipoles and quadrupoles comes from the observed dependence of Tr on \ r\ and L. Log-log plots for variation of these two parameters give straight lines with slopes of 4.05 ±0.01 and —5.09±0.01, respectively, and correlation coefficients of 0.999 99 for both. It should be noted here, that this does not imply that higher terms (V6 and beyond) are small. They are, in fact, large in magnitude but display only a small dependence on rotation. The population of BPTI rotational states at 300 K (averaged over 4>) is displayed in Fig. 3 for two different relative permittivities. Above a relative permittivity of 10.0, the dif- 180 - 90 90 180 270 360 180 90 90 180 270 360 FIG. 3. Population contours (%) for BPTI at 300 K for relative permittivities of (a) 1.0 and (b) 10.0 as a function of the Euler angles 0 and if/. Results have been averaged over . In (a) contours are at 0.05 (dotted) and 0.10 (dashed) and in (b) contours are at 0.005 (dotted) and 0.01 (dashed). ference in population between the different rotational states is very small. This is also illustrated in the Cv curve for BPTI shown in Fig. 4. While the decay to zero is very slow for low values of the permittivity, this is increased dramatically with larger permittivities. For a relative permittivity of 1.0 Tr is found to be 135 K. The difference in energy between the maximum and minimum on the BPTI rotational potential energy surface is 50.4 kj/mol which translates to a temperature for free rotation of 8084 K. However, for a relative permittivity of 80.0 the values of Tr, the energy difference, and the free rotation temperature are reduced to 2 K, 0.6 kj/mol and 101 K, respectively. V. CONCLUSIONS We have shown that rotational Ewald artifacts due to artificial periodicity in liquid simulations can be observed and modeled for simple charge distributions. Above a certain temperature these artifacts become negligible. For a simple dipole and quadrupole this temperature varies as q2r4L~5/er. The first three variables are all known before a simulation is attempted. However, the value of er, which is of major importance, can only be approximated. For a molecular dipole/quadrupole in solution the solvent molecules will, on average, align to oppose the solute dipole. This imparts dielectric screening between the molecule in the central cell and images in neighboring cells. This screening will be related to the dielectric constant of the solvent. For water as solvent a relative permittivity of 80.0 is more appropriate and leads to the conclusion that, for both the model systems and BPTI, Ewald artifacts are completely negligible for most properties of interest. While we expect more sophisticated J. Chem. Phys., Vol. 105, No. 10, 8 September 1996 P. E. Smith and B. M. Pettitt: Ewald artifacts in molecular dynamics simulations 4293 200 300 T(K) FIG. 4. Constant volume heat capacity curves for BPTI as a function of temperature and relative permittivity. Curves are for relative permittivities of 1.0 (solid), 10.0 (dotted), and 80.0 (dashed). approaches for modeling the dielectric screening effect to change this result slightly, the general conclusion should still hold. ACKNOWLEDGMENTS We would also like to thank the Robert A. Welch Foundation The National Institutes of Health and the National Science Foundation for partial support. 'P. E. Smith and B. M. Pettitt, J. Phys. Chem. 98, 9700 (1994). 2C. L. Brooks III, B. M. Pettitt, and M. Karplus, J. Chem. Phys. 83, 5897 (1985). 3P. E. Smith and B. M. Pettitt, J. Chem. Phys. 95, 8430 (1991). 4H. Schreiber and O. Steinhauser, Biochemistry 31, 5856 (1992). 5M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford, 1987). 6J. Andersen, J. J. Ullo, and S. Yip, J. Chem. Phys. 87, 1726 (1987). 7H. J. Strauch and P. T. Cummings, Mol. Simul. 2, 89 (1989). 8T. R. Forester and I. R. McDonald, Mol. Phys. 72, 643 (1991). 9 P. E. Smith, L. X. Dang, and B. M. Pettitt, J. Am. Chem. Soc. 113, 67 (1991). 10D. N. Card and J. P. Valleau, J. Chem. Phys. 52, 6232 (1970). 11 J. P. Valleau and S. G. Whittington, in Statistical Mechanics A. Modern Theoretical Chemistry, edited by B. J. Berne (Plenum, New York, 1977), pp. 137-168. 12W. F. van Gunsteren, H. J. C. Berendsen, and J. A. C. Rullmann, Faraday Discuss Chem. Soc. 66, 58 (1978). 13T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 14D. Fincham, Mol. Simul. 13, 1 (1994). 15P. E. Smith and B. M. Pettitt, Comput. Phys. Commun. 91, 339 (1995). 16F. Figueirido, G. S. Del Buono, and R. M. Levy, J. Chem. Phys. 103, 6133 (1995). 17M. P. Tosi, in Solid State Physics, edited by F. Seitz and D. Turnbull (Academic, New York, 1964), Vol. 16, pp. 107-113. 18J. P. Hansen, Phys. Rev. A. 8, 3096 (1973). 19W. L. Slattery, G. D. Doolen, and H. E. DeWitt, Phys. Rev. A. 21, 2087 (1980). 20D. J. Adams and G. S. Dubey, J. Comput. Phys. 72, 156 (1987). 21K. D. Bemdt, P. Güntert, L. Orbons, and K. Wüthrich, J. Mol. Biol. 227, 757 (1992). 22W. F. van Gunsteren and H. J. C. Berendsen, Groningen Molecular Simulation (GROMOS) Library Manual (Biomos, Groningen, 1987). J. Chem. Phys., Vol. 105, No. 10, 8 September 1996 The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp Copyright of Journal of Chemical Physics is the property of American Institute of Physics and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.