Fragment Quantum Mechanical Method for Large-Sized Ion–Water

Apr 5, 2017 - Department of Chemistry, New York University, New York, New York 10003, United ... Jinfeng Liu , Xiao He , John Z. H. Zhang , Lian-Wen Q...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JCTC

Fragment Quantum Mechanical Method for Large-Sized Ion−Water Clusters Jinfeng Liu,† Lian-Wen Qi,*,† John Z. H. Zhang,‡,§,∥ and Xiao He*,‡,§ †

State Key Laboratory of Natural Medicines, Department of Basic Medicine and Clinical Pharmacy, China Pharmaceutical University, Nanjing, 210009, China ‡ School of Chemistry and Molecular Engineering, East China Normal University, Shanghai, 200062, China § NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai, 200062, China ∥ Department of Chemistry, New York University, New York, New York 10003, United States S Supporting Information *

ABSTRACT: Fragmentation methods have been widely studied for computing quantum mechanical (QM) energy of medium-sized water clusters, but less attention has been paid to large-sized ion−water clusters, in which many-body QM interaction is more significant, because of the charge-transfer effect between ions and water molecules. In this study, we utilized electrostatically embedded generalized molecular fractionation (EEGMF) method for full QM calculation of the large-sized ion−water clusters (up to 15 Na+ and 15 Cl− ions solvated with 119 water molecules). Through systematic validation using different fragment sizes, we show that, by using distance thresholds of 6 Å for both the two-body and three-body QM interactions, the EE-GMF method is capable of providing accurate ground-state energies of large-sized ion−water clusters at different ab initio levels (including HF, B3LYP, M06-2X, and MP2) with significantly reduced computational cost. The deviations of EE-GMF from full system calculations are within a few kcal/mol. The result clearly shows that the calculated energies of the ion−water clusters using EEGMF are close to converge after the distance thresholds are larger than 6 Å for both the two-body and three-body QM interactions. This study underscores the importance of the three-body interactions in ion−water clusters. The EE-GMF method can also accurately reproduce the relative energy profiles of the ion−water clusters. caps (MFCC) approach,3,7−12 the systematic fragmentation method (SFM),13−15 the adjustable density matrix assembler (ADMA) approach,16−18 the molecular tailoring approach (MTA),19−21 the generalized energy-based fragmentation (GEBF) method,22,23 the electrostatically embedded manybody (EE-MB) expansion approach,24−26 the explicit polarization (X-Pol) potential,27,28 the combined fragmentation method (CFM),29 the molecules-in-molecules (MIM) approach,30 the many-overlapping-body (MOB) expansion approach,31 the generalized many-body expansion (GMBE) approach,32,33 and the automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) method.34−36 A comprehensive review of the fragmentation QM methods can be found in two recent review articles, written by Gordon and co-workers2 and Collins and co-workers,1 respectively. The electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method is developed based on the framework of the MFCC approach for accurate calculation of protein energy.3 Numerical tests for

1. INTRODUCTION At present, it is still not practical to apply high-level wave function theories to macromolecules, which contain hundreds of atoms or more. In order to extend ab initio QM calculations to those large systems, various linear-scaling or fragmentation methods have been developed. Because of the computational efficiency and easy implementation, the development of fragmentation methods has attracted significant interest over the past decade.1,2 The fragmentation approach is based on the chemical locality of large molecules, which assumes that the local chemical environment of each subsystem is less influenced by the regions that are far away from the central fragment of interest. Based on this chemical intuition, a large system can be divided into several subsystems (fragments), according to the specific fragmentation scheme, and energies of the fragments can be separately computed at the QM level. The total energy or other properties of the entire system is then approximately derived by taking proper combinations of the corresponding terms of individual subsystems.3 A range of fragmentation QM methods for large systems have been proposed, including the fragment molecular orbital (FMO) method,4−6 the molecular fractionation with conjugate © XXXX American Chemical Society

Received: February 13, 2017 Published: April 5, 2017 A

DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

the conventional full system energies. Finally, conclusions are presented in section 4.

