Accurate Prediction of NMR Chemical Shifts in Macromolecular and

Oct 4, 2017 - Researcher's 'Star Wars' Parody Video Produces Robust New Modeling Tool. The “Star Wars” anthology is one of the most famous movie s...
0 downloads 16 Views 695KB Size
Subscriber access provided by LAURENTIAN UNIV

Article

Accurate Prediction of NMR Chemical Shifts in Macromolecular and Condensed-phase Systems with the Generalized Energy-based Fragmentation Method Dongbo Zhao, Ruiheng Song, Wei Li, Jing Ma, Hao Dong, and Shuhua Li J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00380 • Publication Date (Web): 04 Oct 2017 Downloaded from http://pubs.acs.org on October 7, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Accurate Prediction of NMR Chemical Shifts in Macromolecular and Condensed-phase Systems with the Generalized Energy-based Fragmentation Method Dongbo Zhao1,2, Ruiheng Song2, Wei Li1, Jing Ma1, Hao Dong2,*, Shuhua Li1,* 1

Key Laboratory of Mesoscopic Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, School of Chemistry and Chemical Engineering, Nanjing University, 210023, P. R. China; 2 Kuang Yaming Honors School, Nanjing University, 210023, P. R. China. ABSTRACT: The generalized energy-based fragmentation (GEBF) method is extended to allow calculations of nuclear magnetic resonance (NMR) chemical shifts of macromolecular and condensed-phase systems feasible at a low computational cost. In this approach, NMR shielding constants in a large system are evaluated as a linear combination of the corresponding quantities from a series of small ʺelectrostatically embeddedʺ subsystems. Comparison of NMR shielding constants from the GEBF-X [X is an electronic structure method like Hartree-Fock (HF), Density Functional Theory (DFT), ···] method with those from the conventional quantum chemistry method for two representative systems verifies that the GEBF approach can reproduce the results of the conventional quantum chemistry method very well. This procedure has further been applied to compute NMR shielding constants of a large foldamer and a supramolecular aggregate, and the 15N shielding constant for CH3CN in the CHCl3 solvent. For the former two systems, the predicted 1H chemical shifts are in good agreement with the experimental data. For the CH3CN/CHCl3 solution, the 15N shielding constant of CH3CN is evaluated as the ensemble average of up to 200 sufficiently large CH3CN/CHCl3 clusters from either classical or QM/MM (quantum mechanics/molecular mechanics) molecular dynamics (MD) simulations. Our results reveal that the gas-to-solution shift of 15N (from an isolated CH3CN to the CH3CN/CHCl3 solution) based on PM6-DH+/MM MD simulation is in good accord with the experimental value, outperforming those based on classical MD simulation and the previous IEF-PCM (polarizable continuum model using integral equation formalism) study. This study unravels that the generation of representative liquid structures is critical in evaluating the NMR shielding constants of condensed phase systems.

1. INTRODUCTION Nuclear magnetic resonance (NMR) spectroscopy has an unparalleled impact on the structure determination1 and dynamics2 of various chemical and biological systems. Nonetheless, accurate prediction of NMR chemical shifts with conventional quantum mechanics (QM) methods in macromolecular and condensedphase systems is still a daunting task, as the computational cost skyrockets (at least the fourth power) with molecular size. To overcome this difficulty, lower or linear-scaling methods have gained considerable attraction during the past decades. Progresses revolving round this topic have been summarized in a recent review (Ref. 3). A variety of fragment-based methods4−21 have been applied to compute the NMR shielding constants, which include fragment molecular orbitals (FMO),4−6 adjustable density matrix assembler (ADMA),7−9 combined fragmentation methods (CFM),10 and so on.11−21 All the aforementioned methods share some common characteristics that, in general, they can nicely reproduce the conventional computations of chemical shielding constants for a large molecule, such as a protein where intramolecular interactions dominate, either in gas-phase4,5−8,10−15,20,21 or implicit Poisson–Boltzmann15 and polarization continuum8 solvent models; in addition, some studies have made comparisons with experiment data.4,7−9,14−20 Nevertheless, only quite a few works take explicit solvent effect into consideration.9,20,21 Actually, it is well documented that

Figure 1. Schematic diagram of the GEBF method for computing the NMR shielding constants of a solute in a solution system. The NMR shielding constants in a large system (left) could be approximately computed as a linear combination of the corresponding quantities from a series of small ʺelectrostatically embeddedʺ subsystems (right). For the solute molecule M, only those subsystems containing M are needed to compute the NMR shielding constants in GEBF calculations. The background point charges for atoms outside a subsystem are neglected for clarity.

The combined quantum mechanics/molecular mechanics (QM/MM) method provides an alternative choice to the full QM treatment for evaluating NMR chemical shieldings.26 However, polarization effect of the environment cannot be fully addressed using fixed partial charges in conventional force field methods.

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

For example, Steinmann et al.27 has revealed that ab initio-derived polarizable embedding potentials are superior to the conventional force fields, reducing the size of the QM region at least 2 Å smaller without compromising the quality. But how large should a QM region sufficiently be? Ochsenfeld et al.28 has proven that QM regions have to contain neighboring atoms within at least 610 Å from the considered nucleus/nuclei to reach size convergence. In this regard, linear or lower-scaling methods are good choices for treating such a large QM region.

Page 2 of 10

mentation details of the GEBF-NMR method, and the validation result of the GEBF method as well as the details of the classical and QM/MM MD simulations for systems studied. In Section 3, we apply the GEBF method to three realistic systems, and compare the calculated results with the experimental data. In Section 4, a brief summary is included. 2. METHODOLOGY 2.1 The GEBF method Within the framework of the GEBF method, the ground-state energy of a large system can be evaluated as a linear combination of a series of ″electrostatically embedded″ subsystems, as expressed in eq 1. M M QQ ~ Etot = ∑ Cm Em − [(∑ Cm ) − 1]∑ ∑ A B m m A B > A RAB

(a)

(b)

