Accurate Relative Energies and Binding Energies ... - ACS Publications

Apr 17, 2017 - ωB97XD optimized water 50-mers (Figure 1b), the local modes and .... set as 4.0 Å, and λ is chosen as 12 and 8 [6 for CCSD(T)-. F12b...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCA

Accurate Relative Energies and Binding Energies of Large Ice−Liquid Water Clusters and Periodic Structures Lei Zhang, Wei Li,* Tao Fang, and Shuhua Li Institute of Theoretical and Computational Chemistry, Key Laboratory of Mesoscopic Chemistry of MOE, School of Chemistry and Chemical Engineering, Nanjing University, Nanjing 210023, China S Supporting Information *

ABSTRACT: Relative energies and binding energies are crucial quantities that determine various molecular properties of ice and water. We developed a new effective method to compute those energies of bulk ice−liquid water systems. In this work, ten ice− liquid 144-mers and ten periodic ice−liquid (H2O)64 systems are taken from the molecular dynamics simulations in the melting process of ice-Ih crystals. They are investigated at the levels of density functional theory (DFT), explicitly correlated secondorder Møller−Plesset perturbation theory (MP2-F12), and coupled-cluster singles and doubles with noniterative triples corrections [CCSD(T)-F12b] in the framework of generalized energy-based fragmentation approach. Our results show that the changing of noncovalent interactions significantly influences the performances of DFT and electron correlation methods for those systems in the melting process of ice. Various DFT methods predict quite different results for ice and mixed ice−liquid structures but give similar results for pure liquid ones. It also explains why many DFT-based simulations lead to inaccurate densities of ice and liquid water. The CCSD(T)-F12b results suggest that the MP2-F12 method provides satisfactory results and is expected to be employed to simulate the phase transitions of ice crystal.

1. INTRODUCTION The studies of water behaviors, especially the phase transition from liquid to ice, are of great importance for our deep understanding of the complicated structures and interactions in liquid water and ice.1−3 In the past decade, various advanced experimental efforts have been made to shed light on this process. For example, ultrafast X-ray diffraction, low-temperature scanning tunneling microscopy, and glancing-angle Raman spectroscopy have been applied to water droplets, the surface, and the interface of liquid water and ice.4−6 In recent years, the requirement for the molecular level understanding is urgent. The classical molecular dynamics (MD) simulations with force fields (FFs) for phase transitions of water are still great challenges.7−9 The ab initio MD (AIMD) simulations based on density functional theory (DFT) are also difficult to describe the structural and dynamical properties, such as density and melting temperature (Tm), of liquid water and ice.10−12 The simulations based on more accurate electron structure methods have been performed on liquid water and solid−liquid water interface.13−15 For example, the simulations using embedded-fragment second-order Møller−Plesset perturbation theory (MP2) could predict various properties of liquid water.15 In high-level AIMD simulations, the fast and accurate electron structure calculations of large water systems are required. However, even for water clusters, it is still a challenge to accurately predict the relative stabilities or binding energies with ab initio and DFT methods.16−19 Traditionally, the © 2017 American Chemical Society

coupled-cluster singles and doubles with noniterative triples corrections [CCSD(T)] method20 is employed as benchmark for small- and medium-sized water clusters.17,21−27 For larger water clusters, currently, it is impossible to perform conventional CCSD(T) calculations due to the steep scaling. Therefore, new efficient methods are urgently required for accurate describing the relative energies or binding energies of large water systems. In the past decade, fragmentation approaches have been developed to treat the total energy and molecular properties of large systems at various electron structure methods.28−55 By combining the generalized energy-based fragmentation (GEBF) approach developed by our group29,30 and explicitly correlated F12 method,56−58 the GEBF-X-F12 [X = MP2, CCSD, CCSD(T)] method19,59 has been employed for the binding energies and relative energies of large clusters, such as water 20mers (Figure 1a), at the complete basis set (CBS) limit, which provide accurate values. In addition, the GEBF approach is found to be applicable for the calculations of medium-sized water clusters at any basis set.60 Very recently, for four GEBFωB97XD optimized water 50-mers (Figure 1b), the local modes and hydrogen-bonding analyses explained why does warm water freeze faster than cold water.61 In the melting process of ice crystal, due to the changes of hydrogen-bond length and strength, the electron correlation Received: April 10, 2017 Published: April 17, 2017 4030

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A

Figure 1. Structures of (a) 20-mer, (b) 50-mer, and (c) 144-mer water clusters and (d) the periodic (H2O)64 system.

Figure 2. Ten structures of ice−liquid 144-mers.

and RAB is the distance between atoms A and B. For the periodic systems, the total energy per unit cell could be expressed as31

methods are required to treat large water systems under the periodic boundary conditions (PBC). Recently, the GEBF approach under the PBC (PBC-GEBF) has been developed and used for the structures and vibrational spectra of molecular crystals.31,32 Therefore, we employed the MP2-F12 and CCSD(T)-F12b58 methods in the PBC-GEBF calculations. Using the present method, we can extend the medium-sized water clusters to water 144-mers and the periodic (H2O)64 systems (Figure 1c,d) in the melting process of ice crystal. This paper is organized as follows. In section 2, the GEBF and PBC-GEBF approaches are briefly introduced and the computational details are described. In section 3, the GEBF and PBC-GEBF relative energies and binding energies of ice−liquid clusters and periodic ice−liquid systems at the DFT, MP2-F12, and CCSD(T)-F12b levels are compared, respectively. Finally, a brief summary is provided in section 4.

M PBC ‐ GEBF Eunit cell

k

+ E Ewald sum

N

N

∑ CkEk̃ − (∑ Ck − 1) ∑ k

k



A B>A

Q AQ B RAB

k

∑ ∑ A∈m B>A∈m

Q AQ B RAB (2)

where Ẽ k is the energy of the kth subsystem embedded in the field of point charges of a super cell and EEwald sum is the classical charge−charge interaction energies in the central unit (denoted as m) and those between the central and all image cells obtained using the Ewald summation method.62 The basic steps of constructing GEBF (or PBC-GEBF) subsystems29−31 mainly include: (1) Divide the target system (or the unit cell) into various fragments. (2) Construct primitive subsystem (with its coefficient being one) for each fragment (or fragment in the central cell) by adding its neighboring fragments with a distance threshold ζ. The total number of fragments within any primitive subsystem is limited by another parameter λ. (3) Determine derivative subsystems and their coefficients to remove the overcounting of k-fragment (k < λ) terms that arise from the overlapping of primitive subsystems. For periodic systems, the translation symmetry is employed to determine if two k-fragment terms in different cells are identical. (4) Construct additional two-fragment subsystems (if not previously existed) if a fragment−fragment distance is less than 2ζ. (5) Place each subsystem in the field of natural

2. METHODOLOGY AND COMPUTATIONAL DETAILS GEBF and PBC-GEBF Approaches. In the GEBF approach, the total ground-state energy of a targe system can be approximately obtained from the corresponding energies of a series of electrostatically embedded subsystems as29,30 GEBF Etot =

M

= ∑ CkEk̃ − (∑ Ck)

(1)

where Ck and Ẽ k are the coefficient and the energy (including the self-energy of point charges) of the kth embedded subsystem, respectively, QA is the point charge on atom A, 4031

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A

Figure 3. Ten structures of periodic ice−liquid (H2O)64 systems.

charges63,64 on all other atoms in target system (or super cell) to take the long-range electrostatic interaction and polarization effects into account. The GEBF and PBC-GEBF approaches have been implemented in the LSQC program65,66 at the HF, DFT, and MP2 levels for energy calculations, geometry optimizations, and frequency calculations, and at the CCSD, CCSD(T), MP2-F12, and CCSD(T)-F12a levels for energy calculations, of a broad range of large systems and molecular crystals.19,29−32,59,67,68 Computational Details. In this work, we have employed GEBF- and PBC-GEBF-CCSD(T)-F12b methods to benchmark the relative energies (REs) and binding energies (BEs) of large water clusters and periodic water systems. Ten water 144mers (Figure 2) and ten periodic (H2O)64 systems (Figure 3) are taken from the MD simulations of melting ice-Ih crystals with 768 and 64 water molecules, respectively. The simulations were carried in the NPT ensemble at 270 K with the TIP3P force field.69 In each series, structure 1 is (or is taken from) the initial ice-Ih crystal, structures 2−6 (mixed ice−liquid systems) are taken from the melting area, and structures 7−10 (liquid systems) are taken from the equilibrium area. Here structure 6 is chosen as the reference point of the ten systems for the REs because it can be approximately considered as a turning point from ice to liquid. For example, the first six periodic (H2O)64 structures are taken from the first 50 ps and the remainding ones are taken from the remaining trajectory in the 1 ns MD simulation of the ice-Ih (H2O)64 crystal (Figure S1). The GEBF and PBC-GEBF REs and BEs are computed at various DFT, MP2-F12, and CCSD(T)-F12b levels. Those DFT functionals include BLYP,70,71 B3LYP,70,71 and PBE72 (with and without empirical D3 dispersion corrections73), ωB97X,74 ωB97X-D,75 and M06-2X.76 The BE of (H2O)n is defined as EBE = nEwater − Etot, where Etot is the total energy and Ewater is the energy of optimized water monomer. Due to the large values of REs and BEs, we employ the REs per molecule (REPMs) (relative to structure 6) and BEs per molecule (BEPMs) in our analyses. In the GEBF and PBC-GEBF calculations, each water molecule is defined as a fragment, ζ is set as 4.0 Å, and λ is chosen as 12 and 8 [6 for CCSD(T)F12b], respectively. In the DFT calculations, we employ a popular Pople’s triple-ζ basis set, 6-311++G(d,p), which

balances the computational costs and accuracy. In the electron correlation calculations, we chose the correlation consistent ccpVDZ-F12 basis set, which is specially optimized for F12 methods to yield the results near the CBS limit.77 To characterize the noncovalent interactions, we have also performed reduced density gradient (RDG)78,79 analyses for water 144-mers and hydrogen-bonding analyses for all systems. The RDG (dimensionless), defining as s(ρ) = (1/2)(3π2)1/3|∇ρ|/ρ4/3, is usually plotted as the function of sign (λ2)ρ, where ρ is electron density, λ2 is the second eigenvalue of Hessian matrix, and “sign” is symbolic function. It can be used to distinguish the types of weak interactions. In a RDG figure, the left, middle, and right regions qualitatively represent the stabilizing, van der Waals (vdW), and steric interactions, respectively,78 and the low-density, low-gradient spikes at negative values (left, bottom) represent the hydrogen-bonding interactions. In the hydrogen-bonding analyses, a hydrogen bond (HB) (X−H···Y) is assigned if the H···Y length is between 1.5 and 2.8 Å, and the X−H···Y angle is larger than 145.0°. Specifically, a strong HB is defined if the H···Y length is less than 2.0 Å, instead of 2.8 Å. The MD simulations were performed using the TINKER 6.1 package.80 The GEBF and PBC-GEBF calculations were executed by the LSQC program.65,66 The DFT and electron correlation calculations of subsystems were obtained with the GAUSSIAN0981 and MOLPRO 201282,83 packages, respectively. The RDG analyses are performed using the Multiwfn program.84 All calculations were carried out on workstations equipped with Intel Xeon CPUs with two dodeca-core E5-2692 v2 2.20 GHz or E5-2680 v3 2.50 GHz.

3. RESULTS AND DISCUSSION The REs even for small- and medium-sized water clusters are difficult to be predicted. The benchmark CCSD(T) results usually recommend different functionals for water clusters of different sizes. For example, the best functionals are PWB6K,85 MPWB1K,86 and M05-2X87 for water hexamers,17 ωB97X-D75 for 16-mers,18 LC-ωPBE-D388 for 17-mers,18 and LC-ωPBED3 and B97-D89 for water 20-mers,19 respectively. For bulk water, especially those in the melting process, the choices of 4032

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A

Table 1. Comparisons of the Number of Hydrogen Bonds (n), Average Coordinated Numbers (⟨k⟩), and the Number of kCoordinated Water Molecules (k = 0, 1, ..., 5) of the Ten Ice−Liquid 144-mersa number of k-coordinated molecules (k) n

structure 1 2 3 4 5 6 7 8 9 10 a

230 233 233 225 194 182 150 155 162 151

(230) (177) (177) (167) (140) (124) (98) (98) (116) (105)

⟨k⟩ 3.19 3.24 3.24 3.12 2.69 2.53 2.08 2.15 2.25 2.10

0

(3.19) (2.46) (2.46) (2.32) (1.94) (1.72) (1.36) (1.36) (1.61) (1.46)

0 0 0 0 2 5 7 8 6 12

1

(0) (0) (0) (4) (9) (12) (27) (27) (16) (26)

7 6 6 7 17 22 38 29 35 31

2

(7) (22) (22) (24) (42) (52) (61) (53) (56) (49)

28 24 24 31 44 39 46 57 38 53

(28) (54) (54) (56) (51) (47) (37) (49) (44) (48)

3 39 44 44 43 42 48 43 34 48 30

(39) (48) (48) (42) (32) (30) (15) (15) (24) (19)

4 70 70 70 63 38 30 9 15 16 15

(70) (20) (20) (18) (10) (3) (4) (0) (4) (2)

5 0 0 0 0 1 0 1 1 1 3

(0) (0) (0) (0) (0) (0) (0) (0) (0) (0)

The corresponding values for strong hydrogen bonds included in parentheses.

Table 2. Comparisons of the GEBF REPMs of the Ten Ice−Liquid 144-mers at Various DFT Levels with the 6-311++G(d,p) Basis Set and at MP2-F12 and CCSD(T)-F12b Levels with the cc-pVDZ-F12 Basis Seta structure 1 2 3 4 5 6 7 8 9 10 NPEc MUEd

PBE-D3b −4.16 −3.12 −2.88 −1.82 −0.22 0.00 1.03 1.04 0.94 0.87 2.24 2.03

(−4.18) (−3.29) (−3.04) (−1.96) (−0.33) (0.00) (1.08) (1.09) (0.97) (0.93) (2.32) (2.05)

BLYP-D3b −4.02 −3.07 −2.83 −1.82 −0.24 0.00 1.01 1.01 0.92 0.85 2.08 1.89

(−4.05) (−3.37) (−3.03) (−1.99) (−0.37) (0.00) (1.07) (1.07) (0.94) (0.92) (2.18) (1.92)

B3LYP-D3b −3.14 −3.07 −2.83 −1.81 −0.22 0.00 1.03 1.02 0.93 0.86 1.21 1.01

(−3.17) (−3.24) (−3.00) (−1.96) (−0.33) (0.00) (1.08) (1.07) (0.96) (0.92) (1.30) (1.04)

ωB97X-Db −2.68 −3.01 −2.77 −1.76 −0.19 0.00 1.01 0.99 0.91 0.84 0.73 0.55

(−2.65) (−2.97) (−2.74) (−1.75) (−0.18) (0.00) (1.01) (0.95) (0.89) (0.79) (0.66) (0.52)

M06-2X

MP2-F12

CCSD(T)-F12b

−2.29 −2.68 −2.46 −1.57 −0.10 0.00 0.91 0.85 0.80 0.67 0.23 0.16

−2.45 −2.81 −2.57 −1.60 −0.18 0.00 0.90 0.87 0.84 0.69 0.36 0.32

−2.13 −2.74 −2.51 −1.56 −0.17 0.00 0.87 0.83 0.80 0.66 0.00 0.00

a

All energies are in kcal/mol. bThe value without the empirical dispersion correction included in parentheses. cNonparallelity error relative to the CCSD(T)-F12b results. dMaximum unsigned error relative to the CCSD(T)-F12b results.

Relative Energies of Ice−Liquid 144-mers. We have compared GEBF REPMs of the ten 144-mers at PBE-D3, BLYP-D3, B3LYP-D3, ωB97XD, MP2-F12, M06-2X, and CCSD(T)-F12b levels (Table 2 and Figure 4). Table 2 shows that M06-2X, MP2-F12, ωB97X, and ωB97XD methods provide satisfactory or reasonable REPMs with maximum unsigned errors (MUEs) being only 0.16, 0.32, 0.52, and 0.55 kcal/mol, respectively. The MUEs for B3LYP-D3, BLYP-D3,

functionals or other ab initio methods presents a greater challenge. Hydrogen-Bonding and RDG Analyses of Ice−Liquid 144-mers. The hydrogen-bonding analyses of the ten ice− liquid 144-mers are listed in Table 1. That shows the numbers of HBs decrease from 233 to 183 in the ice-Ih and mixed ice− liquid structures 1−6 but remain around 155 in the liquid structures 7−10. Those numbers in the first three structures are similar, ranging from 230 to 233. However, all the 230 HBs in ice-Ih structure are strong, whereas only three-quarters of HBs in the structure 2 and 3 are strong. Due to surfaces, the average coordinated numbers of ice-Ih structure (3.2) are less than the expected value (4.0). Only 20 of 70 interior (four-coordinated) water molecules are still coordinated by strong HBs in structures 2 and 3. In the remaining structures, more 1,2,3coordinated water molecules spring up. In the liquid structures 7−10, the k-coordinated (k = 1−5) water molecules for each k are similar. Thus, the decreasing (strong) HBs and average coordinated numbers in the ice-Ih and mixed ice−liquid structures reduces the hydrogen-bonding interactions. The RDG analyses in Figures S2 and S3 also indicate the decreasing of hydrogen-bonding interactions in the ice-Ih and mixed ice− liquid structures (through the increasing of spike values) and the similar hydrogen-bonding interactions in the liquid ones. Additionally, they also suggest that both the vdW and steric interactions also change gradually in structures 1−6 but remain similar in structures 7−10.

Figure 4. Comparisons of the GEBF REPMs of the ten ice−liquid 144-mers at PBE-D3, BLYP-D3, B3LYP-D3, ωB97XD, MP2-F12, M06-2X, and CCSD(T)-F12b levels. 4033

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A

6 shows that the PBE functional significantly overestimates the CCSD(T)-F12b BEPMs of ice-Ih (for 2.7 kcal/mol) and mixed

and PBE-D3 increase from 1.01 to 2.03 kcal/mol. However, if only the liquid structures 7−10 are taken into account, the corresponding MUEs for the four functionals range only from 0.18 to 0.21 kcal/mol. Thus, the large deviations of REPMs with different methods come from the ice and mixed ice−liquid structures. Table 2 and Figure 4 also show that the REPMs predicted by each method increase for structures 1−5 but remain similar for 7−10. For example, the corresponding changes of PBE-D3 REPMs are 3.94 and 0.16 kcal/mol, respectively. The similar REs of liquid clusters may arise from the statistically averaging of the intermolecular interaction energies. Thus, the GEBF-CCSD(T)-F12b method could benchmark different methods for those large ice−liquid clusters in the melting process. Furthermore, Table 2 and Figure S4 suggest that the empirical dispersion corrections only slightly improve the DFT REPMs. For example, the D3 correction reduces the MUE of B3LYP REPM for only 0.03 kcal/mol. Additionally, Figure S5 shows that FFs (SPC, TIP3P, TIP4P, and TIP5P) provide a similar changing trend of the REPMs as DFT and electron correlation levels but give larger REPM differences. Thus, we can conclude that the changing of noncovalent interactions significantly influences the performances of DFT, electron correlation, and FF methods for the REs in the melting process of ice-Ih. Binding Energies of Ice−Liquid 144-mers. Table S1 shows that the behaviors of BEPMs of the ten ice−liquid 144mers are similar as the corresponding REPMs. It indicates that the large deviations of BEPMs with different methods come from the ice and mixed ice−liquid structures, the BEPMs predicted by each method decrease for structures 1−5 but remain similar for 7−10, and the MP2-F12 method provides satisfactory BEPMs with MUE being only 0.30 kcal/mol. However, almost all functionals except B3LYP give MUEs larger than 1.5 kcal/mol, which indicates that the best functionals for the REs may be inconsistent with those for the BEs. Also, the empirical dispersion corrections only slightly improve the ωB97X BEPMs. The dispersion corrections are not recommended for PBE, BLYP, and B3LYP BEs because that those functionals without corrections already overestimate the BEs. Additionally, Figure 5 and Table S2 show that two basis sets, cc-pVDZ-F12 and cc-pVTZ-F12, give very similar MP2-F12 BEPMs with the MUE being only 0.09 kcal/mol. It indicates that the explicitly correlated F12 methods could yield nearly converged results with the cc-pVDZ-F12 basis set. Figure

Figure 6. Comparisons of the GEBF BEPMs of the ten ice−liquid 144-mers at PBE, BLYP, and MP2-F12 levels relative to GEBFCCSD(T)-F12b results.

ice−liquid structures for 2.7−1.6 kcal/mol, whereas only slightly overestimates those of liquid structures (7−10) for around 0.4 kcal/mol. It could explain why the PBE-based AIMD simulations overestimate the density of ice-Ih as 0.97 g/ mL (by comparing with experimental value, 0.93 g/mL) but underestimate that of liquid water as 0.88 g/mL.90−92 In addition, Figure 6 also shows that BLYP functional underestimates the CCSD(T)-F12b BEPMs of mixed ice−liquid structures and significantly underestimates those of liquid structures (7−10) for about 1.5 kcal/mol. It could illustrate why the BLYP-based AIMD simulations slightly underestimate the density of ice Ih as 0.91 g/mL but seriously underestimate that of liquid water as 0.75 g/L.92,93 Charges, Dipole Moments, and Forces of Ice−Liquid 144-mers. Table S3 shows that the average natural charges on O atoms do not change for each DFT or MP2 method in structures 7−10, whereas the absolute values decrease gradually in structures 1−6, probably due to the changes of hydrogenbonding interactions. The fluctuations of charges indicates that the classical MD simulations with fixed charges are difficult for describing the melting process of the ice crystal. By comparing with the MP2 natural charges, the M06-2X functional almost gives identical results, whereas the other functionals slightly underestimate the absolute values of charges. It is consistent with the performances of those functionals for REPMs. In addition, Table S3 also indicates that the ωB97XD and M06-2X functionals could reproduce the MP2 dipole moments, whereas the PBE functional overestimates the corresponding results, with the 6-311++G(d,p) basis set. Finally, Table S4 shows that GEBF could reproduce the conventional M06-2X forces with the root-mean-square deviations (RMSDs) being only about 0.0001 hartree/bohr with either the 6-31G(d) or 6-311G(d) basis sets. It indicates the feasibility of the GEBF approach for AIMD simulations. Periodic Ice−Liquid (H2O)64 Systems. To evaluate the performances of periodic water model for the melting process of ice crystal, we have investigated ten periodic (H2O)64 systems in a similar way. The hydrogen-bonding analyses in Table 3 also concludes that the decreasing (strong) HBs and average coordinated numbers in the ice-Ih crystal and mixed ice−liquid systems reduces the hydrogen-bonding interactions.

Figure 5. Comparisons of the GEBF BEPMs of the ten ice−liquid 144-mers at MP2-F12 level with the cc-pVDZ-F12 and cc-pVTZ-F12 basis sets. 4034

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A

Table 3. Comparisons of the Number of Hydrogen Bonds (n), Average Coordinated Numbers (⟨k⟩), and the Number of kCoordinated Water Molecules (k = 0, 1, ..., 5) of the Ten Periodic Ice−Liquid (H2O)64 Systemsa number of k-coordinated molecules (k) n

structure 1 2 3 4 5 6 7 8 9 10 a

128 124 127 114 118 110 103 100 99 97

(128) (108) (106) (97) (96) (86) (78) (72) (69) (74)

⟨k⟩ 4.00 3.88 3.97 3.56 3.69 3.44 3.22 3.12 3.09 3.03

(4.00) (3.38) (3.31) (3.03) (3.00) (2.69) (2.44) (2.25) (2.16) (2.31)

0 0 0 0 0 0 0 0 0 0 0

(0) (0) (0) (0) (0) (1) (1) (4) (3) (3)

1 0 0 0 1 1 1 0 2 1 4

2

(0) (1) (1) (3) (4) (8) (10) (12) (12) (10)

0 0 0 3 3 6 11 15 15 11

3

(0) (5) (7) (11) (14) (15) (22) (24) (26) (23)

0 8 5 19 15 25 30 22 27 28

4

(0) (27) (27) (31) (25) (27) (22) (13) (18) (20)

64 56 56 41 41 28 21 23 19 21

5

(64) (31) (29) (19) (20) (12) (9) (10) (5) (8)

0 0 3 0 4 4 2 2 2 0

(0) (0) (0) (0) (1) (1) (0) (1) (0) (0)

The corresponding values for strong hydrogen bonds included in parentheses.

shows similar behaviors of REPMs (or BEPMs) as those for 144-mers, including the following: the large deviations of REPMs (or BEPMs) with different methods usually come from the ice and mixed ice−liquid structures; the REPMs (or BEPMs) predicted by each method increase (or decrease) gradually for structures 1−5 but remain similar for 7−10; the MP2-F12 method provides satisfactory REPMs (or BEPMs) with MUE being only 0.51 kcal/mol (or 0.53 kcal/mol); the ωB97X and M06-2X functionals could provide reasonable REPMs. Thus, both the large ice−liquid 144-mers and periodic (H2O)64 systems could be approximately employed as models for bulk water. Furthermore, Figure 8 could also be employed to explain the densities deviations (from the experimental values) of ice-Ih and liquid water predicted by PBE- or BLYPbased AIMD simulations. Thus, in the melting process of ice, the PBC-GEBF calculations of periodic ice−liquid systems could qualitatively give conclusions similar to those of the GEBF calculations of large ice−liquid clusters.

Here, the average coordinated number of ice-Ih crystal is 4.0 due to its periodic structure. The PBC-GEBF-X [X = CCSD(T), MP2-F12, DFT] calculations are carried out to investigate these periodic systems. The empirical dispersion corrections are not employed because they do not apparently improve the results of 144-mers. Figure 7 and Tables 4 and S5

4. CONCLUSIONS With the GEBF and PBC-GEBF approaches, we have investigated various DFT and electron correlation methods for the relative energies and binding energies in the melting process of ice-Ih crystal. Both the large cluster model and periodic model are employed to approximately describe bulk water. The studied systems, including ten ice−liquid 144-mers

Figure 7. Comparisons of the PBC-GEBF REPMs of the ten periodic ice−liquid (H2O)64 systems at PBE, BLYP, B3LYP, ωB97X, M06-2X, MP2-F12, and CCSD(T)-F12b levels.

Table 4. Comparisons of the PBC-GEBF REPMs of the Ten Periodic Ice−Liquid (H2O)64 Systems at Various DFT Levels with the 6-311++G(d,p) Basis Set, and at MP2-F12 and CCSD(T)-F12b Levels with the cc-pVDZ-F12 Basis Seta structure

PBE

BLYP

B3LYP

ωB97X

M06-2X

MP2-F12

CCSD(T)-F12b

1 2 3 4 5 6 7 8 9 10 NPEb MUEc

−4.68 −2.65 −1.97 −1.23 −0.57 0.00 0.63 0.53 0.66 0.79 1.78 1.16

−4.77 −2.66 −2.09 −1.26 −0.54 0.00 0.55 0.48 0.52 0.66 1.73 1.25

−4.68 −2.45 −1.86 −1.12 −0.48 0.00 0.66 0.52 0.60 0.76 1.75 1.15

−3.83 −1.77 −1.32 −0.88 −0.38 0.00 0.62 0.47 0.56 0.72 0.86 0.54

−3.30 −1.35 −0.99 −0.74 −0.28 0.00 0.50 0.39 0.42 0.62 0.94 0.82

−3.53 −1.83 −1.37 −0.90 −0.43 0.00 0.61 0.54 0.55 0.69 0.52 0.51

−3.52 −2.17 −1.39 −0.94 −0.87 0.00 0.59 0.51 0.49 0.17 0.00 0.00

a

All energies are in kcal/mol. bNonparallelity error relative to the CCSD(T)-F12b results. cMaximum unsigned error relative to the CCSD(T)-F12b results. 4035

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A ORCID

Shuhua Li: 0000-0001-6756-057X Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (Grant Nos. 21473087, 21333004, and 21673110) and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase). We are grateful to the IBM Blade cluster system provided by the High Performance Computing Center at Nanjing University for some calculations.



