A LAMMPS implementation of volume–temperature replica exchange molecular dynamics

Computer Physics Communications 189 (2015) 119–127

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

A LAMMPS implementation of volume–temperature replica exchange molecular dynamics✩ Liang-Chun Liu ∗ , Jer-Lai Kuo Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei, 10607, Taiwan

article

info

Article history: Received 17 December 2013 Received in revised form 18 November 2014 Accepted 23 November 2014 Available online 18 December 2014 Keywords: Volume and temperature replica exchange Molecular dynamics LAMMPS

abstract A driver module for executing volume–temperature replica exchange molecular dynamics (VTREMD) was developed for the LAMMPS package. As a patch code, the VTREMD module performs classical molecular dynamics (MD) with Monte Carlo (MC) decisions between MD runs. The goal of inserting the MC step was to increase the breadth of sampled configurational space. In this method, states receive better sampling by making temperature or density swaps with their neighboring states. As an accelerated sampling method, VTREMD is particularly useful to explore states at low temperatures, where systems are easily trapped in local potential wells. As functional examples, TIP4P/Ew and TIP4P/2005 water models were analyzed using VTREMD. The phase diagram in this study covered the deeply supercooled regime, and this test served as a suitable demonstration of the usefulness of VTREMD in overcoming the slow dynamics problem. To facilitate using the current code, attention was also paid on how to optimize the exchange efficiency by using grid allocation. VTREMD was useful for studying systems with rough energy landscapes, such as those with numerous local minima or multiple characteristic time scales. Program summary Program title: vttemper Catalogue identifier: AEVB_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEVB_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: GNU General Public License, version 3 No. of lines in distributed program, including test data, etc.: 49706 No. of bytes in distributed program, including test data, etc.: 1249424 Distribution format: tar.gz Programming language: C++/MPI. Computer: Tested on Intel X86 and AMD64 architectures. Should be operational on any computer with a C++ compiler. Operating system: Tested on Linux. Should run on any platform with C++ and MPI library. Has the code been vectorized or parallelized?: Yes. 1 to N processors may be used. RAM: Depends on the system size and how the program is partitioned. Classification: 16.13. External routines: LAMMPS (http://lammps.sandia.gov) Nature of problem: The code implements volume–temperature replica exchange molecular dynamics for LAMMPS. Each replica has its own temperature and density. A two-dimensional temperature–density grid is constructed.

✩ This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). ∗ Corresponding author. E-mail addresses: [email protected] (L.-C. Liu), [email protected] (J.-L. Kuo).

http://dx.doi.org/10.1016/j.cpc.2014.11.021 0010-4655/© 2014 Elsevier B.V. All rights reserved.

120

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

Solution method: Extending the ability of the existing temper command in LAMMPS to handle replicas of two state variables. Restrictions: The code is based on LAMMPS. A well-designed grid on temperature and density is required. Running time: Running time depends on the system size and the complexity of the problem. Using 8 processors the test run takes approximately 1 minute. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Atomistic-level classical molecular dynamics (MD) is widely used in calculating the dynamic, thermodynamic, and transport properties of numerous physical and biological systems [1]. Although advances in high-performance computing have enabled simulations of micrometers and microseconds [2], for complex materials such as water, crucial physical processes occur on the time scales that are still not attainable by using the conventional numerical methods. MD exhibits computational limitations when it is used to investigate systems with rough energy landscapes. For example, when a rich metastable phase behavior is hidden in the deeply supercooled regime, preventing the system from being trapped in one of the local minima is a challenging task, especially when the thermal energy is too low compared with the energy barriers [3]. Various techniques, of which the objective is to encourage the system to sample a wider configuration space, have been proposed to address this problem. The replica exchange (also called parallel tempering) is one of the most commonly used techniques [4]. In this method, independent runs of M replicas of different temperatures in the canonical ensemble are performed in parallel, and state swaps are attempted at every time interval. The Metropolis criterion is used between replica pairs during the state swap process [5]. This method is designed to enable low-temperature replicas to receive more extensive sampling by exchanging configurations with their high-temperature counterparts. In addition, this method is computationally economical because performing M replicas in parallel accelerates the convergence of all states in a simulation [4]. Furthermore, all thermodynamic properties collected from the replicas are valid throughout the simulation; this property is desirable in constructing a wideranged phase diagram. In this study, we implemented a volume–temperature replica exchange molecular dynamics (VTREMD) method for the classical molecular dynamics package LAMMPS [6]. LAMMPS was developed at Sandia National Laboratories and is freely distributed under the GNU general public license. Because it is highly parallelized through a message passing interface, LAMMPS can handle very large-scale systems. VTREMD is an extension of the existing temper command in LAMMPS, and it is able to perform both temperature and density swaps. With minimum modification to the main body of the source code, this driver program is designed as a command in LAMMPS and can be easily added to the package. In addition, multiple cores can be used by each replica, enabling large-scale MD simulation. Thus, a wide-ranged VTREMD simulation is realized. The rest of this paper is organized as follows. Section 2 presents a brief review of the theoretical background of VTREMD and a description of its implementation in LAMMPS. Section 3 discusses the methods used to perform test runs for TIP4P/Ew [7] and TIP4P/2005 [8]. Analysis of the potential energy (PE) landscape was also presented to determine a favorable grid allocation. In addition, this section describes the test results of two grid numbers used to evaluate the effect of the size of the phase diagram on system convergence. Section 4 presents a brief summary.