Figure 2. Two large molecular systems used for validation of the GEBF-NMR protocol. Left (a): a positively (+1) charged Trp-cage mini-protein (304 atoms) in gas phase. The backbone is shown in cartoon model, and side chains are in stick model. Right (b): a truncated cluster (266 atoms) from a snapshot of QM/MM MD simulations, with the central solute molecule (CH3CN) represented in vdW surface and surrounding CHCl3 molecules (within a radius of 10 Å from the solute molecule) in magenta.

In the present work, we calculated the NMR shielding constants of macromolecular and condensed-phase systems within the framework of the generalized energy-based fragmentation (GEBF) method,29−31 denoted as GEBF-NMR. The GEBF method has been demonstrated to be very successful for calculating ground-state energies and some molecular properties (such as dipole moment and polarizability) of very large systems, due to its unique way of constructing “electrostatically embedded” subsystems. In this work, we will first validate the accuracy of the GEBF method in computing the NMR shielding constants for two representative systems (system 1 and system 2). Here system 1 is a mini-protein Trp-cage (PDB ID: 1L2Y32, 304 atoms) in gas phase and system 2 is a set of ten CH3CN/CHCl3 molecular clusters33 (truncated from MD snapshots, where CHCl3 is a solvent as shown in Figure 2). Then, we applied the GEBF-NMR protocol to three real chemical systems, including an aromatic oligoamide βsheet foldamer34 (system 3), a supramolecular aggregate35 (system 4), and a solution system 5 (including 200 CH3CN/CHCl3 clusters, different from system 2). For this solution system, up to 200 large CH3CN/CHCl3 molecular clusters generated from classical molecular dynamics (MD) or QM/MM MD simulations are utilized to obtain the NMR shielding constants for the solute molecule (CH3CN). By comparing the calculated results with the corresponding experimental data, we verified that the GEBF approach can achieve high accuracy in predicting the NMR shielding constants of both the macromolecular and the condensedphase systems. This study reveals that the configuration sampling of the solution system with the QM/MM method can give much better results than that with the force field method. In addition, for fragment-based methods (like the GEBF method) within the density functional theory (DFT) framework, we found that the high quality of the integration grid for the exchange-correlation potential should be employed in NMR calculations of all subsystems. The remainder of this paper is organized as follows. In Section 2, we describe the theoretical foundations of GEBF, the imple-

(1)

where Ẽm stands for the total energy of the mth embedded subsystem (with electrostatic interactions between background point charges), Cm represents the coefficient of the mth subsystem, QA is the net atomic charge on atom A, and M is the total number of subsystems. The determination of Cm and QA has been discussed in our previous work30. Here, we will just outline some key points of the procedure. Let us take a molecular cluster CH3CN/CHCl3 as an example. Here each single molecule (CH3CN or CHCl3) is considered as a fragment. Next, we build primitive subsystems centered at each fragment by a distance threshold. The coefficients for these primitive subsystems are set to be +1. Then, smaller derivative subsystems have to be generated since primitive subsystems are generally overlapping with each other. The coefficients of these derivative subsystems can be automatically determined according to the principle of inclusion and exclusion. The same strategy is also applied to covalently bonded systems. For these systems, the first step is to partition the whole system into many fragments by cutting single bonds. In these cases, additional H atoms should be added in the broken bond in constructing subsystems so that each subsystem has a closed-shell electronic structure. To further reduce the computational cost without much loss of the accuracy, we introduce another parameter (the maximum number of fragments) to control the size of a subsystem. Finally, we obtain the atomic charges for all atoms from calculations of all the primitive subsystems. Since small subsystems can be calculated with existing quantum chemistry programs, the GEBF method allows electronic structure calculations of very large systems feasible. Recent developments and applications of the GEBF method can be found elsewhere.36−41 The first-order derivatives of the total energy with respect to nuclear displacements (or other quantities) of a large system can also be obtained in a similar way as the total energy described above. Thus, the GEBF method is applicable to ab initio geometry optimizations of very large systems.42 Furthermore, the secondorder or even higher derivatives of the total energy (or the corresponding molecular properties) can be approximately obtained using the following formula.42 M ~ Ωtot ≈ ∑CmΩm

(Ω = µi ,αij ,L)

(2)

m

where

Ωtot

represents the total molecular property and

~ Ωm

cor-

responds to the same quantity in small subsystems. In this work, the above formula (eq 2) will be used to evaluate the chemical shielding tensor σ (eq 3) on each nucleus in a large system from the corresponding quantities of small subsystems. Here σ for a given nucleus (N) is defined as a second-order de-

ACS Paragon Plus Environment

Page 3 of 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

rivative of the total energy with respect to external magnetic field (B) and magnetic moment (m),

σ

N ij

∂2E = ∂Bi ∂m N j

(3)

Normally, only the isotropic shielding constant, defined as onethird of the trace of the atomic shielding tensor (eq 4), is considered.

1 3

N σ iso = tr[σ ]

(4)

If the shielding constant of a reference compound is available, the chemical shift (relative to this reference compound) could be calculated as follows, N δ N = σ iso − σ ref

(5)

During the implementation of the GEBF method, a large system will be divided into smaller fragments automatically or manually. Based on these fragments, primitive and derivative subsystems can be automatically constructed for subsequent calculations. Two key parameters are used to control the accuracy and the computational cost of a GEBF calculation. A distance threshold (4 Å) is chosen for the construction of primitive subsystems, and the largest number of the fragments in a primitive subsystem is set to be 8 in all GEBF calculations. It should be mentioned that the GEBF calculation of a large system can be further simplified if only the chemical shifts of a subset of atoms in a large system are needed. According to eq 2, the chemical shift of the atom M in a large molecule is taken as a linear combination of the chemical shift on the corresponding atom in those subsystems containing M as real atoms, as schematically shown in Figure 1. For a solution system, we are usually interested in only chemical shifts of atoms in the solute molecule. Thus, the GEBF method should be a very cost-effective method for evaluating the chemical shifts of condensed phase systems. In addition, our GEBF-NMR computations for all subsystems are efficiently parallelized with the MPI (message passing interface) technique across different computing nodes (since GEBF-NMR calculations on different subsystems can be done on different nodes). All quantum chemical computations were performed with the GEBF method as implemented in our LSQC (low scaling quantum chemistry) package,43 interfaced with the Gaussian 09 package44 for NMR shielding constants calculations on subsystems.