Figure 8. Comparisons of the PBC-GEBF BEPMs of the ten periodic ice−liquid (H2O)64 systems at PBE, BLYP, and MP2-F12 levels relative to PBC-GEBF-CCSD(T)-F12b results.

and ten periodic ice−liquid (H2O)64 systems, are taken from MD simulations of the melting process of ice-Ih crystals. Our results show that both the large cluster model and periodic model could provide consistent results to approximately describe the bulk ice−liquid water. The REPMs and BEPMs predicted by various methods are extremely different and different for the ice-Ih structure (1) and mixed ice−liquid structures (2−6), respectively, but similar for pure liquid structures (7−10). The REPMs (or BEPMs) predicted by each method increase (or decrease) gradually for structures 1−6 but remain similar for 7−10. The hydrogen-bonding and RDG analyses suggest that the changing of noncovalent interactions significantly influences the performances of DFT, electron correlation, and FF methods for the REs and BEs in the melting process of ice-Ih. The benchmark CCSD(T)-F12b results indicate that the MP2-F12 method provides satisfactory REs and BEs, and the M06-2X, ωB97XD, and ωB97X functionals could also provide reasonable REs. The PBE and BLYP BEPMs (relative to CCSD(T)-F12b ones) could be employed to explain why the PBE-based AIMD simulations overestimate and underestimate the densities of ice and liquid water, respectively, and the BLYP-based simulations slightly and severely underestimate the densities of ice and liquid water, respectively. The PBC-GEBF-MP2-F12-based AIMD simulations are expected to be applied to the phase transitions of ice crystal.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.7b03376. The GEBF BEPMs of the ten 144-mers and ten periodic (H2O)64 systems at different levels, the GEBF-MP2-F12 BEPMs with the cc-pVXZ-F12 (X = D, T) basis sets, the charges on O atoms, the RMSD forces, graph of the total energy in the functions of time in the simulations of iceIh (H2O)64, the RDG figures, and bar charts of the REPMs at GEBF-DFT (with and without empirical dispersion corrections) and FFs of the ten 144-mers (PDF)