2. Theoretical background 2.1. Replica exchange Replica exchange molecular dynamics (REMD) is an accelerated MD method that is particularly useful in addressing systems with rough energy landscapes. The purpose of REMD is to encourage a system to sample more states by inserting a Monte Carlo (MC) decision between regular MD moves [9]. Thus, instead of passively waiting for the occurrence of spontaneous rare events, the system can be manually coaxed toward a direction of lower energy. In this approach, replicas at adjacent temperatures exchange their atomic configurations, and a wider range of transition sampling is therefore realized. Unlike other methods, which boost infrequent events by raising the system temperature, the REMD process retains the system’s thermodynamic information because this swap satisfies the detailed balance condition. In conventional REMD, only temperature is used as the control parameter in determining an exchange move. The proposed VTREMD method, which is adapted from Paschek’s method [10,11], is an extension of traditional REMD [9], in which replicas swap their densities and temperatures. The proposed method samples states more efficiently than traditional REMD because an additional coordinate is added to the MC decision. In this implementation, M by N replicas are activated in parallel, and a two-dimensional grid of M temperatures and N densities is constructed [10,11]. After a certain time interval, replica pairs are chosen to make swap attempts of their temperatures or volumes. In principle, a replica can attempt a T - or V -swap with any member in the M by N system. However, because a successful trial is almost impossible for two replicas with very different states, we can simplify the process by considering only the swap between neighboring pairs. The swap criterion between replicas i and j is determined by using the Metropolis MC algorithm. The acceptance probability of a VTREMD attempt between replicas i and j is given by Eq. (1) [9–11]: Pacc = min 1, exp βi Ui,i − Uj,i + βj Uj,j − Ui,j



 







,

(1)

where Pacc denotes the acceptance rate, βi = 1/kB Ti , and Ti represents the target temperature of state i. Moreover, U represents the instantaneous PE; the first subscript (j) indicates the system configuration, and the second subscript (i) denotes the state density. That is, Uj,i represents the resulting PE of configuration j when its density is scaled to that of state i. To simplify the computation, Uj,i is approximated as a first-order expansion in volume:

 Uj,i ≈ Uj,j −

Pj −

A Vj βj



Vi − Vj ,





(2)

where P denotes the instantaneous pressure, A represents the molecular number in the box, and V is the box size. Because this relation is a first-order approximation, the difference in density between adjacent density grids cannot be too high. After a successful T -swap, temperatures of the replica pairs are exchanged by using the velocity rescaling method [9]. Similarly, replica pairs exchange

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

121