2.2 NMR Shielding Constant Calculations Density functional theory calculations with tight self-consistent field (SCF) convergence criteria and ultrafine integration grids (unless otherwise stated) were carried out to calculate the NMR shielding constants for systems (and subsystems in GEBF calculations) under study. The gauge invariance was ensured by using the gauge-including atomic orbitals (GIAOs).45−49 For all systems (except systems 2 and 5) under study, their optimized structures were obtained with the GEBF-B97-250 method at the 6-31G(d,p)51 level in gas phase, and the NMR shielding constants were evaluated with the same method at the locally dense basis set pcSseg-2, which was demonstrated to be suitable for NMR shielding constants calculations.52 In GEBF calculations, for molecular clusters 2 and 5, each molecule is taken as a fragment. For other systems 1, 3 and 4, the fragmentation schemes of partitioning each system into different fragments are shown in the Supporting Information (SI), Figures S1-S3, respectively. No solvation effect was considered for systems 1, 3 and 4. The structure of the reference com-

pound, tetramethylsilane (TMS), and its NMR shielding constants were calculated at the B97-2/pcSseg-2 level of theory. For the solution system 5, representative snapshot clusters were extracted from the trajectory of well-equilibrated MD simulations (as described below), each of which incorporates a solute molecule (CH3CN) and its surrounding solvent molecules (CHCl3) within 10 Å from this solute molecule (see Figure 2b). For each snapshot structure, we have checked the convergence of the calculated 15N chemical shielding constants over the size of CH3CN/CHCl3 cluster (as shown in Figure S4), and found that a cutoff distance of 10 Å can give almost convergent results. It should be mentioned that similar results were reported in the literature.28 In subsequent GEBF calculations, the CH3CN/CHCl3 clusters were calculated at the GEBF-B97-2/pcSseg-2 theory level. To calibrate the GEBF-NMR protocol, a sequence of 10 clusters was used, and the calculated 15N NMR values for the CH3CN solute is compared with the reference data from conventional QM calculations (B97-2/pcSseg-2). To evaluate the 15N NMR value in system 5, ensemble average from up to 200 CH3CN/CHCl3 clusters was used and a convergent result is harvested. For an isolated CH3CN molecule, its optimized structure and the 15N NMR value were computed at the B97-2/pcSseg-2 level of theory. Then the difference between the 15N shielding constant of an isolated CH3CN molecule and that in the CH3CN/CHCl3 solution was calculated to get the gas-to-solution chemical shift of 15N. For the solution phase, the error bars and uncertainties were calculated as standard error of the mean (SEM) of the shielding constants for a set of clusters. 2.3 Molecular Dynamics Simulations for the CH3CN/CHCl3 Solution For the CH3CN/CHCl3 solution, both classical and QM/MM MD simulations were used to provide representative snapshot structures, from which up to 200 CH3CN/CHCl3 clusters can be extracted for NMR calculations. The general Amber force field (GAFF) method53 and restraint electrostatic potential (RESP)54,55 were used to generate force field parameters for CH3CN and CHCl3. In classical MD simulations, an NPT (P = 1 atm, T = 300 K) ensemble was used.56−58 Periodic boundary conditions were applied, with the Particle Mesh Ewald method used for treating long-range electrostatic interactions. All bonds involving hydrogen atoms were constrained and an integration time step of 2 fs was used. The trajectory was accumulated for 150 ns, where the first 100 ns was for equilibration, and the last 50 ns was used to generate clusters for subsequent NMR computations. All classical MD simulations were carried out with AMBER16.59 In QM/MM MD simulations, the solute molecule was treated as the QM region, and all the solvent molecules were in the MM region. For the QM region, either the semi-empirical PM6-DH+60 or the BLYP-D261−63 method was used. The former was carried out with AMBER16,59 and a trajectory of 35 ns (the last 10 ns was used for production) was generated. The latter was run with the QUICKSTEP module in CP2K64,65 for 15 ps (the last 5 ps was for production). In BLYP-D2 simulations, the triple-ζ basis set with two polarization functions,66 and the Goedecker, Teter, and Hutter (GTH) pseudopotential67,68 were used. An energy cutoff of 320 Ry was employed in the plane wave representation of the density. As only a CH3CN molecule was treated with QM method, the box size for the QM region is 8×8×8 Å3. The temperature was controlled at 300 K under an NVT ensemble, and the time step of 0.5 fs was used. 2.4 Validation of the GEBF method

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

To validate the reliability of the implemented GEBF-NMR protocol, systems 1 and 2 were selected as representative examples of large systems and liquid phase structures (Figure 2), respectively. For these two systems, we will compare the calculated chemical shifts obtained with the GEBF-DFT (or GEBF-HF) method and the conventional DFT (or HF) data. Three statistical parameters, root-mean-squared deviation (RMSD), mean signed error (MSE) and mean absolute error (MAE) were used to characterize the performance of our GEBF-NMR method. Table 1. RMSD and Maximum Error of Shielding Constants (in ppm) from GEBF-NMR Computations of a Positively Charged Trp-cage Mini-protein, Compared to Those Obtained from Conventional QM Methods. The pcSseg-2 Basis Set is Used for All Calculations.

HF

all H C N O

CAM-B3LYP

RMSD

max error

0.16 0.01 0.05 0.10 0.51

1.66 0.05 0.20 0.26 1.66

M06-2X

RMSD

max error

RMSD

max error

0.36 0.04 0.19 0.48 1.02

2.43 0.12 0.66 1.24 2.43

0.23 0.02 0.07 0.15 0.70

2.14 0.06 0.18 0.41 2.14