various benchmark systems have demonstrated the efficacy of the EE-GMFCC method at different levels of theory, as compared to the full system calculations for proteins. Furthermore, the EE-GMFCC approach has been successfully utilized for protein structure optimization,37 protein vibrational spectrum calculation,37 protein−ligand binding affinity prediction,38 and ab initio molecular dynamics simulation on proteins in explicit solvent.39 The EE-GMFCC method is an efficient and accurate fragmentation approach for studying the properties of large biomolecules at the QM level. Water clusters are often selected as benchmark systems for testing the accuracy of fragmentation methods by comparing with the conventional full system QM calculations. Truhlar and co-workers have extensively applied the EE-MB method on QM calculations of small- to medium-sized water clusters or mixed clusters containing both ions and water molecules.24−26 However, the size of the water cluster used in their studies is relatively small, which normally contains fewer than 30 water molecules. The fragmentation methods, which show good performance in these relatively small water clusters, are usually difficult to guarantee the same performance for large-sized water clusters which contain hundreds of water molecules because the errors may accumulate in the fragmentation methods.40−42 To the best of our knowledge, a systematic assessment of the performance of fragmentation method on large-sized water clusters (with over 100 water molecules) has not been reported previously. Furthermore, it is also worth investigating the accuracy of the fragmentation methods on large-sized ion−water clusters with high ion concentrations. The ion−water interactions are regarded as the fundamental components involved in various chemical, biological, and environmental processes.43 Since most of the relevant chemical reactions take place in solution, reliable methods for accurate theoretical treatment on chemical species in solution have always been in high demand.44 The crucial role played by the ion−water complex in biology underlines the importance of an accurate and detailed knowledge of the structure and dynamics of these species in solution.45 Because of the long-range electrostatic polarization effect arising from charged ions and charge transfer between ions and water molecules, the manybody QM effect in the ion−water cluster is more significant than that of the pure water cluster, which presents a grand challenge for fragmentation methods.46,47 Therefore, the development of a robust and broadly applicable method, that can accurately calculate the energy of large-sized ion−water clusters with a reduced cost, is of great importance for highlevel ab initio studies on aqueous solutions. In the present work, based on the EE-GMFCC scheme, we develop an electrostatically embedded generalized molecular fractionation (EE-GMF) method to treat molecular clusters and apply it for the total energy calculations of the large-sized water and ion−water clusters. The performance of the EEGMF method is systematically evaluated, including the dependence of EE-GMF energy with the fragment size and the use of background charges. The convergence of the EEGMF total energy versus the distance threshold for the local two-body and three-body QM interactions is also investigated exclusively. This paper is organized as follows. Section 2 provides the detailed description of the EE-GMF scheme for total energy calculations of clusters and the computational details. In section 3, we systematically evaluated the performance of the EE-GMF method on large-sized water and ion− water clusters at various QM levels and compared them with

2. THEORETICAL APPROACHES The EE-GMFCC approach was proposed for calculating the total energy of proteins. In the framework of the EE-GMFCC method, a protein is decomposed to a series of properly capped amino acid fragments by cutting through the peptide bond, and then the subsystem-based energies of neighboring units and two-body interaction energies between non-neighboring units, which are not sequentially connected but spatially in close contact, are taken together to approximately reproduce the total energy of the protein. The QM calculation of each fragment is carried out with the embedding field, which represents the electrostatic environment of the rest of the molecule.3 The detailed description of the EE-GMFCC method can be found in ref 3. For water and ion−water clusters, each water molecule or ion (grouped with its nearby water molecules or not) naturally becomes a single fragment, and it does not need to cut chemical bond or add molecular caps. In this work, we developed the EEGMF method, which removes the conjugate cap (CC) treatment from EE-GMFCC, for specifically dealing with molecular clusters. The energy calculation in the EE-GMF scheme is similar to that in the EE-GMFCC approach, except for the contributions of conjugate caps (overlapping fragments), because there is no molecular cap for each fragment in the molecular cluster. The energy of each fragment and the interaction energy between two fragments that are spatially in close contact are computed by QM, while the two-body interaction energies between distant fragment pairs are described by classical Coulomb interactions. All QM calculations are embedded in the electrostatic field of the point charges representing the remaining system to account for the environmental effect. The EE-GMFCC scheme for protein is truncated at the two-body interactions, while the three- and higher-body interactions are not explicitly considered, which are implicitly incorporated in the electrical embedding scheme. However, for accurate description of the intermolecular interaction in molecular clusters, three-body interactions sometimes have non-negligible contributions, especially for the ion−water cluster. The total EE-GMF energy of the water and ion−water clusters can be expressed using the following formula (the derivation of this term is provided in the Appendix): N EE ‐ GMF Emolecular cluster