REFERENCES

(1) Meng, S.; Greenlee, L. F.; Shen, Y. R.; Wang, E. Basic science of water: Challenges and current status towards a molecular picture. Nano Res. 2015, 8, 3085−3110. (2) Zhang, X.; Sun, P.; Yan, T.; Huang, Y.; Ma, Z.; Zou, B.; Zheng, W.; Zhou, J.; Gong, Y.; Sun, C. Q. Water’s phase diagram: From the notion of thermodynamics to hydrogen-bond cooperativity. Prog. Solid State Chem. 2015, 43, 71−81. (3) Huang, Y.; Zhang, X.; Ma, Z.; Zhou, Y.; Zheng, W.; Zhou, J.; Sun, C. Q. Hydrogen-bond relaxation dynamics: Resolving mysteries of water ice. Coord. Chem. Rev. 2015, 285, 109−165. (4) Sellberg, J. A.; Huang, C.; McQueen, T. A.; Loh, N. D.; Laksmono, H.; Schlesinger, D.; Sierra, R. G.; Nordlund, D.; Hampton, C. Y.; Starodub, D.; et al. Ultrafast X-ray probing of water structure below the homogeneous ice nucleation temperature. Nature 2014, 510, 381−384. (5) Drechsel-Grau, C.; Marx, D. Tunnelling in chiral water clusters: Protons in concert. Nat. Phys. 2015, 11, 216−218. (6) Wren, S. N.; Donaldson, D. J. Glancing-angle Raman spectroscopic probe for reaction kinetics at water surfaces. Phys. Chem. Chem. Phys. 2010, 12, 2648−2654. (7) Sun, C. Q.; Sun, Y. Springer Ser. Chem. Phys. 2016, 113, 1−512. (8) Wang, J.; Yoo, S.; Bai, J.; Morris, J. R.; Zeng, X. C. Melting temperature of ice Ih calculated from coexisting solid-liquid phases. J. Chem. Phys. 2005, 123, 036101. (9) Vega, C.; Sanz, E.; Abascal, J. L. F. The melting temperature of the most common models of water. J. Chem. Phys. 2005, 122, 114507. (10) Gaiduk, A. P.; Gygi, F.; Galli, G. Density and Compressibility of Liquid Water and Ice from First-Principles Simulations with Hybrid Functionals. J. Phys. Chem. Lett. 2015, 6, 2902−2908. (11) Yoo, S.; Zeng, X. C.; Xantheas, S. S. On the phase diagram of water with density functional theory potentials: The melting temperature of ice Ih with the Perdew-Burke-Ernzerhof and BeckeLee-Yang-Parr functionals. J. Chem. Phys. 2009, 130, 221102. (12) Yoo, S.; Xantheas, S. S. Communication: The effect of dispersion corrections on the melting temperature of liquid water. J. Chem. Phys. 2011, 134, 121105. (13) Ben, M. D.; Schonherr, M.; Hutter, J.; VandeVondele, J. Bulk Liquid Water at Ambient Temperature and Pressure from MP2 Theory. J. Phys. Chem. Lett. 2013, 4, 3753−3759. (14) Brorsen, K. R.; Willow, S. Y.; Xantheas, S. S.; Gordon, M. S. The Melting Temperature of Liquid Water with the Effective Fragment Potential. J. Phys. Chem. Lett. 2015, 6, 3555−3559. (15) Willow, S. Y.; Salim, M. A.; Kim, K. S.; Hirata, S. Ab initio molecular dynamics of liquid water using embedded-fragment secondorder many-body perturbation theory towards its accurate property prediction. Sci. Rep. 2015, 5, 14358. (16) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. J. Chem. Phys. 2004, 120, 300− 311. (17) Dahlke, E. E.; Olson, R. M.; Leverentz, H. R.; Truhlar, D. G. Assessment of the Accuracy of Density Functionals for Prediction of

