Available online at www.sciencedirect.com
. CHEMICAL 8CieHCB@DIRecT9 PHYSICS
LETTERS
ELSEVIER Chemical Physics Letters 406 (2005) 49-53
www. else vier. com/locate/cplett
A severe artifact in simulation of liquid water using a long cut-off length: Appearance of a strange layer structure
Yoshiteru Yonetani
Center for Promotion of Computational Science and Engineering, Japan Atomic Energy Research Institute, 8-1 Umemidai, Kizu-cho, Soraku-gun, Kyoto 619-0215, Japan
Received 27 December 2004; in final form 16 February 2005 Available online 17 March 2005
Abstract
We report that a severe artifact appeared in molecular dynamics simulation of bulk water using the long cut-off length 18 A. Our result shows that increasing the cut-off length does not always improve the simulation result. Moreover, the use of the long cut-off length can lead to a spurious result. It is suggested that the simulation of solvated biomolecules using such a long cut-off length, which has been often performed, may contain an unexpected artifact. © 2005 Elsevier B.V. All rights reserved.
Many artifacts caused by truncating long-ranged intermolecular interactions in molecular simulations, i.e., electrostatic interactions, at a finite distance have been observed for various systems such as pure liquids [1-7], ionic solutions [8,9], and solvated biomolecules (e.g., peptides [10], proteins [11-13], DNA [14] and membranes [15]). In spite of this fact, such a cut-off scheme is still used [16-21], especially when large systems composed of biomolecules and water are simulated. When adopting the cut-off scheme, one would expect that natural behavior of the molecular system is well reproduced by increasing a cut-off length to a distance where a correlation between molecules disappears [22] (e.g., a distance where any peaks are not seen in a radial distribution function).
Various problems of treating the long-ranged interactions in simulations of biomolecular systems containing water are reviewed by Lounnas [22J. According to his statement, the cut-off length 9-12 A, which has been used conventionally, is too short to reproduce real behavior of the system, and thus, it is supposed that a longer length of 15-20 A will be required for obtaining
E-mail address: yonetani@apr.jaeri.go.jp
a reliable result. In fact, several examples [23-25] in which favorable results were obtained by use of such a long cut-off length have been reported. However, the result we obtained from molecular dynamics (MD) simulations of bulk water using the long cut-off length 18 A showed an opposite tendency; increasing the cut-off length did not improve the simulation result, and to make matters worse, it led to an unfavorable, more artificial result. In this Letter, we show this result, and point out a problem that one can encounter when using such a long cut-off length.
We performed MD simulations using a cubic MD cell under a three-dimensional periodic boundary condition. The water model used was TIP3P [26], which is the same as that of the previous studies [16,23,24] using the long cut-off length of 18 A. The water molecule was treated as a rigid body by the SHAKE constraint [27]. The equations of motion were numerically solved with the time-step 2 fs. Initial structures were prepared by the LEaP module of the Amber7 software [28]. After the energy minimization of the initial structures, MD simulations were started. The temperature was initially set to 10 K, which was changed to 100, 200, 300, 350, and 298 K every 10 ps. Thus, the temperature reached
0009-2614/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.02.073
50
Y. Yonetani I Chemical Physics Letters 406 (2005) 49-53
298 K when 50 ps passed. After that, the simulations were continued for 1 ns under the condition. The pressure was kept at 1 atm through the simulations. The temperature and pressure were controlled by the weak coupling method [29] with the coupling constants of 0.1 and 0.2 ps, respectively. The group-based cut-off scheme was adopted to avoid the splitting of the water dipole [30]. The cut-off length of 18 A was used for both Coulomb and van der Waals (vdW) interactions. For comparison, we performed two other simulations using different treatment of the interactions, i.e., 9 A cut-off and the particle mesh Ewald (PME) method [31] with the real space cut-off of 9 A. For these three simulations, 2201 water molecules were used. The 18 A cut-off simulation was also performed for larger systems of 3316 and 4736 molecules to examine the system size dependence of the simulation results. These MD simulations were performed by the Gibbs module of the Amber7.
Fig. 1 shows the results of the three kinds of simulations with the 2201 molecules, i.e., 18 A cut-off, 9 A cut-off, and the PME. In the figure, 9 is the angle between the dipoles of water molecules, and the angular bracket (...) denotes the average over the all molecule-pairs in the MD cell at a moment. The orienta-tional order of the water molecules is quantified by (|cos0|). For example, this quantity becomes 0.5 when the molecular orientations are at random. In the simulation with 18 A cut-off (Fig. la), (|cos0|) was suddenly increased at about 100 ps, showing that a phase transition accompanied by a drastic change of the molecular orientations occurred. On the other hand, in the other simulations (Fig. lb,c), such an event was not observed and the value of (|cos0|) was almost 0.50 through the simulations.
The change of the vdW and Coulomb potential energies in the above phase transition is shown in Fig. 2. A large drop was observed for the Coulomb potential, but not observed for the vdW one, suggesting that the transition is caused by the reduction in the Coulomb potential.
The structure obtained by the transition is shown in Fig. 3a. In the figure, we can see a structure that is significantly different from the real water structure. The structure is made of two kinds of layers; one is almost occupied by oxygen atoms, and another is almost occupied by hydrogen atoms. The width of the double layer is about 18 A. This can be approximately measured by comparing with the cell edge length of ^40 A. The width of the double layer is almost identical with the cut-off length 18 A, implying that any relation between them exists. For comparison, structures from the other simulations (i.e., 9 A cut-off and the PME) are also shown in Fig. 3b,c. In these cases, no layer structure was observed. Therefore, the appearance of the artificial layer structure in the 18 A cut-off simulation is owing to the long cut-off length.
0.53
0.52
0.51
0.49
0.53
0.52
0.51 -
0.49
0.53
0.52
0.51
o _o
V
400 600 Time [ps]
1000
400 600 Time [ps]
1000
0.49
400 600 Time [ps]
1000
Fig. 1. Time evolution of (|cosň|) in the simulations with 2201 water molecules: (a) 18 ä cut-off; (b) 9 ä cut-off; (c) PME.
As shown above, the layer structure has the transla-tional periodicity of ^18 A, which is close to half of the cell edge length ^40 A. In general, when MD simulations are performed under a periodic boundary condition, some restrictions are imposed on the resulting structures; even if a structure is energetically favorable, it is not allowed to exist without agreement between the periodicity of the resulting structure and that of the MD cell. Thus, whether the layer structure appears or not is also influenced by the cell size. This is confirmed by the following results of 18 A cut-off simulations with 3316 and 4736 molecules. In each case, the cell edge length obtained after the equilibration at
Y. Yonetani I Chemical Physics Letters 406 (2005) 49-53
51
100 120 Time [ps]
o E
"co
o
c
B o
Q.
"D
>
100 120 Time [ps]
Fig. 2. Time evolution of the Coulomb and vdW potential energy in the 18 A cut-off simulation with 2201 water molecules. The arrow indicates the change of the Coulomb potential in the phase transition.
298 K and 1 atm was 46 and 52 Ä. Time evolution of the quantity (|cos f?|) is shown in Fig. 4. In the simulation with the 3316 molecules (Fig. 4a), (|cos0|) exhibits no large deviation from 0.50. A layer structure, which should be obtained by the use of 18 A cut-oiT, did not appear because of the mismatch in periodicity between the layer structure and the MD cell. On the other hand, in the simulation with the 4736 molecules (Fig. 4b), several large peaks were observed, implying an appearance of a layer structure. Two typical configurations X and Y are shown in Fig. 5. A layer pattern is clearly seen in Y, but not in X. The increase and decrease of (|cos0|) in Fig. 4b indicate the growth and decline of the layer pattern. A complete transition to the layer structure was not realized during the simulation time ~1 ns. Probably, a small mismatch in the periodicity, which can hinder the transition, also exists in this case.
In order to check the reproducibility of the above results, we performed additional 18 A cut-off simulations using a different initial molecular configuration and another software Gromacs3.2 [32,33]. There is no fundamental difference between the simulation algorithms adopted in the Gibbs and the Gromacs3.2. The change of these conditions did not change our results, and the same artificial layer structure appeared when a proper
Fig. 3. Structures of bulk water resulting from the simulations with 2201 water molecules: (a) 18 A cut-off; (b) 9 A cut-off; (c) PME. The black spheres are oxygen atoms, and the grey spheres are hydrogen atoms.
system size was chosen. Therefore, our results do not depend on the initial configuration or software used.
The CPU time used for each calculation is listed in Table 1. As the cut-off length becomes longer, the CPU time increases drastically. The use of the cut-off length 18 A requires a higher computational cost than the PME. This is the generally observed tendency. Thus, when such a long cut-off length is used, in most cases a modification of the conventional cut-off scheme has been made to reduce the computing time. The twin-range method [34] is very useful for computing the
52
Y. Yonetani I Chemical Physics Letters 406 (2005) 49-53
Time [ps]
Fig. 5. Structures of bulk water encountered in the 18 A cut-off Fig. 4. Time evolution of <|cos0|> in the 18 A cut-off simulations: the simulation with 4736 water molecules. X and Y are the structures at
number of water molecules is (a) 3316 and (b) 4736. Structures at the the time indicated by the arrows in Fig. 4.
arrowed points are given in Fig. 5.
long-range contribution to interaction forces fast. In this method, two different cut-off lengths rcl and rc2 {rc\ < r&) are used. Interaction force of molecule-pairs whose distance is less than rcl is updated every time step, whereas interaction force of molecule-pairs whose distance is between rcl and rc2 is updated at a longer time-interval. Although this twin-range method can save a large amount of computing time, it cannot avoid the influence of the cut-off. In fact, this was confirmed by an additional MD simulation we performed using the twin-range method with rcl = 9 A and rc2= 18 A. The contribution from the outer-range between rcl and rc2 was updated every 10 steps. As shown in Table 1, the computing time was largely reduced; the CPU time in the twin-range calculation is about four times shorter than that in the conventional 18 A cut-off calculation. However, the resulting structure was not changed by adopting the twin-range method, that is, the artificial layer structure was obtained.
As a severe artifact caused by cut-off, it is well-known that in radial distribution functions of a solution system containing ions or charged molecules, a remarkable peak appears at a distance corresponding to the cut-off length [9]. On the other hand, as for systems composed of neutral molecules with dipole moment such as water,
Table 1
CPU time for the 1050 ps simulation of a bulk water system of 2201 molecules
CPU time (h)a	
9 A cut-off	26.9
18 A cut-off	173.1
PME	106.3
9 A/18 A twin-range cut-off	42.6
a Measured by a 1.3 GHz Itanium 2 processor.
it has been pointed out that the orientational order of the molecules is enhanced by cut-off treatment [5-7]. The latter artifact is not as severe as the former one (i.e., appearance of the artificial peak), and it is believed that the artificially enhanced order can be reduced by using a longer cut-off length [22]. Thus, many attempts [16,18-20,25] have been done to increase the accuracy of simulation results using the long cut-off length, and then it has been thought that the cut-off length about 18 A is needed for simulating biomolecular systems containing water [22,23,25]. However, it is not assured that such a system leads to a correct result, even if the system does not exhibit strange behavior explicitly. This is because bulk-like water molecules around the solute have a possibility of making an artificial layer structure. A
Y. Yonetani I Chemical Physics Letters 406 (2005) 49-53
53
possible reason that the artificial structure has not been observed for the previous studies is that the number of water molecules around the solute was not enough to form the artificial layer. If we introduce a larger number of water molecules so as to obtain a more realistic result, by contraries, the artificial water layer may appear.
The phenomenon we reported here is the most severe artifact among various cases which have ever been reported for the systems composed of only neutral molecules. This result shows that the use of cut-off can lead to a serious problem, even if the system does not contain ions or charged molecules. Use of the increased cut-off length cannot always improve results of simulations. Moreover, the use of the long cut-off can give a spurious result.
Acknowledgments
The author thanks to Dr. H. Kono for helpful discussions, Professor N. Go for valuable comments, and Dr. K. Yura for constant encouragement. This work was supported by Protein 3000 Project.
References
[1] C.L. Brooks III, B.M. Pettitt, M. Karplus, J. Chem. Phys. 83 (1985) 5897.
[2] HE. Alper, R.M. Levy, J. Chem. Phys. 91 (1989) 1242. [3] K. Tasaki, S. McDonald, J.W. Brady, J. Comput. Chem. 14 (1993) 278.
[4] S.E. Feller, R.W. Pastor, A. Rojnuckarin, S. Bogusz, BR.
Brooks, J. Phys. Chem. 100 (1996) 17011. [5] D. van der Spoel, P.J. van Maaren, H.J.C. Berendsen, J. Chem.
Phys. 108 (1998) 10220. [6] T.M. Nymand, P. Linse, J. Chem. Phys. 112 (2000) 6386. [7] P. Mark, L. Nilsson, J. Comput. Chem. 23 (2002) 1211. [8] J.D. Madura, B.M. Pettitt, Chem. Phys. Lett. 150 (1988) 105. [9] P. Auffinger, DL. Beveridge, Chem. Phys. Lett. 234 (1995) 413. [10] H. Schreiber, O. Steinhauser, Biochemistry 31 (1992) 5856.
[11] R.J. Loncharich, BR. Brooks, Proteins 6 (1989) 32.
[12] M. Saito, Mol. Simul. 8 (1992) 321.
[13] R. Garemyr, A. Elofsson, Proteins 37 (1999) 417.
[14] J. Norberg, L. Nilsson, Biophys. J. 79 (2000) 1537.
[15] M. Patra, M. Karttunen, M.T. Hyvönen, E. Falck, P. Lindqvist,
I. Vattulainen, Biophys. J. 84 (2003) 3636. [16] J. Higo, H. Kono, H. Nakamura, A. Sarai, Proteins 40 (2000)
193.
[17] D. Roccatano, A.E. Mark, S. Hayward, J. Mol. Biol. 310 (2001) 1039.
[18] LH. Shrivastava, D.P. Tieleman, P.C. Biggin, M.S.P. Sansom,
Biophys. J. 83 (2002) 633. [19] CM. Shepherd, H.J. Vogel, D.P. Tieleman, Biochem. J. 370
(2003) 233.
[20] E. Mombelli, R. Morris, W. Taylor, F. Fraternali, Biophys. J. 84 (2003) 1507.
[21] CM. Soares, V.H. Teixeira, A.M. Baptista, Biophys. J. 84 (2003) 1628.
[22] V. Lounnas, in: M.-C. Bellissent-Funel (Ed.), Hydration Process in Biology: Theoretical and Experimental Approaches, IOS Press, Amsterdam, 1999, p. 261.
[23] Y. Komeiji, M. Uebayasi, J. Someya, I. Yamato, Protein Eng. 4 (1991) 871.
[24] D.M. York, T.A. Darden, L.G. Pedersen, J. Chem. Phys. 99 (1993) 8345.
[25] G.E. Arnold, R.L. Ornstein, Proteins 18 (1994) 19.
[26] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey,
M.L. Klein, J. Chem. Phys. 79 (1983) 926. [27] J.-P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys.
23 (1977) 327.
[28] DA. Case, DA. Pearlman, J.W. Caldwell, T.E. Cheatham III, J.
Wang, W.S. Ross, C.L. Simmerling, T.A. Darden, K.M. Merz,
R.V. Stanton, A.L. Cheng, J.J. Vincent, M. Crowley, V. Tsui,
H. Gohlke, R.J. Radmer, Y. Duan, J. Pitera, I. Massova, G.L.
Seibel, U.C  Singh, P.K. Weiner, P.A. Kollman, Amber7,
University of California, San Francisco, 2002. [29] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A.
DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [30] A.R. Leach, Molecular Modeling: Principles and Applications,
Pearson Education Limited, Harlow, England, 2001, p. 327. [31] T. Darden, D. York, L. Pedersen, J. Chem. Phys. 98 (1993) 10089. [32] E. Lindahl, B. Hess, D. van der Spoel, J. Mol. Mod. 7 (2001) 306. [33] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comp. Phys.
Comm. 91 (1995) 43. [34] J. de Vlieg, H.J.C. Berendsen, W.F. van Gunsteren, Proteins 6
(1989) 104.