N−1

= ∑ Eĩ +

∑ ∑

ΔEij̃

i=1 j=i+1 |R ij|≤ λ1

i=1 N−2

+

N

N−1

∑ ∑

N



̃ ΔEijk

i=1 j=i+1 k=j+1 |R ij|≤ λ 2 |R ik|≤ λ 2 |R jk|≤ λ 2 N−1



N

∑ ∑ ∑∑ i=1 j=i+1 m∈i n∈j |R ij|> λ1

qm(i)qn(j) R m(i)n(j) (1)

where N is the number of the fragments. The local two-body QM interaction is defined by ΔEij̃ = Eij̃ − Eĩ − Ej̃ B

(2) DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Constructed (a) water cluster (140 water molecules), (b) NaCl−water cluster from 3 M NaCl solution (136 water molecules, 8 Na+ cations, and 7 Cl− anions), and (c) NaCl−water cluster from 6 M NaCl solution (119 water molecules, 15 Na+ cations, and 15 Cl− anions) for HF and DFT calculations with 6-31G* basis set, and the reduced (d) water cluster (79 water molecules), (e) NaCl−water cluster from 3 M NaCl solution (70 water molecules, 5 Na+ cations, and 5 Cl− anions), (f) NaCl−water cluster from 6 M NaCl solution (61 water molecules, 10 Na+ cations, and 10 Cl− anions) for MP2 calculations with the 6-31G* basis set and HF, DFT calculations with large basis sets, and the reduced (g) water cluster (50 water molecules), (h) NaCl−water cluster from 3 M NaCl solution (46 water molecules, 3 Na+ cations, and 3 Cl− anions), (i) NaCl− water cluster from 6 M NaCl solution (43 water molecules, 6 Na+ cations, and 6 Cl− anions) for MP2 calculations with large basis sets. The molecular clusters were extracted from corresponding equilibrated solutions, and the Cartesian coordinates are provided in the Supporting Information.

charge−charge Coulomb interactions. qn(j) denotes the point charge of the nth atom in the jth subsystem. The EE-GMF approach for dealing with the water cluster is similar to the EEMB method.24,25,48 The difference from EE-MB is that, in EEGMF, the two-body and three-body QM interactions are truncated by using distance thresholds. In addition, the distant two-body interaction is treated by classical Coulomb interactions to achieve linear scale for high-level ab initio calculation. In the framework of the EE-GMF scheme, the short-range interaction is treated by the QM method, while the long-range interaction is approximated using the classical Coulomb interactions. In this study, the convergence of the total EE-GMF energy, as a function of the distance threshold, is investigated for molecular clusters. We also tested different

in which the minimal distance Rij between two non-hydrogen atoms of fragments i and j is no larger than a predefined distance threshold λ1, and the local three-body QM interaction energy is defined by ̃ = Eijk ̃ − ΔEij̃ − ΔEik̃ − ΔEjk̃ − Eĩ − Ej̃ − Ek̃ ΔEijk

(3)

where the distance of heavy atoms among any two fragments from i, j, and k is less than or equal to the distance threshold λ2. Ẽ denotes the sum of the self-energy of the fragment (Ẽ i, Ẽ ij, and Ẽ ijk for monomer, dimer, and trimer, respectively), along with the interaction energy between the fragment and background charges of the remaining system. The doubly counted two-body interaction energies between distant fragments is deducted in the last term of eq 1, using pairwise C

DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. Total EE-GMF energy for the large-sized water cluster (140 water molecules) calculated at the levels of HF/6-31G* (panels (a) and (b)) and HF/6-311+G* (panels (c) and (d)), respectively, as a function of the distance threshold λ for two-body interactions (panels (a) and (c)) and three-body interactions (panels (b) and (d)). In panels (a) and (c), three-body interactions were not included. The distance threshold λ was varied only for two-body interactions. In panels (b) and (d), the distance threshold for two-body interactions was fixed at 7.0 Å. The distance threshold λ was varied only for three-body interactions. “Fragmentation (1)” takes each water molecule as a fragment, while “Fragmentation (2)” takes two nearest water molecules as a fragment.