AUTHOR INFORMATION

Corresponding Author

*W. Li. E-mail: [email protected]. 4036

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A Relative Energies and Geometries of Low-Lying Isomers of Water Hexamers. J. Phys. Chem. A 2008, 112, 3976−3984. (18) Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. Assessing the Accuracy of Density Functional and Semiempirical Wave Function Methods for Water Nanoparticles: Comparing Binding and Relative Energies of (H2O)16 and (H2O)17 to CCSD(T) Results. J. Chem. Theory Comput. 2013, 9, 995−1006. (19) Wang, K.; Li, W.; Li, S. Generalized Energy Based Fragmentation CCSD(T)-F12 Method and Application to Water Clusters (H2O)20. J. Chem. Theory Comput. 2014, 10, 1546−1553. (20) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479−483. (21) Bates, D. M.; Tschumper, G. S. CCSD(T) Complete Basis Set Limit Relative Energies for Low-Lying Water Hexamer Structures. J. Phys. Chem. A 2009, 113, 3555−3559. (22) Bryantsev, V. S.; Diallo, M. S.; van Duin, A. C. T.; Goddard, W. A. Evaluation of B3LYP, X3LYP, and M06-Class Density Functionals for Predicting the Binding Energies of Neutral, Protonated, and Deprotonated Water Clusters. J. Chem. Theory Comput. 2009, 5, 1016−1026. (23) Yoo, S.; Aprà, E.; Zeng, X. C.; Xantheas, S. S. High-Level Ab Initio Electronic Structure Calculations of Water Clusters (H2O)16 and (H2O)17: A New Global Minimum for (H2O)16. J. Phys. Chem. Lett. 2010, 1, 3122−3127. (24) Temelso, B.; Archer, K. A.; Shields, G. C. Benchmark Structures and Binding Energies of Small Water Clusters with Anharmonicity Corrections. J. Phys. Chem. A 2011, 115, 12034−12046. (25) Temelso, B.; Renner, C. R.; Shields, G. C. Importance and Reliability of Small Basis Set CCSD(T) Corrections to MP2 Binding and Relative Energies of Water Clusters. J. Chem. Theory Comput. 2015, 11, 1439−1448. (26) Miliordos, E.; Aprà, E.; Xantheas, S. S. Optimal geometries and harmonic vibrational frequencies of the global minima of water clusters (H2O)n, n = 2−6, and several hexamer local minima at the CCSD(T) level of theory. J. Chem. Phys. 2013, 139, 114302. (27) Miliordos, E.; Xantheas, S. S. An accurate and efficient computational protocol for obtaining the complete basis set limits of the binding energies of water clusters at the MP2 and CCSD(T) levels of theory: Application to (H2O)m, m = 2−6, 8, 11, 16, and 17. J. Chem. Phys. 2015, 142, 234303. (28) Li, S.; Li, W.; Fang, T. An Efficient Fragment-Based Approach for Predicting the Ground-State Energies and Structures of Large Molecules. J. Am. Chem. Soc. 2005, 127, 7215−7226. (29) Li, W.; Li, S.; Jiang, Y. Generalized Energy-Based Fragmentation Approach for Computing the Ground-State Energies and Properties of Large Molecules. J. Phys. Chem. A 2007, 111, 2193−2199. (30) Li, S.; Li, W.; Ma, J. Generalized Energy-Based Fragmentation Approach and Its Applications to Macromolecules and Molecular Aggregates. Acc. Chem. Res. 2014, 47, 2712−2720. (31) Fang, T.; Li, W.; Gu, F.; Li, S. Accurate Prediction of Lattice Energies and Structures of Molecular Crystals with Molecular Quantum Chemistry Methods. J. Chem. Theory Comput. 2015, 11, 91−98. (32) Fang, T.; Jia, J.; Li, S. Vibrational Spectra of Molecular Crystals with the Generalized Energy-Based Fragmentation Approach. J. Phys. Chem. A 2016, 120, 2700−2711. (33) Li, W.; Li, Y.; Lin, R.; Li, S. Generalized Energy-Based Fragmentation Approach for Localized Excited States of Large Systems. J. Phys. Chem. A 2016, 120, 9667−9677. (34) Deev, V.; Collins, M. A. Approximate ab initio energies by systematic molecular fragmentation. J. Chem. Phys. 2005, 122, 154102. (35) Collins, M. A.; Deev, V. A. Accuracy and efficiency of electronic energies from systematic molecular fragmentation. J. Chem. Phys. 2006, 125, 104104. (36) Addicoat, M. A.; Collins, M. A. Accurate treatment of nonbonded interactions within systematic molecular fragmentation. J. Chem. Phys. 2009, 131, 104103.