For system 1 (a Trp-cage mini-protein), the deviations between the calculated shielding constants obtained with the GEBF-X and conventional X methods (X = HF69, CAM-B3LYP70, and M062X71) and Jensen's pcSseg-2 basis set are listed in Table 1. It is unequivocal that the GEBF-QM method can well reproduce the conventional QM method in the prediction of NMR shielding constants, with RMSDs being less than 1.0 ppm for four different atoms at three different theoretical levels considered. It is worth pointing out that the accuracy of the GEBF-QM method is independent of the levels of theory and basis sets adopted, as seen from the comparison of GEBF-NMR results with conventional NMR results in Table 1 and Table S1 (more information is given). 15

Table 2. The N Shielding Constants in 10 CH3CN/CHCl3 Clusters Calculated with the GEBF-B97-2/pcSseg-2 and B97-2/pcSseg2 methods. The Results Obtained with Different Integration Grids are Shown for Comparison.

Ultrafine Grida

Fine Gridb

Clusters

GEBFB97-2

B97-2

Diff.

1 2 3 4 5 6 7 8 9 10 MSE MAE

–34.8 –21.3 –23.8 –10.4 –16.7 –18.2 –34.2 –36.3 –41.8 –16.8

–34.4 –21.3 –23.5 –10.1 –16.3 –18.0 –34.2 –36.6 –40.8 –17.7

–0.4 0.0 –0.3 –0.3 –0.4 –0.2 0.0 0.3 –1.0 0.9 –0.2 0.4

GEBFB97-2

B97-2

–28.8 –16.8 –18.6 –3.2 –11.7 –12.4 –29.9 –33.6 –33.2 –11.4

–34.4 –21.2 –23.4 –10.1 –16.2 –18.0 –34.1 –34.4 –40.6 –21.2

Diff. 5.6 4.4 4.8 6.9 4.5 5.6 4.2 0.8 7.4 9.8 5.4 5.4

To evaluate the performance of the GEBF-QM method in predicting NMR chemical shifts of the solution phase, we take 10 CH3CN/CHCl3 clusters from the trajectory of QM/MM MD simulations (where the QM region is treated with PM6-DH+, as described above) and calculate the 15N NMR shielding constants with the GEBF-B97-2 and conventional B97-2 methods at the pcSseg-2 basis set. The coordinates of these 10 clusters can be found in the SI. The results obtained with fine and ultrafine integration grids are collected in Table 2. Surprisingly, the results obtained with GEBF-B97-2 strongly depend on the accuracy of the integration grid adopted, but the results from the conventional B97-2 method are much less affected. Once an ultrafine integration grid is used, the GEBF-B97-2 NMR shielding constants are very close to the results from B97-2 calculations, with the maximum deviation of 1 ppm. However, this is not the case for a fine integration grid as shown in Table 2. For example, the deviation is 9.8 ppm for cluster 10. It is found that the 15N shielding constant from the subsystem centered at the CH3CN molecule is ‒17.6 ppm, being very close to the value from the whole system (‒17.7 ppm). Nevertheless, the 15N shielding constant from the GEBF approach is only ‒11.4 ppm, indicating that the error accumulation occurs in the GEBF approach when a fine integration grid is employed in NMR calculations of subsystems. With an ultrafine integration grid, the GEBF-QM method is demonstrated to work well for predicting the NMR chemical shift of condensed phase systems. This issue has also been reported by Beran et al. for NMR studies of molecular crystals with a manybody expansion fragment approach.72,73 In addition, with the same computing resources [Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50 GHz 24 CPU cores with 128 GB memory], our calculations show that for system 2 the conventional NMR calculation takes 9.3 hours and the GEBF-NMR calculation needs only 2.5 hours. Hence, the GEBF approach is almost three times faster than its conventional counterpart. One can expect that for larger systems the GEBF approach is even more cost-effective because the computational cost of the GEBF approach scales linearly with the system size. 3. RESULTS AND DISCUSSION 3.1 Macromolecules (A) Aromatic Oligoamide β-Sheet Foldamer. System 3 is a five-stranded artificial β-sheet with 4,6-dinitro-1,3phenylenediamine as linkage building blocks, stabilized by intramolecular aromatic π-π stacking interactions. We took its crystal structure34 as the initial structure and fully optimized it at the GEBF-B97-2/6-31G(d,p) level. The optimized structure is shown in Figure 3. Four types of drastically different H atoms are marked with Hext, Hint, H(NH), and H(NHCO), respectively. Chemical shifts (with TMS as the reference compound) were computed at the GEBF-B97-2/pcSseg-2 (10812 basis functions) level and the mean chemical shifts of each type of H atoms, are displayed in Table 3 in comparison with the experimental data obtained from the CHCl3 solution.34

a Ultrafine grid corresponds to a pruned (99, 590) grid with 99 radial shells and 590 angular points per shell.

b

Page 4 of 10

Fine grid is a pruned (75, 302) grid.

ACS Paragon Plus Environment

Page 5 of 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 3. Schematic representation of a multi-turn aromatic oligoamide β-sheet foldamer (system 3, 414 atoms). Four types of H atoms (Hext, Hint, H(NH) and H(CONH)) are marked. Some hydrogen atoms are omitted for clarity. Color code: H in white, C in golden, N in blue, and O in red. It is clear that the calculated 1H NMR chemical shifts from the GEBF-B97-2 method are in good agreement with those obtained from experiment, as evidenced by the relatively small MSE and MAE values, 0.7 and 0.8 ppm, respectively. The remaining discrepancy is probably due to the neglect of the rovibrational effect and ensemble average. Moreover, we have checked the impact of the solvent CHCl3 by using the implicit conductor polarizable continuum model (CPCM), and found that the solvent effect leads to negligible changes for the NMR chemical shifts of H atoms (results not shown). Table 3. Comparison between the Calculated 1H Chemical Shifts (in ppm, at the GEBF-B97-2/pcSseg-2 Level) and the Experimental Data for a β-Sheet Foldamer. 1

Type Hint Hext H(NH) H(NHCO) MSE MAE

H Chemical Shifts (in ppm)

GEBF-B97-2

Expt.

4.5 10.1 10.3 10.4