their densities (box size) if a V -swap is invoked, and atoms are then moved to the rescaled location in the new box. For the current rigid water systems, only the intermolecular distances are rescaled, and intramolecular distances and angles are kept constant. 2.2. Application of VTREMD to LAMMPS LAMMPS addresses the replica exchange process in the following manner. First, ‘‘replica’’ and ‘‘state’’ must be defined. The term replica refers to a ‘‘container’’ of atoms, and the term state indicates the thermodynamic properties that a group of atoms have. In other REMD engines, atoms are precisely moved to a corresponding replica during a swap action. Hence, the states are constant for all replicas throughout the simulation. In this case, the terms state and replica are identical. This method is straightforward, and the thermodynamic properties can be easily extracted. However, the method may be problematic when the number of particles is high because moving atoms to another replica invokes processes such as memory allocation and grid initiation. To reduce the computational overhead, LAMMPS adopts another approach by retaining atoms in a replica and changing only the state labels. After a successful swap, the temperature and density of the replica are modified to conform to the new state. A swap table is then generated to record accordingly the updated state of each replica. Fig. 1 shows a schematic of the operating procedure of the proposed VTREMD method. In the beginning of the simulation, a matrix of states, which are aligned according to their densities and temperatures, must be constructed. The states are labeled with an array of numbers. Fig. 1(b) and (c) shows the matrix after a sample temperature and volume swap, respectively. In the LAMMPS implementation, the two-dimensional matrix collapses into a onedimensional array (Fig. 2(a)). The arrows indicate possible swap pathways. Fig. 2(b) shows an example LAMMPS log file of replica swapping. In practice, the boundary replicas receive only half of the chance to exchange with their neighbors; therefore, they may be discarded when thermodynamic properties are computed. 2.3. Using the code 2.3.1. Input The proposed VTREMD is designed as a command in LAMMPS and functions only with the canonical ensemble. Possible thermostats are Nosé–Hoover, Berendsen, Langevin, and direct temperature rescaling [6]. Initial configurations of each replica must be carefully prepared according to the density and temperature grid designs. The command line ‘‘vttemper’’ must be added to initialize the VTREMD process. Its syntax is as follows:

Fig. 1. Schematic example of the evolution of state matrix under a VTREMD procedure. (a) Initial state matrix of the replicas. (b) State matrix after one temperature swap. Replicas of the same density are paired up for a T -swap procedure. In the current methodology, swap attempts are made between the closest neighbors. The pair groups are split into odd and even pairs. The choice of a random selection of pair group or an alternative selection is made by the user. In this example, odd pairs of T -swap are chosen, and successful swaps are highlighted in red. Hence, in this example pair (0, l) experiences a successful swap and their new temperatures are introduced by rescaling the velocities of atoms. (c) State matrix after one volume swap. Successful swaps are highlighted in blue. The densities of swapped replicas are then scaled to their current states. Hence, the matrix table will look very different from the initial one after several steps of VTREMD. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

lis criterion. The term O denotes the number of molecules in the system. 2.3.2. Output In addition to the thermodynamic, dynamic, and configurational information that can be stored simultaneously with the simulation, a swap log file is generated to reflect the current state of each replica (as explained in Section 2.2). This information is required for recovering thermodynamic outputs and restart files because those results are recorded by replicas but not states.

vttemper NM NT NV NP temp fix-ID seed1 seed2 O NM represents the total number of time steps. NT and NV mean that one attempts a T -swap every NT time step and attempts a V -swap every NV T -swap. For example, (NM , NT , NV ) = (500, 10, 10) denotes that the simulation runs for 500 steps, a T swap attempt is invoked every 10 steps, and a V -swap attempt is invoked every 100 steps. A total of 45 T -swap and 5 V -swap attempts are registered. In addition, the condition NM ≥ NT ≥ NV must be satisfied. The number of T -swaps and V -swaps represents the quotient values of NM/NT and NT/NV. NP indicates the size of the temperature grid, and this number must be a factor of the total grid number. That is, if NP equals the number of total grids, then V -swap is not used, and the simulation is simply an REMD simulation. The target temperature and thermostat used in the simulation are temp and fix-ID, respectively. seed1 and seed2 are the random number seeds used in the simulation; seed1 determines the pairing mode, and seed2 represents the random number seed for the Boltzmann weighted Metropo-

2.4. Optimal allocation of temperature and volume grids To achieve an optimal replica exchange, which entails performing a uniform flux of swap rates among different configurations, proper temperature and volume grid allocations must be carefully selected. Rathore et al. [12] and Earl and Deem [4] reported that in the canonical ensemble, the acceptance rate is proportional to the overlap area between the PE distributions of replica pairs. Therefore, an optimal temperature grid allocation can be designed using Eq. (3):

 1E = erfc √ , 2 2σm 

Aov erlap

(3)

where Aov erlap denotes the overlap area of the PE distribution between replica pairs, 1E represents the difference of the mean PE of the distributions, and σm is the average value of the two distributions’ standard deviations. Grid allocation can be fine-tuned

122

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

Fig. 2. (a) Schematic of replica swap in LAMMPS. In this example, replica sets (6 to 11), (12 to 17), and (18 to 23) are grouped since the elements are of the same density. Temperature grids are repeated every 6 replicas. Temperature swap is performed between adjacent replicas. (b) An example swap log table of a VTREMD run in LAMMPS. The entire evolution of the matrix table is needed for post-processing to compute the thermodynamic properties of each ensemble.