charge model is more computationally efficient for large system calculation, and analytical energy gradients can be readily derived. Moreover, in our previous studies,3,51 we have demonstrated that utilizing the Amber charge model is a simpler way to include the electrostatic polarization with acceptable accuracy. Therefore, in this study, we employed the Amber charge model to approximately describe the electrostatic potential within the molecule. Full system and EE-GMF calculations were performed at HF, B3LYP, M06-2X, and MP2 levels, respectively, with various basis sets, including 6-31G*, 6311+G*, 6-311++G**, cc-pVDZ, and aug-cc-pVDZ. All QM calculations were carried out using the Gaussian09 package.52

fragmentation schemes in the EE-GMF approach for these clusters. To test the performance of the EE-GMF method for the water and ion−water clusters, we constructed a large-sized water cluster (140 water molecules), a NaCl-water cluster chosen from 3 M NaCl solution (136 water molecules, 8 Na+ cations, and 7 Cl− anions) and a NaCl-water cluster selected from 6 M NaCl solution (119 water molecules, 15 Na+ cations, and 15 Cl− anions), as shown in Figure 1. Cartesian coordinates of these three clusters are given in the Supporting Information. We first investigated the dependence of the EE-GMF total energy on the fragment size and the energy convergence versus the distance threshold for the local two-body and three-body QM interactions. The deviations of EE-GMF from the conventional full system energies are analyzed to assess the accuracy of this method. The contributions of the two-body and three-body interactions in the water and ion−water clusters are analyzed. Twenty configurations for each system were selected from 500 ps classical molecular dynamics (MD) simulation (at 300 K), using the SPC water model.49 The relative energies of these conformations are evaluated. In this study, the fixed charge models of the SPC water model and the Amber ff14SB force field50 is utilized to describe the embedding field. In this study, we did not use the self-consistently determined point changes derived from the same level of ab initio theory as the embedding field. The use of the fixed Amber

3. RESULTS AND DISCUSSION 3.1. Accuracy of the EE-GMF Approach for LargeSized Water Clusters. The objective of the EE-GMF method is to carry out the QM calculation for the large system as efficiently as possible while achieving a desired accuracy with reference to the full system result. Hence, the distance threshold for the two-body and three-body QM interactions must be rigorously validated to strike a compromise between the computational efficiency and attained accuracy. We first apply the EE-GMF approach to calculate the total energy of the large-sized water cluster at the HF/6-31G* level. The fragment size changes from one water molecule to two water molecules as a fragment to test the accuracy of the fragmentation scheme. D

DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