4.6 9.3 8.9 9.5

Dev. (ppm) –0.1 0.8 1.4 0.9 0.7 0.8

(B) Supramolecular Aggregate. System 4 is a supramolecular aggregate CA3·M3 (CA means cyanuric acid and M is melamine) rosette, stabilized by 18 hydrogen bonds. Its optimized structure obtained at the GEBF-B97-2/6-31G(d,p) level is shown in Figure 4.

Figure 4. A supramolecular aggregate of the CA3·M3 rosette (system 4, 252 atoms), where CA means cyanuric acid and M is melamine. In our case, the building blocks are the CA and M derivatives. The target structural motif is stabilized by a total of 18 hydrogen bonds shown in pink lines. Characters a-h mark different hydrogen atoms with NMR activity. Color code: H in white, C in golden, N in blue, and O in red. Table 4. Comparison between the Calculated 1H Chemical Shifts (in ppm, at the GEBF-B97-2/pcSseg-2 Level) and the Experimental Data for a Complex CA3M3. 1

Type

H Chemical Shifts (in ppm)

GEBF-B97-2

Expt.

Dev. (ppm)

Ha 7.6 7.1 0.5 Hb 10.3 9.4 0.9 Hc 8.5 7.7 0.7 Hd 7.8 7.4 0.4 He 1.5 1.4 0.1 Hf 1.2 1.0 0.2 Hg 2.4 2.2 0.2 Hh 14.8 13.5 1.3 MSE 0.6 MAE 0.6 The 1H chemical shifts obtained at the GEBF-B97-2/pcSseg-2 level, together with the experimental data measured in CDCl3 solution, are collected in Table 4. The GEBF-B97-2 method was found to give satisfactory predictions for

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 10

(a) (b) Figure 5. The calculated gas-to-solution shift of 15N NMR in system 5. (a) Convergence check of the number of clusters used for the GEBF-NMR calculations, where the structures were generated by using three different protocols of MD simulations. (b) Distributions of the calculated shifts by using clusters from three protocols, where red, green and blue histograms are for classical MD, PM6-DH+/MM MD and BLYP-D2/MM MD, respectively. 1 H NMR chemical shifts, with both MSE and MAE values of 0.6 ppm, respectively. It is worthy to point out that among different H atoms, the largest change of 1H chemical shifts from an isolated molecule to the supramolecular aggregate occurs in Hh of cyanuric acid (from ~7.0 ppm (results not shown) to 14.8 ppm). This huge downfield shift indicates that the self-assembly process take place in solution. Therefore, our results corroborate the experimental finding that the hetero-hexamer aggregate is predominant in solution. Overall, we have demonstrated that the GEBF-B97-2 method is capable of reproducing the measured 1H NMR chemical shifts in large chemical systems. Combined with properly chosen DFT methods, the GEBF method is expected to be a powerful tool for computing the NMR chemical shifts of large complicated systems.

3.2 The CH3CN/CHCl3 Solution The prediction of NMR shielding constants for molecules in the solution phase is not an easy job,33 because the NMR shielding constants are dramatically affected by the environment. Here we use the GEBF-B97-2/pcSseg-2 method to calculate the 15N shielding constants in the CH3CN/CHCl3 solution (system 5). The 15N shielding constants of CH3CN in the CHCl3 solvent is evaluated as the average of up to 200 CH3CN/CHCl3 clusters extracted from the trajectory of well-equilibrated MD simulations. Here three different MD simulations were employed to generate the trajectories: classical MD, PM6-DH+/MM MD, and BLYP-D2/MM MD. The details of MD simulations were described in the preceding section. In Figure 5a, we have plotted the calculated gas-to-solution shift of 15N of CH3CN in system 5 against the number of CH3CN/CHCl3 clusters. First, one can see that the calculated gasto-solution shift of 15N from three different MD simulations is almost converged when more than 120 clusters were used to calculate the ensemble average. As shown in Figure 5b, the calculated shifts of 15N from different MD simulations exhibit a wide distribution. Second, the calculated gas-to-solution shifts of 15N from QM/MM MD simulations are much closer to the experimental value (−6.2 ppm),74 with a deviation of 2-4 ppm, while the result from classical MD simulations underestimates the shift by ~8 ppm (Figure 5a and Table 5). Therefore, the QM/MM treatment of the CH3CN/CHCl3 clusters is critical for reasonably de

scribing the liquid environment, and the satisfactory prediction of NMR shielding constants. Table 5. Comparison of Gas-to-solvent Shift (in ppm) for 15N in CH3CN in the CHCl3 Solvent Calculated with the GEBFB97-2/pcSseg-2 Method. Experimental Data and the Result from IEF-PCM Calculations are also Collected for Comparison. Method a

exp. IEF-PCMb classical MD PM6-DH+/MM MD BLYP-D2/MM MD

δ(15N) −6.2 −17.7 −14.2 ± 0.8 −4.0 ± 0.7 −2.1 ± 0.6

a

Ref. 74.

b

Ref. 33, IEF-PCM-B3LYP/6-311+G(d,p).

With regard to the QM method, the semi-empirical PM6-DH+ results seem to have a slightly better agreement with the experimental data than the more expensive BLYP-D2 results. A possible explanation on the less accuracy of BLYP-D2 results is that QM/MM MD with BLYP-D2 can only sample much smaller configuration space than PM6-DH+, due to its computational cost. For BLYP-D2, representative structures are generated from 15 ps MD simulations, whereas for PM6-DH+ structures are from 35 ns MD simulations. To test this hypothesis, another QM/MM MD simulations with PM6-DH+ were run for only 15 ps. The subsequent GEBF-NMR calculations show that the average result from PM6-DH+ structures has a larger deviation to the experiment result than that from BLYP-D2 structures (Figure S5), if both structures are from 15 ps MD simulations. Therefore, both the accuracy of the QM method and the efficiency of MD simulations (for configurational sampling) play important roles in generating suitable solute/solvent clusters for NMR calculations. Our calculations show that QM/MM MD with PM6-DH+ is a good choice for generating snapshot structures for subsequent NMR calculations. We should mention that the accuracy of the GEBF-QM method for the prediction of NMR shielding constants in condensed-phase systems may be further improved by adopting more accurate electronic structure methods (for NMR calculations) and