by iteratively searching for the next temperature or volume to ensure that the overlapping area is the same as that of the previous pair [4]. In practice, because the PE distribution is not known a priori, a test run can be performed to collect a rough PE landscape and, accordingly, estimate a proper grid allocation. A formal VTREMD simulation can then be designed according to the aforementioned information. However, the PE landscape may evolve continuously as the simulation proceeds, and the grid allocation must be ‘‘massaged’’ occasionally to reach a uniform acceptance rate. 3. Applications 3.1. Example 1: TIP4P/Ew water 3.1.1. Problem setup Water exhibits numerous property anomalies [13], and most of these anomalies emerge in the supercooled or deeply supercooled regimes. However, experimental investigation of supercooled water is impossible because water crystallizes before the system enters the temperature range of interest [14]. A numerical study of super- and deeply supercooled water is therefore critical, and we adopted such a study as a working example for the proposed VTREMD method. In this section, an analysis of TIP4P/Ew [7] water is first presented. Two types of volume–temperature grid designs were employed; one grid comprised 576 states (24 temperatures and 24 densities) and the second grid comprised 720 states (24 temperatures and 30 densities). 512 water molecules were in each replica. Periodic boundary conditions were used in all directions to model the bulk condition. The time step was 2 fs. Temperature was controlled through NVT Nosé-Hoover thermostat with damping parameter as 0.1 ps. The SHAKE algorithm was used to fix the internal OH bonds and H–O–H angles of the water. Lennard-Jones forces were discontinued after 9 Å, and its long-range effect was corrected for energy and pressure [1]. Short-range electrostatic interaction was evaluated using a cutoff distance of 9 Å, whereas long-range force was calculated using the PPPM method [1], and

the accuracy was set as 10−4 . For the case with a large phase diagram, the temperature ranged from 145 to 362.8 K and the density ranged from 0.865 to 1.3g/cm3 . Density grids were uniformly separated every 0.015 g/cm3 . Initial temperature grids were allocated according to the recommendations of Paschek et al. [11]. Replicas were labeled from state 0 to 719, and their initial thermodynamic properties were allocated as described in Section 2.2. Replicas 0 to 23, for example, had the same density (0.865 g/cm3 ) and different temperatures (145 and 362.8 K). The same setting was applied to other replica sets. The case with a small phase diagram was also processed using these settings, except that the lower limit for density was set to 0.955 g/cm3 . Water was prepared from ice Ic through a conventional melt–quench procedure; the cubic phase of ice was heated from 1 to 600 K, equilibrated at 600 K, cooled to the desired temperature, and then equilibrated at the target temperatures. The duration of each stage was 1 ns. After this preparation, a VTREMD run was initiated with T -swap attempts every 100 steps and V -swap attempts every 1000 steps. After a successful V -swap, water molecules were relocated according to the new volume. Both TIP4P/Ew and TIP4P/2005 are rigid water models; therefore, only the intermolecular distances were changed, and the OH bond length and H–O–H angles were kept constant. The molecular topology of the replica was thus kept the same as that before V -swap. 3.1.2. State swap in the current TIP4P/Ew water Fig. 3 shows a sample swap history for visualizing how state swapping is performed. The time evolutions of states in four replicas (Replicas 216, 222, 227, and 235) are illustrated in the figure. The initial density of the replicas was 1 g/cm3 , and the initial temperatures were 145, 183.3, 224.6, and 306.4 K, respectively. As described in Section 2.2, the ‘‘state’’ within a ‘‘replica’’ continuously changes along with a VTREMD simulation; therefore, we monitored the instantaneous state numbers within the aforementioned four replicas. In Fig. 3(a), a steep jump implies a change in the density state, and minor variations indicate temperature swaps. Fig. 3(b) shows a longer time evolution of the same test case. The replica with the lowest initial temperature (Replica 216)

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

123

Fig. 3. (a) A sample swap history in the TIP4P/Ew simulation of four replicas. Their initial temperatures are 145 K, 183 K, 224.6 K and 306.4 K and initial density is 1 g/cm3 . State number in the y-axis indicates the current state of each replica after multiple swaps. (b) A longer time evolution profile of swap history.