the electrostatic embedding scheme. Therefore, we conclude that including the environmental polarization effect in fragmentation approach significantly improves the convergence and accuracy of the calculated total energy. For the compromise between the computation time and calculation accuracy, we set the distance thresholds of the twobody and three-body QM interactions at λ1 = λ2 = 6.0 Å in the following calculations on the water cluster. As shown in Table 1, the deviation of the EE-GMF energy from the full system energy at the HF/6-31G* level is 6.31 kcal/mol when one water molecule is taken as a fragment, and this deviation can be decreased to 3.14 kcal/mol when two nearest water molecules are grouped as a fragment. We further applied the EE-GMF approach at DFT (B3LYP and M06-2X) and MP2 levels using the 6-31G* basis set. Table 2 shows that the energies of the water cluster calculated by EEGMF with two different fragment sizes are also in excellent agreement with the conventional full system B3LYP/6-31G* result, where the larger deviation in the case of one water molecule per fragment is just 3.39 kcal/mol. Table 3 shows the results obtained from M06-2X/6-31G* calculations. The deviation of 10.37 kcal/mol for the EE-GMF energy calculated with one water molecule per fragment is slightly larger. However, the deviation of EE-GMF energy will be largely decreased to 5.65 kcal/mol when two nearest water molecules are assigned as a fragment. Because the full system MP2 calculation on such a large water cluster is very time-consuming, we trimmed the water cluster to 79 water molecules to make the conventional full system MP2 calculation more affordable (see Figure 1d). Table 4 shows the full system MP2 result of that water cluster, together with the deviations of EE-GMF energies. For the fragmentation using two water molecules as a fragment, the error of EE-GMF energy is as small as 2.89 kcal/mol, and it becomes slightly larger (6.90 kcal/mol) when one water molecule is assigned as a fragment. We also evaluated the performance of the EE-GMF method with larger basis sets. Four other basis sets (namely, 6-311+G*, 6-311++G**, cc-pVDZ, and aug-cc-pVDZ) were chosen in this study. The convergence test as a function of the distance threshold was carried out at the HF/6-311+G* level on the large-sized water cluster (140 water molecules; see Figure 1a). HF calculations with 6-311++G**, cc-pVDZ, and aug-cc-pVDZ basis sets were performed on the reduced water cluster (79 water molecules; see Figure 1d) considering the affordable computation cost for the full system QM calculations. The variation of the EE-GMF energy at the HF/6-311+G* level as a function of the distance threshold is shown in Figures 2c and 2d. As can be seen from the figure, the total energies are all nearly converged at λ1 = 6.0 Å (Figure 2c) and λ2 = 6.0 Å (Figure 2d), respectively, for the two-body and three-body QM interactions. For the calculated energy which includes only the two-body QM interactions (see Figure 2c), it converges to −10644.5034 a.u. with two water molecules as a fragment. When the local three-body interactions are included (see Figure 2d), the EE-GMF calculated energy converges around −10644.520 a.u., which is close to the conventional full system energy of −10644.5290 a.u. The EE-GMF energy calculated using λ1 = λ2 = 6.0 Å has a deviation of 6.31 kcal/mol from the full system result (see Table 1). When only one water molecule is taken as a fragment, both of the converged energies, including two- and three-body interactions, show much larger deviations from the full system energy. Here, we show that the

We also investigate the convergence of the EE-GMF energy by varying the two-body and three-body distance threshold values of λ1 and λ2, respectively. Figure 2 shows the calculated total energy of the water cluster with the change of the distance thresholds, with respect to the conventional full system result. As can be seen from Figure 2, the total energies are all nearly converged at λ1 = 6.0 Å (Figure 2a) and λ2 = 6.0 Å (Figure 2b), respectively, for the two-body and three-body QM interactions. The distance threshold λ1 for the two-body QM interaction is fixed at 7.0 Å, when the total energy convergence of the water cluster with the change of the distance threshold λ2 for the three-body QM interaction is tested (see Figure 2b). The total energy, obtained from fragmentation using two water molecules as a fragment, seems to converge more rapidly. The accuracy of the total EE-GMF energy is improved when the fragment size increases from one to two water molecules. For the calculated energy, which includes only the two-body QM interaction (see Figure 2a), it converges to −10641.1197 a.u. with one water molecule as a fragment. When two nearest water molecules are taken as a fragment, the converged energy is −10641.1089 a.u., which is closer to the full system energy of −10641.0851 a.u. When the local three-body QM interaction is also included (see Figure 2b), both of the energies obtained using these two fragment sizes are in better agreement with the full system result. The converged energies are −10641.0794 and −10641.0819 a.u., respectively, for fragmentation with one water molecule as a fragment and two water molecules as a fragment. The absolute errors of 3.58 and 2.01 kcal/mol using the two fragmentation schemes are very small, especially for the fragmentation with two water molecules as a fragment. It is worth noting that the local three-body QM interactions play an important role in reproducing the full system energy using the fragmentation method. The three-body QM interaction using one water molecule as a fragment contributes 25.29 kcal/mol to the total energy, and the three-body QM interaction correction is 16.94 kcal/mol if two water molecules are taken as a fragment. Therefore, the three-body QM interaction has a nonnegligible contribution to the total energy for clusters, as also demonstrated by previous studies.53−56 In addition, the fragmentation method by including many-body interactions up to three-body correction is capable of accurately reproducing full system energy (within the error of 0.01 a.u.). To investigate the environmental polarization effect in the fragmentation method, we also carried out the fragment QM calculation without the presence of the background charges. The results are shown in Figure S1 in the Supporting Information. It clearly shows that the neglect of the environmental polarization effect leads to the converged energies deviating largely from the full system energy, and the total energies also converge much slower. The energy calculated with one water molecule as a fragment begins to converge at λ2 = 6.75 Å (λ1 is fixed at 7.0 Å) without electrostatic embedding, while the energy calculated with two water molecule as a fragment begins to converge at λ2 = 6.50 Å without electrostatic embedding. However, the converged energy using one water molecule as a fragment without electrostatic embedding has a much larger error (14.62 kcal/ mol) with reference to the full system energy, compared to the results obtained using the electrostatic embedding scheme. Moreover, the error of the converged energy (without electrostatic embedding) using two water molecules as a fragment shows a similar magnitude as that of the fragmentation using one water molecule as a fragment using E

DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 1. Total Energies of the Water and Ion-Water Clusters Obtained from Conventional Full System HF Calculations with Different Basis Sets and Deviations of EE-GMF from the Corresponding Full System Energies at the Same Level with the Distance Threshold λ1 = λ2 = 6.0 Å for Water and Ion-Water Clusters basis set

full system energy (a.u.)

fragmentation

6-31G*

Water Cluster System −10641.085051 (1) (2)

6-311+G*

−10644.529007

Table 2. Total Energies of the Water and Ion−Water Clusters Obtained from Conventional Full System B3LYP Calculations with Different Basis Sets and Deviations of EEGMF from the Corresponding Full System Energies at the Same Level with the Distance Threshold λ1=λ2 = 6.0 Å basis set

deviation of EE-GMF (kcal/mol) 6.31 3.14

aug-cc-pVDZ

deviation of EE-GMF (kcal/mol)

cc-pVDZ

−8464.541613

aug-cc-pVDZ

−8465.366270

(2, 2) (n(2.5 Å), 2)

7.27 6.14

−14852.810043

(1, 1) 19.45 (1, 2) 12.55 (2, 1) 15.68 (2, 2) 8.96 (n(2.5 Å), 2) 5.90 Reduced NaCl−Water Cluster from 3 M NaCl Solution Systemb 6-311++G** −8430.618633 (2, 2) 2.42 (n(2.5 Å), 2) 3.02 cc-pVDZ

fragmentation

Water Cluster Systema 6-31G* −10698.540794 (1) 3.39 (2) 2.89 Reduced Water Cluster Systemb 6-311+G* −6040.427095 (2) 1.28 cc-pVDZ −6039.065645 (2) 5.43 aug-cc-pVDZ −6040.096325 (2) 4.97 NaCl−Water Cluster from 3 M NaCl Solution Systema 6-31G* −14913.757499 (2, 2) 2.39 (n(2.5 Å), 2) −2.26 Reduced NaCl−Water Cluster from 3 M NaCl Solution Systemb 6-311+G* −8465.651088 (2, 2) −1.23 (n(2.5 Å), 2) −1.25

a

(1) 18.13 (2) 6.31 Reduced Water Cluster Systemb 6-311++G** −6007.524558 (2) 4.12 cc-pVDZ −6007.102332 (2) 3.77 aug-cc-pVDZ −6006.556007 (2) 5.65 NaCl−Water Cluster from 3 M NaCl Solution Systema 6-31G* −14849.320671 (1, 1) 11.55 (1, 2) 6.65 (2, 1) 8.53 (2, 2) 3.51 (n(2.5 Å), 2) −2.59 6-311+G*

full system energy (a.u.)

−8430.369079

(2, 2) (n(2.5 Å), 2)

(2, 2) 9.34 (n(2.5 Å), 2) 6.24 NaCl−Water Cluster from 6 M NaCl Solution Systema 6-31G* −18432.712785 (2, 2) 5.52 (n(2.5 Å), 2) 4.20 Reduced NaCl−Water Cluster from 6 M NaCl Solution Systemb 6-311+G* −10890.792853 (2, 2) −5.24 (n(2.5 Å), 2) −2.53

3.85 3.15

cc-pVDZ

−10889.901582

(2, 2) (n(2.5 Å), 2)

10.96 6.27

aug-cc-pVDZ

−10890.552749

(2, 2) (n(2.5 Å), 2)

8.62 6.47