ACS Paragon Plus Environment

Page 7 of 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

better QM/MM methods (for sampling the solute/solvent cluster structures). We also compared the results from this work with that from other theoretical method. The results are collected in Table 5. Tomasi et al.33 proposed an IEF-PCM (polarizable continuum model using integral equation formalism) method for the same system. In their work, a small cluster with the solute surrounded by solvent molecules in the first coordination shell (with a buffer distance of 3~4 Å) was calculated with the B3LYP/6-311+G(d,p) method, and PCM was used to treat the remaining bulk solvent effect. This protocol leads to a deviation from the experiment of 11.5 ppm.33 As suggested by P. Bouř and M. Dračínský, standard dielectric continuum solvent models may not be sufficiently accurate for some solution systems.75 Our present study demonstrated that the GEBF-B97-2 method, combined with much larger CH3CN/CHCl3 clusters from QM/MM MD simulations, can give fairly good descriptions on the 15N NMR shielding constants for the CH3CN molecule in the CHCl3 solvent. If much longer QM/MM MD simulations with more accurate DFT methods are employed to generate representative liquid snapshot structures, one would expect that the GEBF-B97-2 method based on these snapshot structures lead to even more accurate results. 4. CONCLUSIONS The GEBF approach is extended to predict the NMR shielding constants for macromolecules and condensed phase systems. Comparison of NMR shielding constants from the GEBF-QM method with those from the conventional QM method for two representative systems demonstrates that the GEBF-QM approach can reproduce the results of the conventional QM approach with a much less computational cost. This procedure has been applied to compute NMR shielding constants of a large foldamer and a supramolecular aggregate, and the 15N shielding constant for CH3CN in the CHCl3 solvent. For the foldamer and the supramolecular aggregate, the calculated 1H chemical shifts from GEBFB97-2/pcSseg-2 calculations are in very good agreement with the experimental data. For the CH3CN/CHCl3 solution, we constructed up to 200 sufficiently large CH3CN/CHCl3 clusters from classical MD and QM/MM MD simulations. The ensemble average of the calculated 15N shielding constants for CH3CN in these clusters (with the GEBF-B97-2/pcSseg-2 method) is used to evaluate the gas-to-solution shift of 15N in CH3CN. Our results show that the calculated gas-to-solution shift of 15N from QM/MM MD simulation based on semi-empirical PM6-DH+ is quite consistent with the experimentally measured value, much more accurate than those from classical MD simulation and the previous IEF-PCM study. The result from QM/MM MD simulation based on BLYPD2 is slightly less accurate than that based on PM6-DH+, due to the fact that at the BLYP-D2 level representative liquid structures are obtained in a shorter simulation time. This study reveals that the generation of representative liquid structures (for the ensemble average) plays a critical role in evaluating the NMR shielding constants of condensed phase systems. In addition, our calculations also suggest that for various solute/solvent clusters accurate predictions of NMR shielding constants with the GEBF-DFT method require an ultrafine integration grid for DFT calculations. It is anticipated that the same strategy should be used for other fragment-based methods if they are used to predict NMR shielding constants of condensed phase systems. In the future study, the GEBF method with more advanced electronic structure methods and ab initio MD simulations for generating the trajectory will be employed to investigate the NMR properties of more complicated and challenging chemical systems.

ASSOCIATED CONTENT Supporting Information. Shielding constant convergence over the QM region, basis sets on the shielding constants, comparison of PM6-DH+ and BLYP-D2 results with the same MD trajectory length and coordinates of ten CH3CN/CHCl3 clusters. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author *[email protected]; *[email protected].

ACKNOWLEDGMENT This work was supported by the National Natural Science Foundation of China (Grants No. 21333004, NO. 21673110, NO. 21361140376), the Natural Science Foundation of JiangSu Province (Grant No. SBK2015041570), and the ″Fundamental Research Funds for the Central Universities″. Part of the calculations were performed using computational resources on an IBM Blade cluster system from the High Performance Computing Center (HPCC) of Nanjing University and the Shenzhen Supercomputer Center (SSC, China).