had a limited span of state swapping with its neighbors compared with other replicas. In ideal cases, low-temperature states receive enhanced sampling because swap is attempted with hightemperature states, and such swaps must propagate until all configurations are sampled. However, a high acceptance rate of low temperature states does not guarantee a good sampling to the high-temperature ones. Instead, state swaps may be limited to happen only with their neighbors. Hence, an intermediate region may be necessary to facilitate information exchange. In the current case, Replica 227 had the widest sampling among the states, indicating that it assumes the role of communicating the configurational sampling between high- and low-temperature states. 3.1.3. Design of grid allocation Fig. 4(a) illustrates the selected T -swap acceptance rates collected from 8 to 14 ns of the TIP4P/Ew simulation based on the original temperature grids [11], indicating that states at lower temperatures exhibited higher swap successful rates than those at high temperatures did. Although the difference was not large, the spacing of the temperature grid can be fine-tuned to unify the acceptance rates of replica pairs with different temperatures. Because the proposed method requires a matrix of temperature–density space, the same grid design must be applied to the whole system; hence 1 g/cm3 was used as the model density. Fig. 5(a) shows the PE landscape of the 1 g/cm3 cases collected between 8–14 ns. With this information, the PE is curve-fitted with Gaussian distributions; Fig. 5(b) and (c) shows the average PE and its standard deviation, respectively. The estimated acceptance rate can thus be smoothed by rearranging the temperature grid according to Fig. 5 and Eq. (3).

Fig. 4. (a) T -swap acceptance rate of selected densities with the temperature grids adopted from Paschek et al. [10]. (b) T -swap acceptance rate after grid refinement.

Table 1 lists the original [10] and revised temperature grids. With this modification, VTREMD simulation was performed for another 5 ns by using the new grid allocation. Fig. 4(b) shows the acceptance rate of the new grid design. A more uniform distribution of acceptance rate at 1 g/cm3 , approximately 0.06, was then obtained. The aforementioned protocol can also be used for designing density grids. However, because the potential landscapes of different densities at the same temperature are closely packed, manually allocating the density grids is easier than invoking the discussed protocol. 3.1.4. Region of temperature grid VTREMD is much more expensive than conventional REMD because it is calculated on two coordinates, and the computational cost increases linearly as the size of the temperature–density grid increases. Theoretically, replicas at low temperatures must swap with those at high temperatures to facilitate their escape from PE basins. However, when the high-temperature part has converged, a full-domain simulation may not be required. Therefore, we adopted an approach that involves cascading the computational grids (i.e., stop simulating the converged regime and perform VTREMD only in the grids that are still in transition). The objective of the proposed VTREMD method is to enhance space sampling; therefore, the molecules in the upper limit of temperature should

124

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

Fig. 5. (a) Potential energy probability distribution functions of the 24 temperature grids at 1 g/cm3 . (b) Averaged potential energy and (c) its corresponding standard deviation σ of 1 g/cm3 TIP4P/Ew water. Temperatures span from 145 K to 363 K, and the statistics is collected every 20 fs for 6 ns at the canonical ensemble. Table 1 Original and modified temperature grids in Section 3.1.3. Original Modified

145 145

150 152

156 159

162.4 166.5

169 174.5

176 182.5

183.3 190.5

191 198.5

199 206.5

207.3 214.5

215.8 222.5

224.6 230.5

Original Modified

233.5 238.5

242.7 247

252.1 255.5

261.9 264

272.2 273

283 282.5

294.4 292.5

306.4 302.5

319.2 313.5

332.9 325

347.7 337.5

362.8 350.5

be sufficiently energetic to facilitate good state exploration. As a demonstration, Fig. 6 shows plots of the mean square displacement (MSD) of oxygen atoms with a density of 0.955 g/cm3 and different temperatures. The MSD plots indicate that the time required for systems to enter the diffusive regime increases exponentially as the temperature decreases. For the 261.9 K case, for example, the water in this particular state requires a simulation time of up to 300 ps to migrate across three molecules (approximately 1 nm). This simulation time is computationally affordable and thus taken as the upper limit for the subsequent VTREMD test. Physically, because thermodynamic equilibrium cannot be achieved before the system is dynamically relaxed, the self-diffusion behavior can be used as a reference to determine the upper limit in the temperature grid.