(37) Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Fast electron correlation methods for molecular clusters in the ground and excited states. Mol. Phys. 2005, 103, 2255−2265. (38) Huang, L.; Massa, L.; Karle, J. Kernel energy method illustrated with peptides. Int. J. Quantum Chem. 2005, 103, 808−817. (39) He, X.; Zhang, J. Z. H. The generalized molecular fractionation with conjugate caps/molecular mechanics method for direct calculation of protein energy. J. Chem. Phys. 2006, 124, 184703. (40) Wang, X.; Liu, J.; Zhang, J. Z. H.; He, X. Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method for Full Quantum Mechanical Calculation of Protein Energy. J. Phys. Chem. A 2013, 117, 7149−7161. (41) Bettens, R. P. A.; Lee, A. M. A New Algorithm for Molecular Fragmentation in Quantum Chemical Calculations. J. Phys. Chem. A 2006, 110, 8777−8785. (42) Le, H.-A.; Tan, H.-J.; Ouyang, J. F.; Bettens, R. P. A. Combined Fragmentation Method: A Simple Method for Fragmentation of Large Molecules. J. Chem. Theory Comput. 2012, 8, 469−478. (43) Gadre, S. R.; Ganesh, V. Molecular Tailoring approach: towards PC-based ab initio treatment of large molecules. J. Theor. Comput. Chem. 2006, 05, 835−855. (44) Kavathekar, R.; Khire, S.; Ganesh, V.; Rahalkar, A. P.; Gadre, S. R. WebMTA: A web-interface for ab initio geometry optimization of large molecules using molecular tailoring approach. J. Comput. Chem. 2009, 30, 1167−1173. (45) Isegawa, M.; Wang, B.; Truhlar, D. G. Electrostatically Embedded Molecular Tailoring Approach and Validation for Peptides. J. Chem. Theory Comput. 2013, 9, 1381−1393. (46) Sahu, N.; Singh, G.; Nandi, A.; Gadre, S. R. Toward an Accurate and Inexpensive Estimation of CCSD(T)/CBS Binding Energies of Large Water Clusters. J. Phys. Chem. A 2016, 120, 5706−5714. (47) Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded ManyBody Expansion for Large Systems, with Applications to Water Clusters. J. Chem. Theory Comput. 2007, 3, 46−53. (48) Sorkin, A.; Dahlke, E. E.; Truhlar, D. G. Application of the Electrostatically Embedded Many-Body Expansion to Microsolvation of Ammonia in Water Clusters. J. Chem. Theory Comput. 2008, 4, 683−688. (49) Mayhall, N. J.; Raghavachari, K. Molecules-in-Molecules: An Extrapolated Fragment-Based Approach for Accurate Calculations on Large Molecules and Materials. J. Chem. Theory Comput. 2011, 7, 1336−1343. (50) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646−9663. (51) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632−672. (52) Richard, R. M.; Herbert, J. M. A generalized many-body expansion and a unified view of fragment-based methods in electronic structure theory. J. Chem. Phys. 2012, 137, 064113. (53) Lao, K. U.; Liu, K.-Y.; Richard, R. M.; Herbert, J. M. Understanding the many-body expansion for large systems. II Accuracy considerations. J. Chem. Phys. 2016, 144, 164105. (54) Collins, M. A.; Bettens, R. P. A. Energy-Based Molecular Fragmentation Methods. Chem. Rev. 2015, 115, 5607−5642. (55) Raghavachari, K.; Saha, A. Accurate Composite and FragmentBased Quantum Chemical Models for Large Molecules. Chem. Rev. 2015, 115, 5643−5677. (56) Ten-no, S. Initiation of explicitly correlated Slater-type geminal theory. Chem. Phys. Lett. 2004, 398, 56−61. (57) Werner, H.-J.; Adler, T. B.; Manby, F. R. General orbital invariant MP2-F12 theory. J. Chem. Phys. 2007, 126, 164102. (58) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (59) Li, W. Linear scaling explicitly correlated MP2-F12 and ONIOM methods for the long-range interactions of the nanoscale clusters in methanol aqueous solutions. J. Chem. Phys. 2013, 138, 014106. 4037

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038