REFERENCES (1) Wylie, B. J.; Sperling, L. J.; Nieuwkoop, A. J.; Franks, W. T.; Oldfield, E.; Rienstra, C. M. Ultrahigh resolution protein structures using NMR chemical shift tensors, Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 16974–16979. (2) Robustelli, P.; Stafford, K. A.; Palmer II, A. G. Interpreting Protein Structural Dynamics from NMR Chemical Shifts, J. Am. Chem. Soc. 2012, 134, 6365−6374. (3) Collins, M. A.; Bettens, R. P. A. Energy-Based Molecular Fragmentation Methods, Chem. Rev. 2015, 115, 5607–5642. (4) Gao, Q.; Yokojimab, S.; Kohnoc, T.; Ishidad, T.; Fedorovd, D. G.; Kitaurad, K.; Fujihiraa, M.; Nakamura, S. Ab initio NMR chemical shift calculations on proteins using fragment molecular orbitals with electrostatic environment, Chem. Phys. Lett. 2007, 445, 331–339. (5) Gao, Q.; Yokojima, S.; Fedorov, D. G.; Kitaura, K.; Sakurai, M.; Nakamura, S. Fragment-Molecular-Orbital-Method-Based ab Initio NMR Chemical-Shift Calculations for Large Molecular Systems, J. Chem. Theory Comput. 2010, 6, 1428–1444. (6) Amos, R.; Kobayashi, R. Ab Initio NMR Chemical Shift Calculations Using Fragment Molecular Orbitals and Locally Dense Basis Sets, J. Phys. Chem. A 2016, 120, 8907–8915. (7) Frank, A.; Onila, I.; Möller, H. M.; Exner, T. E. Toward the quantum chemical calculation of nuclear magnetic resonance chemical shifts of proteins, Proteins, 2011, 79, 2189−2202. (8) Frank, A.; Möller, H. M.; Exner, T. E. Toward the Quantum Chemical Calculation of NMR Chemical Shifts of Proteins. 2. Level of Theory, Basis Set, and Solvents Model Dependence, J. Chem. Theory Comput. 2012, 8, 1480−1492. (9) Exner, T. E.; Frank, A.; Onila, I.; Möller, H. M. Toward the Quantum Chemical Calculation of NMR Chemical Shifts of Proteins. 3. Conformational Sampling and Explicit Solvents Model, J. Chem. Theory Comput. 2012, 8, 4818–4827. (10) Tan, H.-J.; Bettens, R. P. A. Ab initio NMR chemical-shift calculations based on the combined fragmentation method, Phys. Chem. Chem. Phys. 2013, 15, 7541–7547. (11) Lee, A. M.; Bettens, R. P. A. First Principles NMR Calculations by Fragmentation, J. Phys. Chem. A 2007, 111, 5111–5115. (12) Reid, D. M.; Kobayashi, R.; Collins, M. A. Systematic Study of Locally Dense Basis Sets for NMR Shielding Constants, J. Chem. Theory Comput. 2014, 10, 146–152. (13) Reid, D. M.; Collins, M. A. Calculating nuclear magnetic resonance shieldings using systematic molecular fragmentation by annihilation, Phys. Chem. Chem. Phys. 2015, 17, 5314–5320.

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) He, X.; Wang, B.; Merz Jr., K. M. Protein NMR Chemical Shift Calculations Based on the Automated Fragmentation QM/MM Approach, J. Phys. Chem. B 2009, 113, 10380–10388. (15) Zhu, T.; He, X.; J. Zhang, Z. H. Fragment density functional theory calculation of NMR chemical shifts for proteins with implicit solvation, Phys. Chem. Chem. Phys. 2012, 14, 7837–7845. (16) Zhu, T.; Zhang, J. Z. H.; He, X. Correction of erroneously packed protein's side chains in the NMR structure based on ab initio chemical shift calculations, Phys. Chem. Chem. Phys. 2014, 16, 18163–18169. (17) Swails, J.; Zhu, T.; He, X.; Case, D. A. AFNMR: automated fragmentation quantum mechanical calculation of NMR chemical shifts for biomolecules, J. Biomol. NMR 2015, 63, 125–139. (18) Jin, X.; Zhu, T.; Zhang, J. Z. H.; He, X. A systematic study on RNA NMR chemical shift calculation based on the automated fragmentation QM/MM approach, RSC Adv. 2016, 6, 108590–108602. (19) Steinmann, C.; Bratholm, L. A.; Olsen, J. M. H.; Kongsted, J. Automated Fragmentation Polarizable Embedding Density Functional Theory (PE-DFT) Calculations of Nuclear Magnetic Resonance (NMR) Shielding Constants of Proteins with Application to Chemical Shift Predictions, J. Chem. Theory Comput. 2017, 13, 525–536. (20) Jose, K. V. J.; Raghavachari, K. Fragment-Based Approach for the Evaluation of NMR Chemical Shifts for Large Biomolecules Incorporating the Effects of the Solvent Environment, J. Chem. Theory Comput. 2017, 13, 1147–1158. (21) Ochsenfeld, C.; Kussmann, J.; Koziol, F. Ab Initio NMR Spectra for Molecular Systems with a Thousand and More Atoms: A Linear-Scaling Method, Angew. Chem. Int. Ed. 2004, 43, 4485–4489. (22) von Philipsborn, W.; Müller, R. 15N-NMR Spectroscopy— New Methods and Applications [New Analytical Methods (28)], Angew. Chem. Int. Ed. 1986, 25, 383–413. (23) Mason, J. Nitrogen nuclear magnetic resonance spectroscopy in inorganic, organometallic, and bioinorganic chemistry, Chem. Rev. 1981, 81, 205–227. (24) Witanowski, M.; Biedrzycka, Z.; Sicinska, W.; Grabowski, Z. A Study of Solvent Polarity and Hydrogen Bonding Effects on the Nitrogen NMR Shielding of Isomeric Tetrazoles and Ab Initio Calculation of the Nitrogen Shielding of Azole Systems, J. Magn. Reson. 1998, 131, 54–60. (25) Solum, M. S.; Altmann, K. L.; Strohmeier, M.; Berges, D. A.; Zhang, Y.; Facelli, J. C.; Pugmire, R. J.; Grant, D. M. 15N Chemical Shift Principal Values in Nitrogen Heterocycles, J. Am. Chem. Soc. 1997, 119, 9804–9809. (26) Cui, Q.; Karplus, M. Molecular Properties from Combined QM/MM Methods. 2. Chemical Shifts in Large Molecules, J. Phys. Chem. B 2000, 104, 3721–3743. (27) Steinmann, C.; J. Olsen, M. H.; Kongsted, J. Nuclear Magnetic Shielding Constants from Quantum Mechanical/Molecular Mechanical Calculations Using Polarizable Embedding: Role of the Embedding Potential, J. Chem. Theory Comput. 2014, 10, 981–988. (28) Flaig, D.; Beerand, M.; Ochsenfeld, C. Convergence of Electronic Structure with the Size of the QM Region: Example of QM/MM NMR Shieldings, J. Chem. Theory Comput. 2012, 8, 2260– 2271. (29) 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. (30) 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. (31) 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. (32) Neidigh, J. W.; Fesinmeyer, R. M.; Andersen, N. H. Designing a 20-residue protein, Nat. Struct. Biol. 2002, 9, 425–430. (33) Mennucci, B.; Martínez, J. M.; Tomasi, J. Solvent Effects on Nuclear Shieldings: Continuum or Discrete Solvation Models to Treat Hydrogen Bond and Polarity Effects? J. Phys. Chem. A 2001, 105, 7287−7296.

Page 8 of 10