−8430.896775

(2, 2) 7.96 (n(2.5 Å), 2) 6.75 NaCl−Water Cluster from 6 M NaCl Solution Systema 6-31G* −18366.354379 (2, 2) 4.33 (n(2.5 Å), 2) 1.51 Reduced NaCl−Water Cluster from 6 M NaCl Solution Systemb 6-311++G** −10853.534302 (2, 2) 1.11 (n(2.5 Å), 2) 1.10

cc-pVDZ

−10853.516035

(2, 2) (n(2.5 Å), 2)

8.79 5.93

aug-cc-pVDZ

−10853.928523

(2, 2) (n(2.5 Å), 2)

9.18 5.57

a The full molecular clusters shown in Figures 1a−c. bThe reduced molecular clusters shown in Figures 1d−f.

pVDZ. As shown in Table 1, the deviation of the EE-GMF energy from the full system energy at the HF/6-311++G** level is 4.12 kcal/mol when two nearest water molecules are taken as a fragment, and the deviations are 3.77 and 5.65 kcal/ mol, respectively, using HF/cc-pVDZ and HF/aug-cc-pVDZ. The EE-GMF approach was also applied at the DFT (namely, B3LYP and M06-2X) and MP2 levels with larger basis sets. As shown in Tables 2 and 3, the total energies of the water cluster calculated by EE-GMF with two water molecules as a fragment are also in good agreement with the full system result at DFT levels, where the largest error is 6.12 kcal/mol in the case of M06-2X/6-311+G* level. For MP2 calculations, we also utilized larger basis sets on the reduced water cluster (Figure 1g). As shown in Table 4, the performance of the EE-GMF method, using two water molecules as a fragment, is also very good at the MP2 level with large basis sets. The largest error of EE-GMF is 6.79 kcal/mol at the MP2/aug-cc-pVDZ level, while for the other two basis sets (6-311+G* and cc-pVDZ), the errors are smaller than 4.90 kcal/mol. Hence, we conclude

a The full molecular clusters shown in Figures 1a−c. bThe reduced molecular clusters shown in Figures 1d−f.

EE-GMF method with the larger basis set can also achieve good convergence and accuracy using the distance thresholds of λ1 = λ2 = 6.0 Å for the two-body and three-body QM interactions, respectively, when two nearest water molecules are grouped into one fragment. The same distance thresholds (λ1 = λ2 = 6.0 Å) are utilized for HF calculations on the reduced water cluster (Figure 1d) with the basis sets of 6-311++G**, cc-pVDZ, and aug-ccF

DOI: 10.1021/acs.jctc.7b00149 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 3. Total Energies of the Water and Ion−Water Clusters Obtained from Conventional Full System M06-2X Calculations with Different Basis Sets and Deviations of EEGMF from the Corresponding Full System Energies at the Same Level with the Distance Threshold (λ1 = λ2 = 6.0 Å) basis set

full system energy (a.u.)

fragmentation

Table 4. Total Energies of the Water and Ion−Water Clusters Obtained from Conventional Full System MP2 Calculations with Different Basis Sets and Deviations of EEGMF from the Corresponding Full System Energies at the Same Level with the Distance Threshold (λ1 = λ2 = 6.0 Å)

deviation of EE-GMF (kcal/mol)

basis set

Water Cluster Systema 6-31G* −10693.646053 (1) 10.37 (2) 5.65 Reduced Water Cluster Systemb 6-311+G* −6037.736597 (2) 6.12 cc-pVDZ −6036.543433 (2) 5.47 NaCl−Water Cluster from 3 M NaCl Solution Systema 6-31G* −14908.500727 (2, 2) 3.26 (n(2.5 Å), 2) 0.50 Reduced NaCl−Water Cluster from 3 M NaCl Solution Systemb 6-311+G* −8462.949262 (2, 2) 8.42 (n(2.5 Å), 2) 6.15 cc-pVDZ

fragmentation

deviation of EE-GMF (kcal/mol)