Article

The Journal of Physical Chemistry A (60) Yuan, D.; Shen, X.; Li, W.; Li, S. Are fragment-based quantum chemistry methods applicable to medium-sized water clusters? Phys. Chem. Chem. Phys. 2016, 18, 16491−16500. (61) Tao, Y.; Zou, W.; Jia, J.; Li, W.; Cremer, D. Different Ways of Hydrogen Bonding in Water - Why Does Warm Water Freeze Faster than Cold Water? J. Chem. Theory Comput. 2017, 13, 55−76. (62) Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys. 1921, 369, 253−287. (63) Foster, J. P.; Weinhold, F. Natural hybrid orbitals. J. Am. Chem. Soc. 1980, 102, 7211−7218. (64) Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural population analysis. J. Chem. Phys. 1985, 83, 735−746. (65) Li, S.; Li, W.; Fang, T.; Ma, J.; Hua, W.; Hua, S.; Jiang, Y. LSQC Program, Version 2.2; Nanjing University: Nanjing, 2012; http://itcc. nju.edu.cn/lsqc. (66) Li, W.; Chen, C.; Zhao, D.; Li, S. LSQC: Low scaling quantum chemistry program. Int. J. Quantum Chem. 2015, 115, 641−646. (67) Hua, W.; Fang, T.; Li, W.; Yu, J.-G.; Li, S. Geometry Optimizations and Vibrational Spectra of Large Molecules from a Generalized Energy-Based Fragmentation Approach. J. Phys. Chem. A 2008, 112, 10864−10872. (68) Hua, S.; Hua, W.; Li, S. An Efficient Implementation of the Generalized Energy-Based Fragmentation Approach for General Large Molecules. J. Phys. Chem. A 2010, 114, 8126−8134. (69) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (70) Becke, A. D. Density-functional thermochemistry. III The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (71) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (72) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (73) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (74) Chai, J.-D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (75) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (76) Zhao, Y.; Truhlar, D. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (77) Peterson, K. A.; Adler, T. B.; Werner, H.-J. Systematically Convergent Basis Sets for Explicitly Correlated Wavefunctions: the Atoms H, He, B-ne, and Al-ar. J. Chem. Phys. 2008, 128, 084102. (78) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498−6506. (79) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792− 794. (80) Ponder, J. W.; Richards, F. M. An efficient newton-like method for molecular mechanics energy minimization of large molecules. J. Comput. Chem. 1987, 8, 1016−1024. (81) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision B.01; Gaussian Inc.: Wallingford, CT, 2009. (82) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2012, 2, 242−253.