(34) Sebaoun, L.; Maurizot, V.; Cranier, T.; Kauffmann, B.; Huc, I. Aromatic Oligoamide β-Sheet Foldamers, J. Am. Chem. Soc. 2014, 1136, 2168–2174. (35) Mathias, J. P.; Simanek, E. E.; Zerkowski, J. A.; Seto, C. T.; Whitesides, G. M. Structural Preferences of Hydrogen-Bonded Networks in Organic Solution−The Cyclic CA3·M3 ″Rosette″, J. Am. Chem. Soc. 1994, 116, 4316–4325. (36) 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. (37) 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. (38) 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. (39) 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. (40) Fang, T.; Li, Y.; Li, S. Generalized energy-based fragmentation approach for modeling condensed phase systems, WIREs Comput. Mol. Sci. 2017, 7, e1297–e1309. (41) Wang, F.; Zhao, D.; Dong, H.; Jiang, L.; Liu, Y.; Li, S. Terahertz spectra of DNA nucleobase crystals: A joint experimental and computational study, Spectrochim. Acta A 2017, 179, 255–260. (42) Hua, W.; Fang, T.; Li, W.; Yu, J.; 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. (43) Li, W.; Chen, C.; Zhao, D.; Li, S. LSQC: Low Scaling Quantum Chemistry Program, Int. J. Quantum Chem. 2015, 115, 641–646. (44) 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski J.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2013. (45) London, F. The quantic theory of inter-atomic currents in aromatic combinations, J. Phys. Radium 1937, 8, 397–409. (46) McWeeny, R. Perturbation Theory for Fock-Dirac Density Matrix, Phys. Rev. 1962, 126, 1028–1034. (47) Ditchfield, R. Self-consistent perturbation theory of diamagnetism. 1. Gauge-invariant LCAO method for N.M.R. chemical shifts, Mol. Phys. 1974, 27, 789–807. (48) Wolinski, K.; Hilton, J. F.; Pulay, P. A Comparison of Models for Calculating Nuclear Magnetic Resonance Shielding Tensors, J. Am. Chem. Soc. 1990, 112, 8251–8260. (49) Cheeseman, J. R.; Trucks, G. W.; Keith, T. A.; Frisch, M. J. A comparison of models for calculating nuclear magnetic resonance shielding tensors, J. Chem. Phys. 1996, 104, 5497–5509. (50) Wilson, P. J.; Bradley, T. J.; Tozer, D. J. Hybrid exchangecorrelation functional determined from thermochemical data and ab initio potentials, J. Chem. Phys. 2001, 115, 9233–9242. (51) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent Molecular Orbital Methods. 25. Supplementary Functions for Gaussian Basis Sets, J. Chem. Phys. 1984, 80, 3265−3269. (52) Jensen, F. Segmented Contracted Basis Sets Optimized for Nuclear Shielding, J. Chem. Theory Comput. 2015, 11, 132−138. (53) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; D. A. Case, Development and testing of a general amber force field, J. Comput. Chem. 2004, 25, 1157–1174.

ACS Paragon Plus Environment

Page 9 of 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(54) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints For Determining Atom-Centered Charges: The RESP Model, J. Phys. Chem. 1993, 97, 10269–10280. (55) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A. Application of RESP charges to calculate conformational energies, hydrogen bond energies and free energies of solvation, J. Am. Chem. Soc. 1993, 115, 9620–9631. (56) Adelman, S. A.; Doll, J. D. Generalized Langevin equation approach for atom/solid-surface scattering: General formulation for classical scattering off harmonic solids, J. Chem. Phys. 1976, 64, 2375–2388. (57) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms, J. Chem. Phys. 1994, 101, 4177– 4189. (58) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: The Langevin piston method, J. Chem. Phys. 1995, 103, 4613–4621. (59) Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An overview of the Amber biomolecular simulation package, WIREs Comput. Mol. Sci. 2013, 3, 198–210. (60) Korth, M. Third-Generation Hydrogen-Bonding Corrections for Semiempirical QM Methods and Force Fields, J. Chem. Theory Comput. 2010, 6, 3808–3816. (61) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 1988, 38, 3098– 3100. (62) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 1988, 37, 785–789. (63) Grimme, S. Semiempirical GGA-type density functional constructed with a long range dispersion correction, J. Comput. Chem. 2006, 27, 1787–1799. (64) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach, Comput. Phys. Commun. 2005, 167, 103–128.

(65) Lippert, G,; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme, Mol. Phys. 1997, 92, 477–487. (66) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J. Chem. Phys. 2007, 127, 114105–114113. (67) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials, Phys. Rev. B 1996, 54, 1703–1710. (68) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn, Phys. Rev. B 1998, 58, 3641–3662. (69) Roothaan, C. C. J. New Developments in Molecular Orbital Theory, Rev. Mod. Phys. 1951, 23, 69–89. (70) Yanai, T.; Tew, D.; Handy, N. A new hybrid exchangecorrelation functional using the Coulomb-attenuating method (CAMB3LYP), Chem. Phys. Lett. 2004, 393, 51–57. (71) Zhao, Y.; Truhlar, D. G. 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. (72) Hartman, J. D.; Beran, G. J. O. Fragment-Based Electronic Structure Approach for Computing Nuclear Magnetic Resonance Chemical Shifts in Molecular Crystals, J. Chem. Theory Comput. 2014, 10, 4862−4872. (73) Hartman, J. D.; Kudla, R. A.; Day, G. M.; Mueller, L. J.; Beran, G. J. O. Benchmark fragment-based 1H, 13C, 15N and 17O chemical shift predictions in molecular crystals, Phys. Chem. Chem. Phys. 2016, 18, 21686−21709. (74) Witanowski, M.; Stefaniak, L.; Webb, G. In Annual Reports on NMR Spectroscopy; Webb, G. A., Ed.; Academic Press: London, 1993; Vol. 25. (75) Bouř, P.; Dračínský, M. Computational Analysis of Solvent Effects in NMR Spectroscopy, J. Chem. Theory Comput. 2010, 6, 288–299.

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 10

Table of Contents

ACS Paragon Plus Environment

10