Reduced Water Cluster Systema 6-31G* −6020.155540 (1) 6.90 (2) 2.89 Reduced Water Cluster Systemb 6-311+G* −3813.037944 (2) 3.50 cc-pVDZ −3812.344181 (2) 4.83 aug-cc-pVDZ −3813.678530 (2) 6.79 Reduced NaCl−Water Cluster from 3 M NaCl Solution Systema 6-31G* −8442.409507 (2, 2) 3.39 (n(2.5 Å), 2) 1.20 Reduced NaCl−Water Cluster from 3 M NaCl Solution Systemb 6-311+G* −5372.917475 (2, 2) 2.49 (n(2.5 Å), 2) −0.85

−8461.955346

(2, 2) 6.53 (n(2.5 Å), 2) 5.10 NaCl−Water Cluster from 6 M NaCl Solution Systema 6-31G* −18428.284049 (2, 2) 6.27 (n(2.5 Å), 2) 4.73 Reduced NaCl−Water Cluster from 6 M NaCl Solution Systemb 6-311+G* −10888.128271 (2, 2) 6.34 (n(2.5 Å), 2) 3.63 cc-pVDZ

full system energy (a.u.)

−10887.305851

(2, 2) (n(2.5 Å), 2)

cc-pVDZ

−5372.324474

aug-cc-pVDZ

−5373.570427

cc-pVDZ

−7008.537902

(2, 2) (n(2.5 Å), 2)

5.93 5.45

aug-cc-pVDZ

−7009.712489

(2, 2) (n(2.5 Å), 2)

12.03 5.20

(2, 2) (n(2.5 Å), 2)

4.21 1.62

(2, 2) 11.66 (n(2.5 Å), 2) 6.56 Reduced NaCl−Water Cluster from 6 M NaCl Solution Systema 6-31G* −10864.501264 (2, 2) 3.48 (n(2.5 Å), 2) 3.45 Reduced NaCl−Water Cluster from 6 M NaCl Solution Systemb 6-311+G* −7009.029723 (2, 2) 5.17 (n(2.5 Å), 2) 4.78

10.41 6.02

a The full molecular clusters shown in Figures 1a−c. bThe reduced molecular clusters shown in Figures 1d−f.

that EE-GMF is an efficient and accurate method for reproducing the full system energy of large-sized water cluster at diverse QM levels. The accuracy of the EE-GMF method can be improved as the fragment size increases. 3.2. Accuracy of the EE-GMF Method for the LargeSized Ion−Water Clusters. The main focus of this study is to investigate the accuracy of the EE-GMF approach in reproducing the full system energy of clusters containing mixed water molecules and ions. We applied the EE-GMF approach to calculate the energy of ion−water cluster from 3 M NaCl solution, which contains 136 water molecules, 8 Na+ cations, and 7 Cl− anions (see Figure 1b). For each ion in the cluster, we combine it with its neighboring water molecules to form an ion−water fragment to treat the charge transfer effect between each ion and its nearby water molecules properly, which is an important factor governing the performance of the EE-GMF method. In this study, we systematically test different fragment sizes for the ion with its neighboring water molecules. Each ion in the cluster can be assigned to its nearest water molecule or several nearest water molecules to form a fragment. For the rest of water molecules in the cluster, one or two water molecules can be taken as a fragment, as we assigned for the pure water cluster. Different fragmentation schemes used for the ion−water cluster are shown in Figure 3. Here, we use “Fragmentation (n1, n2)” to denote the fragmentation scheme, in which each ion is assigned to its n1 nearest water molecules to form a fragment, and for the rest of water molecules, n2 nearest water molecules are taken as a fragment.

a

The reduced molecular clusters shown in Figures 1d−f. reduced molecular clusters shown in Figures 1g−i.

b

The

In a similar manner as the test procedure for the pure water cluster, we first test the total EE-GMF energy convergence as the distance thresholds λ1 and λ2 (for the cutoff of two-body and three-body QM interactions, respectively) vary at the HF/ 6-31G* level, using different fragmentation schemes (namely, Fragmentations (1, 1), (1, 2), (2, 1), (2, 2), and (n(2.5 Å), 2)). In “Fragmentation (n(2.5 Å), 2)”, the neighboring n water molecules within 2.50 Å from the ion center (the distance between the ion and the oxygen atom of water) are assigned to the central ion to form a fragment, and two nearest water molecules are taken as a fragment for the rest water molecules. The case that two ions get close together (