3.1.5. Region of density grid To demonstrate the importance of density grid design, we calculated the TIP4P/Ew phase diagram on the basis of two lower density limits. For the case with a smaller phase-diagram size, a density limit similar to that in [11], which spans from 0.955 to 1.3 g/cm3 every 0.015 g/cm3 , was adopted. Therefore, 24 temperature and 24 density states were obtained (i.e., a total of 576 grid points). Initial structures were derived from preequilibrated SPC/E water through a continual VTREMD simulation up to 30 ns. Fig. 7 shows the temporal behavior of the phase diagrams of TIP4P/Ew water; each diagram was obtained by averaging the pressure profiles over 0.5 ns at t = 0, 7.5, 15, and 22.5 ns, respectively. Isobars are plotted every 50 MPa. The TIP4P/Ew results obtained in [11] are also

Fig. 6. Oxygen mean-square-displacement plots of TIP4P/Ew water with density as 0.955 g/cm3 . Temperatures are 145, 150.0, 156.0, 162.4, 169.0, 176.0, 183.3, 191.0, 199.0, 207.3, 215.8, 224.6, 233.5, 242.7, 252.1, 261.9, 272.2, 283.0, 294.4, 306.4, 319.2, and 332.9 K, from bottom to top.

plotted as a benchmark for our simulation. The circles indicate isobars from [11] at 0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, and 1 GPa, from bottom to top. The solid circles are at a constant pressure of 0 GPa. The transient plots indicate that, despite the low-temperature regime (where slow dynamics occurs), the simulation quickly achieved the same system feature as did that in [11]. The thin black ellipses indicate the region at which the putative liquid–liquid phase transi-

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

125

Fig. 7. Time evolution of TIP4P/Ew phase diagram with lower limit in density at 0.955 g/cm3 . Each phase diagram is averaged over a time frame of 0.5 ns. Isobars are plotted every 50 MPa. Circles indicate VTREMD results by Paschek et al. [11], representing isobars of 0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, and 1 GPa, from bottom to top. Thin ellipses indicate where the slow dynamics emerges, and thick ellipses show the critical region that the two VTREMD tests give slightly different results.

tion occurs; therefore, the states in this region require more time to determine the condition to reside. Furthermore, the results of the test case were highly similar to those of [11], but the test failed to converge to exactly the same pressure in the region depicted by the thick red ellipses. For the isobar at 0 GPa, for example, under the same pressure, the current prediction generated a density that is lower than that of [11] by approximately 0.07 g/cm3 , and this difference persisted even after a long (30 ns) VTREMD simulation. Therefore, we infer that the current test reached its convergence and that the test enables recovering the prediction of Paschek et al. [11] qualitatively. To observe the influence of the density limit of the VTREMD grid on the convergence state, we added two sets of density grid points (0.940 and 0.925 g/cm3 ) below those of the aforementioned test. The highest two temperature grid points were discarded (347.4 and 362.8 K) because the subsequent highest temperature (332.9 K) was sufficiently high to facilitate sampling the phase space of the remaining states. The temperature grids therefore spanned from 145 to 332.9 K (22 temperatures), and the density grids spanned from 0.940 to 1.3 g/cm3 (26 densities), comprising a total of 572 states. To observe the transient behavior, we initiated the VTREMD simulation and recorded its phase diagram for 2 ns, with each phase diagram being averaged at an interval of 0.5 ns. The initial structure of each state was obtained from various snapshots of the test in the previous section. The melting–quenching procedure was used again to prepare the initial condition for the two newly added densities. To examine the influence of the initial structure and reduce the statistical uncertainty, five snapshots (t = 4, 5, 6, 8, and 12 ns) of the discussed test were extracted and used as the initial conditions. We verified that the five VTREMD tests produce the same prediction; therefore, we computed the average of these tests. Fig. 8 shows the phase diagrams with this enlarged phase space. Because the new test inherited the structure of the previous test, we assumed that they behave similarly, at least at the early stages. However, the low-temperature region of the isobar at 0