(83) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; et al. MOLPRO, version 2012.1, a package of ab initio programs 2012; http://www.molpro.net. (84) Lu, T.; Chen, F. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 2012, 33, 580−592. (85) Zhao, Y.; Truhlar, D. G. Design of Density Functionals That Are Broadly Accurate for Thermochemistry, Thermochemical Kinetics, and Nonbonded Interactions. J. Phys. Chem. A 2005, 109, 5656−5667. (86) Zhao, Y.; Truhlar, D. G. Hybrid Meta Density Functional Theory Methods for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions: The MPW1B95 and MPWB1K Models and Comparative Assessments for Hydrogen Bonding and van der Waals Interactions. J. Phys. Chem. A 2004, 108, 6908−6918. (87) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364−382. (88) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (89) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (90) Santra, B.; Klime, J.; Tkatchenko, A.; Alfè, D.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. On the accuracy of van der Waals inclusive density-functional theory exchange-correlation functionals for ice at ambient and high pressures. J. Chem. Phys. 2013, 139, 154702. (91) Whalley, E. Energies of the phases of ice at zero temperature and pressure. J. Chem. Phys. 1984, 81, 4087−4092. (92) Schmidt, J.; VandeVondele, J.; Kuo, I.-F. W.; Sebastiani, D.; Siepmann, J. I.; Hutter, J.; Mundy, C. J. Isobaric-Isothermal Molecular Dynamics Simulations Utilizing Density Functional Theory: An Assessment of the Structure and Density of Water at Near-Ambient Conditions. J. Phys. Chem. B 2009, 113, 11959−11964. (93) Feibelman, P. J. Lattice match in density functional calculations: ice Ih vs. β-AgI. Phys. Chem. Chem. Phys. 2008, 10, 4688−4691.

4038

DOI: 10.1021/acs.jpca.7b03376 J. Phys. Chem. A 2017, 121, 4030−4038