GPa indicated that the results of this study’s calculation converge within 0.5 ns with respect to the results obtained in [11]; this behavior was not observed in the previous long VTREMD test of a smaller phase space. In addition, the calculated isobar at 0 GPa seemed to move ‘‘downward’’ as the time increased. For the 183 K and 0.97 g/cm3 state, for example, the calculated average pressure before the two lower-density grids were included was −31 MPa (Fig. 7). The average pressure of this state soon increased to 24 MPa (Fig. 8, 0–0.5 ns) and then to 66 MPa (Fig. 8, 1.5–2 ns) after the two lower density grids were included. Therefore, for this particular case, a density close to 0.955 g/cm3 and temperature below 215 K were critical. The added density grids may accelerate the convergence by adjusting the configuration of the jammed states. Finally, we compared two VTREMD test results of TIP4P/Ew with different phase-diagram domain sizes. Both tests had the same temperature grids as those of Paschek et al. [11], and the lower limit for density was 0.955 g/cm3 for the small grid and 0.865 g/cm3 for the large grid. The large grid had a lower density limit than that shown in Fig. 8. Fig. 9(a) and (b) illustrates the two cases. The results obtained in [11] and the experimental measurements in [15] were plotted for comparison. As shown in Fig. 9(b), along the 0-GPa isobar and below 215 K, a clear difference was observed between the test results of the proposed VTREMD method and those obtained in the literature. The adjacent isobars were further pushed ‘‘downward’’ compared with those illustrated in Fig. 8. The change in the phase diagram was phenomenal, and the results of this test suggested that density swapping facilitates state exploration and indicated the existence of a pathway that ‘‘loosens up’’ the jammed states caused by the state exchange in the density coordinates. We closely observed the PE profiles of selected states. Fig. 10 illustrates the plotted PE values of the 0.955 g/cm3 states shown in Fig. 9(a) and (b). The PE values were averaged over 0.5 ns after equilibrium. This density was used because it is the ‘‘boundary’’ density of the smaller case, where a V-swap is only possible with higher densities. It is seen that states of a

126

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

Fig. 8. Time evolution of TIP4P/Ew phase diagram with lower limit in density at 0.865 g/cm3 . Isobars are plotted every 50 MPa. Circles indicate VTREMD results by Paschek et al. [11], representing isobars of 0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, and 1 GPa, from bottom to top. The red ellipses indicate the critical region where states gradually evolve.

Fig. 9. Phase diagram of TIP4P/Ew water obtained by two VTREMD tests of different ranges, where (a) is averaged over 7 ns and (b) is averaged over long time (high temperature for 14 ns and low temperature for 28 ns). White circles indicate VTREMD results by Paschek et al. [11], representing isobars of 0, 100, 200, 300, 400, 600, and 800 MPa from bottom to top. Solid white lines are the isobars at 0.1, 50, 100, 200, 500 MPa and 1 GPa (bottom to top) from experimental measurements [15].

larger grid domain process smaller potential energies (i.e., more thermodynamically stable) than their counterparts calculated in a smaller grid. We believe that the two schemes will correspond with each other, provided that a very long simulation is conducted. With the current limited computational capability, it is demonstrated that the VTREMD test with a larger bounds converges faster than the smaller one. Hence, the possible effect of the phase domain size must be considered in designing VTREMD tests. 3.2. Example 2: TIP4P/2005 water Another popular four-point water model, TIP4P/2005 [8], was also calculated (Fig. 11) to demonstrate the effectiveness of the

proposed method. Except the model parameters and cutoff lengths of Lennard-Jones and Coulomb interaction (both were 9.5 Å), all other settings were the same as those used in the TIP4P/Ew case. A VTREMD test was first performed across the entire phase domain. Average pressures were computed for 8.74 ns after convergence. Because the high-temperature regime had converged, the temperature range was reduced (145 K to 262 K) to decrease the computational cost. Another 20-ns VTREMD simulation of this lowtemperature regime was performed. The phase diagram shown in Fig. 11 represents the average value of the two-stage simulation, and it is consistent with the single-run MD results of Pi et al. [16] (isobars at 0.1, 40, 100, and 150 MPa). For the 0.1-MPa isobar, for example, Pi et al. [16] used a method similar to slow cooling.

L.-C. Liu, J.-L. Kuo / Computer Physics Communications 189 (2015) 119–127

127

4. Summary In this study, a VTREMD command for LAMMPS was developed and verified using working examples of two four-site water models. The TIP4P/Ew water example indicated that VTREMD enabled the system to converge within 42 ns, which is much shorter than the duration of classical MD simulations. The thermodynamic properties in the deeply supercooled regime can thus be computed with affordable computational expenses. The derived topologies can also be further investigated to understand their dynamic properties. In this paper, the method for grid design and importance of the size of the phase diagram are described. In conclusion, VTREMD provides an efficient technique for accelerating classical MD simulations and simultaneously covering a large-range phase diagram; it is particularly useful when a system becomes trapped in a glassy state. Acknowledgments Fig. 10. The average potential energy profiles of TIP4P/Ew water at 0.955 g/cm3 , collected from a larger (0.865 to 1.3 g/cm3 ) and smaller (0.955 to 1.3 g/cm3 ) VTREMD range in density. Data is averaged over 0.5 ns of simulation.

This work is supported by Academia Sinica under Research Program on NanoScience and NanoTechnology, the Ministry of Science and Technology of Taiwan (NSC101-2113-M-001-023-MY3 and NSC-101-2911-I-001-501), and the National Center for Theoretical Sciences (South) Physics Division. Computational resources in part are supported by the National Center for High Performance Computing. References

Fig. 11. Phase diagram of TIP4P/2005 water calculated with VTREMD method. White circles indicate the single-run MD results from Tables 1–4 of Pi et al. [16]; isobars are at 0.1, 40, 100, and 150 MPa, from bottom to top. Solid lines are the isobars at 0.1, 50, 100, 200, 500 MPa and 1 GPa (bottom to top) from experimental measurements [15].

The simulations were conducted consecutively from high temperatures to lower temperatures. Each state performs an NPT simulation for 40 ns; therefore, it takes a considerable amount of time to reach the lowest temperature because the NPT simulations must be performed sequentially. Furthermore, because VTREMD performs calculations in parallel, it thus achieves the same accuracy as [16] but with much lower wall-time expenses. The experimental results [15] from NIST were plotted as a benchmark for the model performance. Both water models (TIP4P/Ew and TIP4P/2005) reproduced the experimental values at low densities but deviated from NIST results as the density increased. Of the two models, TIP4P/2005 exhibited better consistency with experiments.

[1] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989. [2] Y. Li, J. Li, F. Wang, Liquid-liquid transition in supercooled water suggested by microsecond simulations, Proc. Natl. Acad. Sci. USA 110 (2013) 12209–12212. http://dx.doi.org/10.1073/pnas.1309042110. [3] F.H. Stillinger, A topographic view of supercooled liquids and glass formation, Science 267 (1995) 1935–1939. http://dx.doi.org/10.1126/science.267.5206. 1935. [4] D.J. Earl, M.W. Deem, Parallel tempering: Theory, applications, and new perspectives, Phys. Chem. Chem. Phys. 7 (2005) 3910–3916. http://dx.doi.org/ 10.1039/B509983H. [5] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (2001). [6] http://lammps.sandia.gov/. [7] H.W. Horn, W.C. Swope, J.W. Pitera, J.D. Madura, T.J. Dick, G.L. Hura, et al., Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew, J. Chem. Phys. 120 (2004) 9665–9678. http://dx.doi.org/10. 1063/1.1683075. [8] J.L.F. Abascal, C. Vega, A general purpose model for the condensed phases of water: TIP4P/2005, J. Chem. Phys. 123 (2005) 234505. http://dx.doi.org/10. 1063/1.2121687. [9] Y. Sugita, Y. Okamoto, Replica-exchange molecular dynamics method for protein folding, Chem. Phys. Lett. 314 (1999) 141–151. http://dx.doi.org/10. 1016/S0009-2614(99)01123-9. [10] D. Paschek, A.E. García, Reversible temperature and pressure denaturation of a protein fragment: a replica exchange molecular dynamics simulation study, Phys. Rev. Lett. 93 (2004) 238105. http://dx.doi.org/10.1103/PhysRevLett.93. 238105. [11] D. Paschek, A. Rüppert, A. Geiger, Thermodynamic and structural characterization of the transformation from a metastable low-density to a very high-density form of supercooled TIP4P-Ew model water, Chem. Phys. Chem. 9 (2008) 2737–2741. http://dx.doi.org/10.1002/cphc.200800539. [12] N. Rathore, M. Chopra, J.J. de Pablo, Optimal allocation of replicas in parallel tempering simulations, J. Chem. Phys. 122 (2004) 024111. http://dx.doi.org/ 10.1063/1.1831273. [13] http://www.lsbu.ac.uk/water/anmlies.html. [14] F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, D. Majolino, V. Venuti, et al., Evidence of the existence of the low-density liquid phase in supercooled, confined water, Proc. Natl. Acad. Sci. 104 (2007) 424–428. http://dx.doi.org/ 10.1073/pnas.0607138104. [15] http://webbook.nist.gov/chemistry/fluid/. [16] H.L. Pi, J.L. Aragones, C. Vega, E.G. Noya, J.L.F. Abascal, M.a. Gonzalez, et al., Anomalies in water as obtained from computer simulations of the TIP4P/2005 model: density maxima, and density, isothermal compressibility and heat capacity minima, Molecular Phys. 107 (2009) 365–374. http://dx.doi.org/10. 1080/00268970902784926.