Vibrational Spectra of Molecular Crystals with the Generalized Energy

Article pubs.acs.org/JPCA

Vibrational Spectra of Molecular Crystals with the Generalized Energy-Based Fragmentation Approach Tao Fang, Junteng Jia, and Shuhua Li* School of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, Nanjing 210093, People’s Republic of China ABSTRACT: The generalized energy-based fragmentation (GEBF) approach for molecular crystals with periodic boundary condition (PBC) (denoted as PBC-GEBF) is extended to allow vibrational spectra of molecular crystals to be easily computed at various theory levels. Within the PBC-GEBF approach, the vibrational frequencies of a molecular crystal can be directly evaluated from molecular quantum chemistry calculations on a series of nonperiodic molecular systems. With this approach, the vibrational spectra of molecular crystals can be calculated with much reduced computational costs at various theory levels, as compared to those required by the methods based on periodic electronic structure theory. By testing the performance of the PBCGEBF method for two molecular crystals (CO2 and imidazole), we demonstrate that the PBC-GEBF approach can reproduce the results of the methods based on periodic electronic structure theory in predicting vibrational spectra of molecular crystals. We apply the PBC-GEBF method at second-order Møller−Plesset perturbation theory (PBC-GEBF-MP2 in short) to investigate the vibrational spectra of the urea and ammonia borane crystals. Our results show that the PBC-GEBF-MP2 method can provide quite accurate descriptions for the observed vibrational spectra of the two systems under study.

I. INTRODUCTION Vibrational infrared (IR) and Raman spectra are important experimental methods to characterize structural and electronic features of molecules and solids. The interpretation of these spectra requires accurate calculations from theoretical approaches. The calculation of vibrational spectra of molecular systems can be routinely performed with existing quantum chemistry packages. However, for periodic systems, such as molecular crystals, accurate prediction of vibrational spectra remains a challenging task. The cluster model provides a simple way of evaluating the structure and vibrational spectra of molecules in the crystal environment.1,2 In this model, a supermolecule composed of a monomer and a dozen of its neighboring monomers is constructed. With existing molecular quantum chemistry packages, one can obtain fairly accurate predictions on the molecular structure and intramolecular vibration modes for some molecular crystals. However, the cluster model is rarely used to provide accurate predictions on the lattice vibration modes and the phonon-dispersion curves of molecular crystals, because the periodic condition for lattice is not taken into account. Moreover, because the cluster model usually consists of more than a dozen monomers, it is still computationally prohibitive to investigate the vibrational spectra with electron correlation methods. At present, the most popular way of investigating vibrational spectra of solids is to employ periodic electronic structure methods and quantum chemistry packages based on periodic boundary conditions. For instance, periodic density functional theory (DFT) and © 2016 American Chemical Society

Hartree−Fock (HF) methods are widely employed in the vibrational calculations of molecular crystals. However, periodic post-Hartree−Fock methods such as second-order Møller− Plesset perturbation theory (MP2) have not been implemented for vibrational spectra calculations of molecular crystals. Different from standard periodic electronic structure methods, several fragment-based approaches for predicting the structures and properties of molecular crystals have been proposed.3−17 In most fragment-based methods, the energy and properties of a molecular crystal are approximately evaluated with molecular quantum chemistry calculations on small molecular clusters (periodic boundary conditions are explicitly employed in constructing molecular clusters). Some of these methods have also been further developed to allow vibrational spectra of molecular crystals to be calculated.5,8,9,12 For instance, the binary interaction method (BIM) proposed by Hirata’s group6 has been demonstrated to give quite satisfactory results for several molecular crystals. In the BIM approach, the unit cell energy is approximated by the linear combination of energies calculated from a series of monomers and molecule dimers, which are embedded in a self-consistently determined electrostatic field (simplified by a huge number of point charges), and a correction term responsible for the long-range electrostatic interaction. Within the BIM method, one can Received: November 7, 2015 Revised: April 13, 2016 Published: April 13, 2016 2700

DOI: 10.1021/acs.jpca.5b10927 J. Phys. Chem. A 2016, 120, 2700−2711

Article

The Journal of Physical Chemistry A easily calculate the ground-state unit cell energy, and the firstand second-order derivatives of unit cell energy with respect to atomic coordinates. This approach has been utilized to predict the crystal structure, lattice cohesive energy, vibrational spectra, and phonon dispersion curves for several systems. Especially, vibrational spectra of several molecular crystals including hydrogen fluoride, Ih phase of ice, and carbon dioxide calculated with the BIM method have been reported.8,9 It should be mentioned that the BIM method may be less accurate for some molecular crystals, in which the trimer and tetramer interactions are important. It should be mentioned that a different way of taking long-range electrostatic interactions of the crystal environment into account has also been proposed for fragment-based methods of periodic systems.18 Recently, we have developed the generalized energy-based fragmentation (GEBF) approach19−24 for molecular crystals with the periodic boundary condition (PBC).25 In the so-called PBC-GEBF method, the energy of a unit cell in a molecular crystal can be expressed as a linear combination of energies of a series of embedded subsystems, which are composed of several monomers. Each subsystem is embedded in the electric field produced by background point charges generated by all other atoms in neighboring cells and compensation charges, which are employed to mimic the rest of background point charges in distant cells. The use of compensation charges allows the longrange polarization effect and electrostatic interactions of the crystal environment to be treated with reasonable accuracy. Illustrative applications demonstrated that the PBC-GEBF method can provide accurate lattice energies and structures for a wide range of molecular crystals. In this article, the calculation of vibrational spectra of molecular crystals within the PBC-GEBF approach is implemented. Within the PBC-GEBF approach, the secondorder derivatives of the unit cell energy of a molecular crystal can be estimated by combining the corresponding derivatives of a series of subsystems. As in the conventional periodic electronic structure methods, we then build the dynamical matrix of a molecular crystal and diagonalize it to yield vibrational frequencies and corresponding normal modes of this crystal. On the other hand, the derivatives of the dipole moment and polarizability of this crystal with respect to the atomic coordinates can also be readily obtained with the corresponding quantities of subsystems and the normal mode components of the crystal. These derivative matrixes could be employed to evaluate the IR and Raman vibrational intensities. Thus, the PBC-GEBF method allows the vibrational frequencies and intensities of molecular crystals to be calculated with molecular quantum chemistry methods. Especially, the calculation of vibrational spectra for various molecular crystals composed of small molecules at the MP2 level is now feasible on ordinary workstations. This article is organized as follows: The methodology and computational details are introduced in section II. In section III, we first apply the PBC-GEBF approach to investigate vibrational spectra in two molecular crystals, and compare the results from PBC-GEBF calculations and those from calculations based on periodic electronic structure theory. The PBC-GEBF method at the MP2 level then is applied to predict the vibrational spectra of the urea and ammonia borane crystals. The predicted vibrational spectra are in good agreement with the experimental spectra. Finally, a brief conclusion is given in section IV.

II. METHODOLOGY II.A. Unit Cell Energy and Force in the PBC-GEBF Method. In the PBC-GEBF method, the unit cell energy of a molecular crystal is evaluated as25 M

Eunit‐cell =

M

∑ CmEm̃ − (∑ Cm) ∑ ∑ m

m

A ∈ K (B > A) ∈ K

+ E Ewald‐sum

Q AQ B RAB (1)

in which Em̃ stands for the ground-state energy of an “electrostatically embedded” subsystem, which is composed of a few monomers and embedded in a point charge grid in a super cell K. Here, the super cell K is an array of the crystallographic unit cell within a cutoff distance, and Cm represents the coefficient of this subsystem. M represents the total number of subsystems. QA denotes the natural population atomic (NPA) charge at atom A, and RAB stands for the interatomic distance. EEwald‑sum is the classical electrostatic interaction within the central cell and between the central cell and all image cells, which can be easily computed with the Ewald summation method:26 E Ewald‐sum =

1 2

∑ ∑ Q AQ B n

A ,B

erfc(α rAB , n) rAB , n

⎤2 1 4π −k2 /4α 2⎡ ⎢∑ Q A exp(i k·rA)⎥ + ∑ e ⎢⎣ A ⎥⎦ 2V k ≠ 0 k 2 α − ∑ Q2 π A A

(2)

in which rAB,n stands for the distance between atom A in the original cell and atom B in cell n = (n1,n2,n3), and k represents a reciprocal space vector. α is a positive parameter that balances the computational cost for evaluating the real space term and the reciprocal space term, and in our calculation we always choose α = 0.37. V stands for the volume of the unit cell. In Appendix A, we will provide some discussions on how eq 1 is derived. The size of the super cell is controlled by a cutoff distance Rcutoff. For a three-dimensional periodic system, the super cell is usually selected to contain all unit cells, whose distance from the central unit cell in at least one dimension is less than the cutoff distance. The number of unit cells can be determined in the following way. Assume that the lattice parameter of a molecular crystal is a, b, and c, the number of unit cells on both sides of the central unit in each dimension should be Na = 1 + int(Rcutoff/a), Nb = 1 + int(Rcutoff/b), and Nc = 1 + int(Rcutoff/c), respectively. Therefore, the total number of unit cells in the super cell should be (2Na + 1)(2Nb + 1)(2Nc + 1). Our calculations show that the default value for the cutoff distance may be set to be 18 Å. In Figure 1, we provide an illustrative picture of a subsystem, which is embedded in the electrostatic field (a finite array of point charges) of a super cell. All atoms in the subsystem are displayed in ball-and-stick format (in the central unit cell), while the point charges are shown in wire frame format. It is worth mentioning that the compensating point charges at the boundaries of the super cell (discussed later in this section) are neglected in Figure 1 for clarity. For molecular crystals, each monomer is naturally defined as a fragment, and a subsystem is defined to be composed of a few fragments. It is well-known that the unit cell energy of a 2701


Article

The Journal of Physical Chemistry A

generating derivative subsystems with n fragments, we then move to construct derivative subsystems with (n − 1) fragments. Repeating this process will automatically generate all derivative subsystems. (3) The periodic boundary condition is employed in generating all primitive and derivative subsystems. It is also worth mentioning the utilization of the compensation field scheme. In deriving the unit cell energy for a molecular crystal, we have assumed that each subsystem is embedded in the presence of a finite array of point charges in a supercell (rather than an infinite array of point charges). This approximation works well only if the dipole moment of each unit cell vanishes. To achieve this goal, we add some compensating point charges at the boundaries of each unit cell to ensure that the overall dipole moment of each unit cell is zero. Because compensating charges from neighboring unit cells naturally cancel each other inside the supercell, only those on the boundaries of the super cell are necessary. With this compensation field, the PBC-GEBF method is demonstrated to be almost independent of the choice of unit cell (or nearly translationally invariant). For instance, the phase-II ice crystal is investigated with the PBC-GEBF method at the B3LYP/631G(d,p) level. The unit cell can be chosen in different ways by translating the atomic coordinates by a different distance. It is demonstrated that the unit cell energy difference is no more than 1.0 millihartree for any two different unit cell choices. Detailed descriptions on the implementation of the field compensation method and the calculation of background point charges can be found in our previous work.25 The gradient of the unit cell energy with respect to atomic coordinates as well as the lattice translation vectors can be obtained as

Figure 1. Illustrative picture of a super cell constructed in a PBCGEBF calculation for molecular crystal of carbon dioxide. The subsystem contains 4 CO2 monomers (shown in ball-and-stick format) in the central unit cell, which are embedded in the array of point charges (shown in wire frame format) in other unit cells of the super cell.

periodic molecular crystal can be expanded as the summation of intrafragment and n-fragment energy terms, and the total energy of a subsystem (or a molecular cluster) has a similar many-body expansion. The basic idea of the PBC-GEBF method is to evaluate the unit cell energy of a periodic molecular crystal as the combination of energies of a set of nonperiodic subsystems (or molecular clusters), which can be calculated with existing molecular electronic structure methods. In Appendix B, we have described the many-body expansion of the total energy of a subsystem or the unit cell energy of a molecular crystal, and some ideas on how to construct subsystems. The detailed procedure on how to construct subsystems and derive their coefficients was described previously.25 Here, we only give a brief introduction on the procedure: (1) For each fragment, we add its neighboring fragments within a given distance threshold (ζ) to build a subsystem. To control the size of a subsystem, the parameter (γmax), which denotes the maximum number of fragments, is introduced so that only the closest γmax fragments are included in the subsystem. In addition, for a pair of fragments within a given distance ζ, we also build a subsystem centered around the middle point of this fragment pair in a similar manner. These subsystems are called primitive subsystems, and their coefficients are all set to be +1. (2) Derivative subsystems are constructed sequentially. At first, we find out those n-fragment interactions that occur in all primitive subsystems. The initial value of n is set to (γmax − 1). If an n-fragment interaction occurs more than once from existing subsystems, a derivative subsystem containing the fragments involved has to be constructed. The principal rule for determining the coefficient of a derivative subsystem is to guarantee that the summation of coefficients of any specific nfragment interaction term from all subsystems equals +1. After

∂Eunit‐cell = ∂rA ,0

⎡

i ∈ Box

⎛ ⎞ + ⎜⎜∑ Cm⎟⎟ ⎝m ⎠ ∂Eunit‐cell = ∂L

⎛

∂Em̃ − Fm , a , iQ a − ⎝ ∂rA , i

∑ ⎢⎢∑ Cm⎜⎜ ⎣

m

⎞ fai , bj⎟⎟ j ∈ Box, b ⎠

∑

⎤ ∂E ∑ fai ,bj⎥⎥ + Ewald‐sum ∂rA ,0 j ∈ Box, b ⎦

∑ i ∈ Box

⎡ ⎛ ̃ ∂E iL⎢∑ Cm ∑ ⎜⎜ m − Fm , a , iQ a ⎢ A ⎝ ∂rA , i ⎣m

−

⎞ ⎛ ⎞ fai , bj⎟⎟ + ⎜⎜∑ Cm⎟⎟ ∑ ⎠ a j ∈ Box, b ⎠ ⎝m

+

∂E Ewald‐sum ∂L

∑

(3)

⎤ fai , bj⎥ ⎥ j ∈ Box, b ⎦

∑

(4)

In the above equations, the gradients of Em̃ with respect to atom coordinates can be obtained from a conventional force calculation on the embedded subsystem m. The capital letter A denotes a real atom in a given subsystem, while a and b denote the point charge centers located at atoms A and B, and i, j denotes the cell index. Qa denotes the NPA charge on atom A. Fm,a,i stands for the electric field generated by the subsystem m on the center a in the ith cell, and fai,bj stands for the electrostatic force between the charge on a in the ith cell and the charge on b in the jth cell. It should be mentioned that the derivatives of the NPA charge with respect to atomic coordinates (or the lattice translation vectors) in the above 2702


Article

The Journal of Physical Chemistry A equations are ignored. Our previous work on nonperiodic systems demonstrated that the neglect of these terms to the gradient is a reasonable approximation.24 All quantities mentioned above can be obtained by calculations on subsystems with existing quantum chemistry packages. The last term of the right-hand side in eq 3 or eq 4 stands for the gradient of the classical electrostatic interaction of the infinite crystal with respect to atomic coordinates or lattice translation vectors, which can be calculated analytically. It should be pointed out that in the GEBF method for molecular systems, the gradient on a given atom can be further approximated as the combination of the corresponding gradients on this atom in those subsystems explicitly containing this atom.24 In a similar way, in the GEBF approach for periodic systems, the gradient of the unit cell energy with respect to the atomic coordinates, as formulated in eq 3, can also be approximately evaluated as ∂Eunit‐cell = ∂rA ,0

∑ ∑ Cm i ∈ Box m

∂Em̃ ∂rA , i

1 mA mB

D(rA, rB, k) =

∑ H(rA ,0, rB,l) e−ikR(l) l

(6)

where k defines a given point in the Brillouin zone, and H(rA,0, rB,l) stands for the second-order derivative of the total energy per unit cell with respect to atom A in the 0th cell and atom B in the lth cell at the equilibrium geometry. In addition, mA and mB denote the mass of atom A and atom B, respectively. The observed vibrational spectrum is closely related to the dynamical matrix at Γ point, where k = (0,0,0). Here, we focus on the Γ point only. For our present purpose, the dynamical matrix is reduced to 1 mA mB

D(rA, rB) =

∑ H(rA ,0, rB,l) l

(7)

The second-order derivative H(rA,0, rB,l) can be evaluated in a similar way as the energy gradient with respect to atomic coordinates, as shown in eq 5. With a similar treatment, we can evaluate the second-order derivative matrix H(rA,0, rB,l) with the corresponding second-order derivatives obtained from calculations on “embedded” subsystems, as shown below:

(5)

Here, the summation over m is limited to those subsystems, in which atom A in cell i is included as a real atom. Thus, the gradient on a given atom can be directly evaluated as a linear combination of the corresponding gradients on this atom (in both central unit cell and other unit cells in the super cell) in those subsystems explicitly containing this atom. Our numerical tests show that the gradients estimated with eq 5 are excellent approximations to those obtained with eq 3. For instance, for 10 randomly selected structures of the imidazole crystal at the HF/6-311G(d,p) level, the root-mean-squared deviations (RMSDs) between the gradients with respect to atomic coordinates obtained with eq 3 and those obtained with eq 5 for all atoms are about 0.0003 au/bohr. Thus, one may employ the simple expression eq 5, instead of eq 3, to optimize the structure of a molecular crystal. Furthermore, on the basis of eq 5, the second-order derivatives of the unit cell energy with respect to the atomic coordinates can be easily evaluated with existing molecular quantum chemistry programs, as discussed in the following subsection. Different from most other fragment-based methods, the PBC-GEBF approach relies on the compensation field to take the effect of the infinite point charge arrays into account. Therefore, we do not need to calculate subsystems in the electric field composed of a huge number of point charges located in a giant super cell (as in the BIM method). It would be interesting to compare the accuracy and computational cost of different fragment-based methods for molecular crystals in the future. Although the PBC-GEBF approach is applied to treat molecular crystals consisting of small molecules in the present study, the PBC-GEBF method has the potential to treat molecular crystals consisting of larger molecules. For molecular crystals composed of large molecules, constructed subsystems described in the preceding subsection may be too large for ab initio (especially post-HF) calculations. In such cases, we may further apply the GEBF approach to simplify the calculation of these subsystems. We will report the application of the PBCGEBF approach to molecular crystals with large monomers in the near future. II.B. Vibrational Frequencies and Intensities of Molecular Crystals. The dynamical matrix of a periodic system can be expressed as

H(rA ,0, rB , l) =

∂ 2Em̃ ∂rA ,0rB , l

∑ Cm m

(8)

in which A and B are real atoms (not point charges) in a given subsystem and l denotes the cell index. The summation over m is limited to those subsystems in which atom A in cell 0 and atom B in cell l are included as real atoms. The second-order derivatives occurring in the right-hand side of eq 8 can be directly obtained from corresponding calculations for each subsystem. Once the dynamical matrix of a periodic system is available, one can solve the following equation to obtain the vibrational frequencies as well as the corresponding normal modes: |D(rA, rB) − ω 2 I| = 0

(9)

in which ω stands for the vibrational frequency and I denotes a unit matrix. It should be mentioned that the above procedure can also be used to compute phonon dispersion curves, that is, vibrational frequencies beyond the Γ point. In this case, an even larger supercell may be needed to obtain accurate results. The calculation of the IR and Raman intensity of each normal mode in a molecular crystal is similar to that described previously by our group for molecular systems.22 As in molecular systems, the derivatives of the dipole moment and the polarizability tensor with respect to the atomic coordinates for a molecular crystal within the PBC-GEBF approach can be evaluated as ∂Ω ≈ ∂rA

∑ Cm m

∂Ωm ∂rA

(A ∈ subsystem m , Ω = μ , α) (10)

in which Ω stands for the dipole moment (denoted as μ) or the polarizability tensor (denoted as α). Thus, the derivative of the dipole moment (or polarizability) with respect to an atomic coordinate in a molecular crystal can be expressed as a linear combination of the corresponding quantities of all subsystems, which can be obtained with conventional molecular electronic structure calculations. Furthermore, the derivative of the dipole moment (or polarizability) with respect to normal modes (Qk) can be derived from the following equation: 2703


Article


Table 1. Comparison of Optimized Lattice Parameters for Two Molecular Crystals from PBC-GEBF Calculations with Those Obtained from the Conventional Calculationsa system

method

a (Å)

b (Å)

c (Å)

α (deg)

β (deg)

γ (deg)

cell volume

imidazole

PBC-GEBF conventional PBC-GEBF conventional

8.648 8.658 5.383 5.384

6.206 6.200 5.383 5.384

10.144 10.137 5.383 5.384

90.0 90.0 90.0 90.0

107.5 107.5 90.0 90.0

90.0 90.0 90.0 90.0

519.19 518.97 155.97 156.07

CO2 a

For the CO2 crystal, the PBE0 functional with Grimme’s D2 dispersion scheme (PBE0-D2 in short) is employed, and for the imidazole crystal, the Hartree−Fock method is employed. The 6-311G(d,p) basis set is used in all calculations.

Figure 2. IR and Raman spectra of the crystalline imidazole at the HF/6-311G(d,p) level from (a) the conventional calculation and (b) the PBCGEBF calculation. Comparison of conventional and PBC-GEBF results using the Lorentz broadening with HWHM = 20 cm−1 is shown in (c).

∂Ω = ∂Q k

3N

∑ a=1

∂Ω ∂rA ∂rA ∂Q k

II.C. Computational Details. To calculate the vibrational spectra of a molecular crystal, we should first optimize its crystal structure. The vibrational frequency calculation then is carried out at the optimized structure with the same method as in the geometry optimization. Finally, the calculated vibrational frequencies and corresponding intensities are transformed into IR and Raman spectra with a Lorentz broadening program with a half-width at half-maximum (HWHM) being 20 cm−1. In our calculations, the PBC-GEBF program is interfaced with the Gaussian 09 program,28 which is employed to perform calculations for each embedded subsystem. In general, our previous calculations demonstrated that two parameters (ζ = 4.0 Å and γmax = 4) are suitable for PBC-GEBF calculations of molecular crystals, and thus will be employed in all PBC-GEBF calculations. To validate the performance of the PBC-GEBF approach, two typical molecular crystals, CO2 and imidazole, are selected for illustrative purpose. The results from PBC-GEBF calculations will be compared to those obtained from conventional periodic calculations, which are performed with the CRYSTAL package.29−31 The results shown in the next section demonstrate that the PBC-GEBF method is capable of reproducing the results from much more expensive conventional periodic calculations.

(k = 1, 2, ..., 3N ) (11)

With these derivatives, the IR intensities Ik and Raman Intensities Rk of a molecular crystal can be evaluated with the equations:27 Ik =

∂μ ⃗ ∂rA

2

(12)

R k = (45αk̅ 2 + 7γk2)

(13)

where αk̅ =

γk2

1 3

∑ a=x ,y,z

∂αaa ∂Q k

⎡ ⎤ ⎛ ∂α ⎞2 1⎢ 3 ∑ ⎜⎜ ab ⎟⎟ − 9αk̅ 2 ⎥ = ⎥ 2 ⎢⎣ a , b = x , y , z ⎝ ∂Q k ⎠ ⎦

(14)

(15)

In the above equations, μ⃗ stands for the dipole moment, and αaa and αab denote the components of the polarizability tensor. 2704


Article


Figure 3. IR and Raman spectra of the crystalline CO2 at the PBE0-D2/6-311G(d,p) level from (a) the conventional calculation and (b) the PBCGEBF calculation. Comparison of conventional and PBC-GEBF results using the Lorentz broadening with HWHM = 20 cm−1 is shown in (c).

PBC-GEBF method provides a reliable tool for predicting the crystal structures of molecular crystals. The IR and Raman spectra for both systems are displayed in Figure 2 and Figure 3, respectively. The computed peaks from the conventional calculations are displayed in Figure 2a and Figure 3a, while the computed peaks from the PBC-GEBF calculations are displayed in Figure 2b and Figure 3b, respectively. The spectra generated with a Lorentz broadening program are displayed in Figure 2c and Figure 3c for comparison. It can be observed from two figures that the shapes of the spectrum from PBC-GEBF calculations are in good agreement with those from conventional calculations. The main peaks from both approaches generally differ from each other by a few wavenumbers, indicating that the PBC-GEBF approach provides satisfactory predictions on vibrational frequencies. The overall intensity of the peak after Lorentz transformation is still in satisfactory agreement with that from conventional calculations, as illustrated in Figure 2 and Figure 3. From this direct comparison with the conventional calculation, one can conclude that the PBC-GEBF method is capable of providing satisfactory predictions on the vibrational IR and Raman spectra. It is worthwhile to analyze the effect of the embedding charges on the calculated vibrational frequencies. We take the crystalline imidazole system as an example to demonstrate this effect. We have performed PBC-GEBF calculations with and without embedding charges to optimize the crystal structure and then computed the vibrational frequencies. The calculated volume of the unit cell decreases from 518 Å3 (with embedding charges) to 480 Å3 (without embedding charges), indicating that the embedding charges are of great significance in the structure of the imidazole crystal. Furthermore, vibrational spectra obtained with and without the embedding charges are displayed in Figure 4. One can observe that the IR and Raman peaks obtained without the embedding charges show significant

For systems under study, space group symmetry is exploited in all calculations. The symmetry of a molecular crystal is maintained throughout the geometry optimization process. Moreover, space group symmetry can also be utilized to identify structurally equivalent subsystems so that the computational cost can be considerably reduced. Calculations are performed only for unique subsystems. For those symmetrically equivalent subsystems, the quantities such as energies, gradient, and Hessian matrix are obtained from transformation of the corresponding quantities of those subsystems to be calculated. The detailed procedure on how to exploit space group symmetry will be described elsewhere.

III. RESULTS AND DISCUSSION III.A. Validation of the PBC-GEBF Method for Computing Vibrational Spectra of Molecular Crystals. The PBC-GEBF approach is first applied to two typical molecular crystals to validate its performance in describing vibrational spectra of molecular crystals. Two theoretical methods are employed in PBC-GEBF calculations. For the van der Waals bonded crystals CO2 (Pa3 phase), calculations are performed with the PBE0 functional32 with Grimme’s D2 dispersion scheme,33 which are implemented in both Gaussian 09 and CRYSTAL packages. For the hydrogen-bonded imidazole crystal (P21/c phase), the calculations are performed at the Hartree−Fock level. The 6-311G(d,p) basis set is employed in all calculations on these two systems. The optimized lattice parameters are displayed in Table 1. It can be observed from the table that the lattice parameters obtained with PBC-GEBF are in excellent agreement with those estimated by conventional periodic calculations. The average deviation for unit cell lengths is less than 1%, and the unit cell angles are almost the same from these two approaches. The deviations of the unit cell volume for CO2 and imidazole crystals are both less than 1%. These results indicate that the 2705


Article


The urea crystal has been investigated extensively in the past decade.1,35−37 Some experimental investigations on the IR spectra of the urea crystal as well as urea−water solution35,36 and some theoretical studies based on cluster models have been reported.1,37 For instance, in ref 1, the authors employed a cluster model with 15 urea molecules and 5312 background charges to reproduce the crystal environment. Because of the size of this cluster model, theoretical calculations were performed only at the RHF/6-31++G** level. Here, we will investigate vibrational spectra of the urea crystal (in its P4̅21m crystal phase) with the PBC-GEBF-MP2 approach. The optimized lattice parameters and the corresponding experimental parameters are listed in Table 2 for comparison, and the computed and observed IR spectra are displayed in Figure 5. At first, one can see that the PBC-GEBF-MP2 method Figure 4. Comparison of IR and Raman spectra of the crystalline imidazole at the HF/6-311G(d,p) level obtained from the PBC-GEBF method with and without embedding charges.

deviations from those obtained with the embedding charges. The main peak in high frequency zone is red-shifted by 110− 120 cm−1 in both spectra. For IR spectra, the main peak is shifted from 3613 to 3503 cm−1, while for Raman spectra, the main peak is shifted from 3605 to 3486 cm−1. The analysis of the corresponding vibrational modes shows that these frequencies correspond to the N−H bond stretching. Examining the interatomic distances between neighboring imidazole molecules in the optimized crystal structures, we find that the distance of the N−H···N hydrogen bond is 3.017 Å (the N−H distance of 1.008 Å) from the calculation with the embedding charges, and 2.976 Å (the N−H distance of 1.012 Å) from the calculation without the embedding charges. Hence, the absence of embedding charges leads to stronger N−H···N hydrogen bonds and somewhat weaker N−H bonds, which are responsible for the red shift of the N−H bond stretching. Therefore, the inclusion of the embedding charges is critical for accurate predictions of the vibrational frequencies in hydrogenbonded molecular crystals. III.B. IR Spectra of the Urea Crystal. In the following, the PBC-GEBF approach at the MP2 level (PBC-GEBF-MP2 in short) will be employed to investigate the vibrational spectra of the urea and ammonia borane crystals. The 6-311++G(d,p) basis set is employed in all calculations. Because MP2 calculations tend to overestimate the vibrational frequencies with respect to the experimental data, scaling factors (which depend on the basis set used) are used to scale the frequencies obtained. Previous work on molecular systems suggested different scaling factors for different combinations of methods and basis sets.34 For the basis set used in this work, the factor is chosen as 0.9523 for MP2.

Figure 5. Comparison of IR spectra of the crystalline urea obtained from PBC-GEBF-MP2 calculations with the experimental spectra. Reused with permission from ref 36, J. Mol. Struct. 2002, 615, 177−189. Copyright 2002 Elsevier Science B.V.

gives very accurate predictions on the crystal structure. The deviations from the experimental data are no more than 3% for all three lattice parameters. From the comparison of calculated and experimental IR spectra in Figure 5, one can find that the calculated IR spectra are in good agreement with the observed IR spectra. The predicted frequencies agree quite well with the experiment values, although the relative height of some peaks is somewhat different from that in the experimental spectra. A detailed comparison between calculated and experimental IR spectra can be employed to assign the experimental peaks with vibrational modes. The NH2 stretching frequencies predicted by PBC-GEBF-MP2 are 3497 and 3322 cm−1, in good accordance with the experiment peaks (3437 cm−1, 3343 cm−1). The NH2 bending frequency predicted by PBC-GEBFMP2 is 1682 cm−1, which should correspond to the experimental peak 1668 cm−1. The calculated peak at 1447 cm−1 can also be clearly assigned to the CN stretching mode observed at 1466 cm−1. Calculated peaks with relatively smaller intensity at 1115 and 1045 cm−1 represent the mode of ρ(NH2), corresponding to the experiment peaks at 1156 and 1057 cm−1, respectively. One relatively broad peak centered at

Table 2. Comparison of Optimized Lattice Parameters Obtained from PBC-GEBF-MP2 Calculations for the Urea and BH3NH3 Crystals with Those from the Experimental Crystal Structuresa system urea BH3NH3

a

method

a (Å)

b (Å)

c (Å)

α (deg)

β (deg)

γ (deg)

cell volume

PBC-GEBF-MP2 X-ray PBC-GEBF-MP2 X-ray

5.471 5.565 5.366 5.395

5.471 5.565 5.152 4.887

4.823 4.684 4.955 4.986

90.0 90.0 90.0 90.0

90.0 90.0 90.0 90.0

90.0 90.0 90.0 90.0

144.3 145.1 137.0 131.4

The 6-311++G(d,p) basis is employed in all calculations. 2706


Article

The Journal of Physical Chemistry A 1598 cm−1 from the PBC-GEBF-MP2 calculation corresponds to the experimental peaks at 1624 cm−1 (another NH2 bending mode) and 1599 cm−1 (the CO stretching mode). Furthermore, it is also worth analyzing how the many-body interaction affects the calculated spectrum. This information can be easily obtained within the PBC-GEBF framework. We optimized the crystal structure and calculated the IR spectrum with the many-body interaction truncated up to three-body (γmax = 3) and two-body (γmax = 2), respectively. The obtained IR spectrum is compared to that from our standard GEBF calculation, where γmax is set to be 4 by default. Calculated IR spectra are displayed in Figure 6. One may find that the

III.C. Raman Spectra of the Ammonia Borane Crystal. The crystalline ammonia borane (BH3NH3) has been studied extensively by both experimental and theoretical chemists,2,38−40 due to its potential applications as a hydrogen storage system. DFT methods are most commonly employed in theoretical studies. Nevertheless, theoretical predictions on vibrational spectra have not yet performed with post-HF theories. In this work, we apply the PBC-GEBF-MP2 method to investigate the crystal structure and Raman spectra of the BH3NH3 crystal, which belongs to space group Pmn21 at low temperature. We use the structure obtained from the neutron diffraction38 as the initial geometry and fully optimize the structure of the ammonia borane crystal. The optimized lattice parameters are summarized in Table 2. One can see that the optimized structure from PBC-GEBF-MP2 calculation is in good accord with the experimental structure. Furthermore, the calculated and experimental Raman spectra for both intramolecular vibrational modes and lattice modes are collected in Table 3 and Table 4, respectively. Table 3. Vibrational Frequencies for Intramolecular Vibrational Modes (in cm−1) of the BH3NH3 Crystala molecule

crystal calc

exp

calc

exp

calc

modes

A1

3337

3330

3255

3254

−82

−76

A1

2340

2363

2280

2304

−60

−59

A1

1343

1291

1375

1370

32

79

A1

1175

1164

1158

1220

−17

56

A1

968/785

646

782

766

E

3386

3449

3316

3342

−70

−107

E

2415

2425

2375

2379

−40

−46

E

1608

1586

1596

1572

−12

−14

E

1186

1168

1185

1149

−1

−19

E

1052

1035

1066

1049

14

14

E

603

624

726

715

123

90

sym N−H stretch sym B−H stretch sym N−H bend sym B−H bend B−N stretch B−N torsion asym N− H stretch asym B−H stretch asym N− H bend asym B−H bend asym N− H + B− H rock asym N− H + B− H rock

sym

exp

b

shift

b

120

A2

a

The corresponding vibrational frequencies of an isolated BH3NH3 molecule are also listed for comparison. bReference 38. Figure 6. Comparison of IR spectra obtained from PBC-GEBF-MP2 calculations with different many-body interaction truncations. Part (a) displays the IR spectra in the region from 500 to 2000 cm−1, and part (b) shows the IR spectra in the region from 3000 to 4000 cm−1.

For the BH3NH3 molecule, there are 18 vibrational modes, 6 of which belong to one-dimensional representation (A1 and A2), while the other 12 vibrational modes belong to twodimensional representation (E), which correspond to 6 doubly degenerate frequencies. In the BH3NH3 crystal, each unit cell consists of two molecules; hence there are 36 intramolecular vibrational modes. It should be noted that the degenerate frequencies may split to some extent because the symmetry of molecule is not maintained in the crystal structure. Therefore, in Table 3, the frequency with larger intensity for each mode is displayed. It can be observed that for most vibrational modes, the deviation of calculated frequencies with respect to the

spectrum calculated with γmax = 3 is very close to that obtained with γmax = 4. However, the IR spectrum calculated with γmax = 2 deviates considerably from the spectrum obtained with γmax = 4, both in low frequency and in high frequency regions, which are shown in Figure 6a and b, respectively. This result indicates that fragment-based approaches have to include the trimer interaction for predicting a satisfactory vibrational spectrum for this system. 2707


Article

The Journal of Physical Chemistry A Table 4. Lattice Mode Vibrational Frequencies (in cm−1) for the BH3NH3 Crystal PBC-PBEa

a

PBC-GEBF-MP2

sym

frequency

intensity

frequency

intensity

A2 A2 A1 B1 B2 A1 B1 A2 B2 B2 A2

77.2 92.9 103.4 143.8 154.2 167.0 173.5 216.4 237.8 320.9 343.3

0.2312 0.1550 0.0312 0.0004 0.1877 0.0849 0.0118 0.2705 0.0572 0.0577 0.0798

97.2 99.2 120.2 119.0 130.5 190.9 193.6 218.1 222.8 337.3 354.1

0.0926 0.1050 0.0063 0.0001 0.3658 0.1172 0.1117 1.2245 0.2390 0.0820 1.1342

which are also available from calculations on subsystems. In comparison with the periodic quantum chemistry method, the PBC-GEBF method provides an appealing alternative approach to compute the vibrational spectra of molecular crystals with much reduced computational costs, especially when calculations are done at the post-Hartree−Fock levels. By comparing the results from the PBC-GEBF method and the methods based on periodic electronic structure theory for two typical molecular crystals (CO2 and imidazole), we demonstrate that the PBC-GEBF approach can reproduce the results in predicting vibrational spectra of molecular crystals. For the hydrogen-bond imidazole crystal, our analysis shows that the inclusion of the long-range electrostatic interaction (from the crystal environment) is critical for accurate descriptions of its vibrational spectra. We then employed the PBC-GEBF approach to investigate the IR spectra of the urea and BH3NH3 crystals at the MP2 level. For the urea crystal, our calculations show that the PBC-GEBF-MP2 method give satisfactory predictions on the observed IR spectra. We have assigned the main peaks in the experimental IR spectra with the calculated results. For the BH3NH3 crystal, the calculated vibrational spectra from PBC-GEBF-MP2 are also in very good agreement with the observed Raman spectra. In conclusion, one can expect that the PBC-GEBF approach, combined with accurate molecular quantum chemistry methods, becomes an efficient theoretical tool for computing vibrational spectra of various types of molecular crystals.

experimenta frequency ≃100 116 150 159 185 212 221 ≃337 ≃337

Reference 39.

experimental values are less than 20 cm−1. In addition, we also performed vibrational frequency calculations for an isolated BH3NH3 at the MP2/6-311++G(d,p) level. The frequency shifts for all vibrational modes in the molecular crystal are also listed. Red-shift can be observed for both B−H and N−H stretching modes. The shifts in N−H stretching modes are larger as compared to shifts in B−H stretching modes, which are in good accord with the experiment values. One may also notice a remarkable blue-shift (90 cm−1) in asymmetrical N−H plus B−H rock mode, which is in accord with the experimental observation (123 cm−1). A comparison of calculated and observed frequency shifts suggests that the PBC-GEBF-MP2 approach is capable of providing a satisfactory description on intramolecular vibrational modes. Except for the intramolecular vibrational modes, there are also 9 intermolecular modes in a unit cell. These 9 intermolecular modes, together with two intramolecular modes with lowest frequencies, are denoted as lattice vibrational modes. For these lattice modes, the calculated values are compared to the experimental data obtained at low temperature (T = 15 K).40 For most lattice modes, the deviations between calculated and experimental vibrational frequencies are a few wavenumbers. Nevertheless, the largest deviation (about 30 cm−1) appears at the vibrational mode at 159 cm−1. For the last mode in the low frequency region, a large deviation (about 17 cm−1) may result from the unclear experimental assignment, because a broad band is observed in the experimental spectra even at low temperature. In summary, the results displayed in Table 3 and Table 4 show that the PBC-GEBF-MP2 method provides satisfactory predictions for vibrational spectra of the crystalline BH3NH3 system.

■

APPENDIX A: THE DERIVATION OF THE FORMULA FOR THE UNIT CELL ENERGY OF A MOLECULAR CRYSTAL The formula for evaluating the unit cell energy of a periodic molecular crystal in the PBC-GEBF method can be derived from the expression of the total energy of a nonperiodic molecular system in the GEBF method, as shown below: E Tot =

M

⎛

M

⎞

m

⎝

m

⎠

∑ CmEm̃ − ⎜⎜∑ Cm − 1⎟⎟ ∑ ∑ A

B>A

Q AQ B RAB

(A1)

Naturally, the unit cell energy can be defined as the total energy of a supercell containing N3 unit cells divided by the number of unit cells (N3) when N goes to infinity, which is expressed as Eunit‐cell

⎡ N3 M ⎛ N3 M ⎞ 1 ⎢ = lim 3 ∑ ∑ Cm , nEm̃ , n − ⎜⎜∑ ∑ Cm , n − 1⎟⎟ N →∞ N ⎢ ⎝ n m ⎠ ⎣ n m

∑ (A > B) ∈ N3

IV. CONCLUSIONS We have implemented the calculation of vibrational spectra of molecular crystals within the PBC-GEBF approach. In the PBC-GEBF approach, the dynamical matrix of a molecular crystal can be evaluated by combing the mass-weighted Hessian matrixes of small embedded subsystems, which can be routinely calculated with existing quantum chemistry softwares. The diagonalization of this dynamical matrix can lead to vibrational frequencies and normal modes of a molecular crystal. In a similar way, the IR and Raman intensities of these normal modes can also be evaluated with the derivatives of the dipole moment and polarizability with respect to atomic coordinates,

⎤ Q AQ B ⎥ RAB ⎥⎦

(A2)

It should be mentioned that in eq A2 the periodic boundary condition is implicitly adopted so that each generated subsystem has N3 replicas due to the translation symmetry. In fact, here a subsystem means an electrostatically embedded subsystem in the presence of an infinite number of point charges generated by atoms in all other unit cells. Obviously, the direct calculation for such a subsystem (including the effect of an infinite number of point charges into the Hamiltonian of a subsystem) is difficult. To avoid this difficulty, we introduce the following relationship (following the idea of the ONION method): 2708


The Journal of Physical Chemistry A Em̃ , n ≈ Em̃ ′ , n −

∑ 3

A ∈ K (B > A) ∈ K

+

∑

∑

A ∈ N3 (B > A) ∈ N3

■

Q AQ B

∑

APPENDIX B: THE MANY-BODY EXPANSION OF THE UNIT CELL ENERGY OF A MOLECULAR CRYSTAL The total energy of a molecular cluster can be decomposed as

RAB

3

Q AQ B RAB

(A3)

N

E=

Here, Em̃ ′ , n stands for the energy of a subsystem embedded in the field of a finite array of point charges. In another word, now a subsystem is only embedded in the presence of a finite array of point charges within a supercell (containing neighboring K3 unit cells, while K is finite). This electrostatically embedded subsystem can be treated with a given quantum chemistry method. The physical approximation behind eq A3 is that the effect of short-range electrostatic interaction to the subsystem (due to point charges within neighboring unit cells) is treated quantum mechanically, and the effect of long-range electrostatic interaction (due to point charges from all distant unit cells) is treated classically. It should be mentioned that the above relationship works well only if the overall dipole moment of each unit cell is zero. To meet this requirement, we can add the compensating charges to each unit cell in the construction of supercell so that the overall dipole moment of each unit cell is zero. With the compensating field, one can use a moderate supercell (as discussed in the text). Combining eq A3 and eq A2, we can get Eunit‐cell

∑

RAB

(A > B) ∈ K3

∑

+

(A > B) ∈ N3

N

Eunit‐cell =

+ (A4)

m

+ lim

N →∞

⎛

M

⎝

m

⎞

Q AQ B

⎛

M

⎝

m

+ E Ewald‐sum

(B3)

N

N

N

Lat

N

N

i

j>i

N

L>0 Lat

N

N

i

j

N

∑ ∑ ∑ Ei0,j0,k0 + ∑ ∑ ∑ ∑ Ei0,j0,kL j>i k>j Lat

N

L>0 N

i

j>i

k

N

∑ ∑ ∑ ∑ ∑ EiM ,j0,kL ... i

j

k

(B4)

Figure 7. An illustrative scheme of a one-dimensional periodic system, in which each unit cell is composed of two fragments.

(A7)

⎞

∑ CmEm̃ − ⎜⎜∑ Cm⎟⎟ ∑ ∑ m

(B2)

⎠ (A > B) ∈ K3 RAB

⎡ Q AQ B ⎤ 1 ⎢ ⎥ ∑ N3 ⎢⎣(A > B) ∈ N3 RAB ⎥⎦

M

(B1)

j>i k>j

in which i, j, k,... represent fragments in a given system, and M (or L) represents the cell index of neighboring cells around the central unit cell. Lattice index 0 refers to the central unit cell. It should be mentioned that the overcounting of all energy terms is avoided by limited summations. To further illustrate the scheme of the energy decomposition, we use a one-dimensional periodic system as an example. Each unit cell is composed of two fragments, which is shown in Figure 7. The approximate unit cell energy, which is truncated

The last term in eq A7 is just the classical electrostatic interaction within the central cell and between the central cell and all image cells, which can be calculated with the Ewald summation method. Eunit‐cell =

i

∑ Ei0 + ∑ ∑ Ei0,j0 + ∑ ∑ ∑ Ei0, jL

M0

Substituting the above relationships into eq A4 and extracting the first two terms out of the bracket, we can obtain M

j>i

and so on and so forth. For a periodic molecular crystal, the unit cell energy can also be expanded as the summation of the following energy terms:

Lat

(A6)

N

E(ijk) = Ei + Ej + Ek + Ei , j + Ei , k + Ej , k + Ei , j , k

+

Em̃ ′ ,1 = Em̃ ′ ,2 = Em̃ ′ ,3 = ... = Em̃ ′ , N3

∑ CmEm̃ ′ − ⎜⎜∑ Cm⎟⎟ ∑

i

i

(A5)

N

Similarly, for a system with three fragments, the energy can be expressed as

N

Cm ,1 = Cm ,2 = Cm ,3 = ... = Cm , N3

N

E(ij) = Ei + Ej + Ei , j

Using the translation symmetry in periodic systems, one may naturally find

Eunit‐cell =

N

in which Ei represents an intrafragment energy for the ith monomer, and Ei,j (or Ei,j,k) denotes the two-fragment (or three-fragment) interaction energy. Moreover, E(ij) (or E(ijk)) denotes the two-fragment (or three-fragment) total energy. For instance, the energy of a system with two fragments can be expressed as

i

⎤ Q AQ B ⎥ RAB ⎥⎦

N

∑ Ei + ∑ ∑ Ei ,j + ∑ ∑ ∑ Ei ,j ,k + ... i=1

⎡ N3 M ⎞ ⎛ N3 M 1 ⎢ = lim 3 ∑ ∑ Cm , nEm̃ ′ , n − ⎜⎜∑ ∑ Cm , n⎟⎟ N →∞ N ⎢ ⎠ ⎝ n m ⎣ n m Q AQ B

Article

⎠

A ∈ K (B > A) ∈ K

here up to three-body interactions among three nearest neighboring fragments, can be expressed as

Q AQ B

Eunit‐cell ≈ E1 + E2 + E1,2 + E2,1 ″ + E1,1 ″ + E2 ′ ,2 + E1,2,1 ″

RAB

+ E2 ′ ,1,2

(A8)

(B5)

According to eq B2 and eq B3, the energies of different molecular clusters can be expanded as

In summary, with the translation symmetry of a periodic molecular crystal, we can extend the GEBF method from molecular systems to periodic molecular crystals and obtain the energy expression for a unit cell as in eq A8.

E(2′12) ≈ E2 ′ + E1 + E2 + E1,2 + E2 ′ ,1 + E2 ′ ,2 + E2 ′ ,1,2 (B6) 2709



■

E(121″) ≈ E1 + E2 + E1 ″ + E1,2 + E2,1 ″ + E1,1 ″ + E1,2,1 ″ (B7)

E(12) ≈ E1 + E2 + E1,2

(B8)

E(2′1) ≈ E2 ′ + E1 + E2 ′ ,1

(B9)

(B10)

By comparing energy components of eq B5 with those of eqs B6−B10, we then can easily find that the unit cell energy can be approximated by Eunit‐cell ≈ ( +1) × E(2′12) + ( +1) × E(121″) + ( −1) × E(12) + ( − 1) × E(2′1)

(B11)

Hence, the unit cell energy of a periodic system can be approximately obtained as the combination of energies of a set of nonperiodic subsystems (or molecular clusters), which can be calculated with existing molecular quantum chemistry methods. According to the definition of subsystem described in part A of section II, two subsystems (2′12) and (121″) are called primitive subsystems with coefficient of +1. The other two subsystems (12) and (2′1) are called derivative subsystems with coefficient of −1, as shown in Figure 8.

Figure 8. An illustrative scheme of primitive and derivative subsystems generated from a one-dimensional periodic system.

For more complex three-dimensional periodic systems, a general procedure of constructing subsystems can be proposed. This procedure is briefly described in part A in section II. It should also be mentioned that long-range interactions are taken into account by the embedding charges in the super cell and compensating point charges at the boundaries of the super cell.

■

REFERENCES

(1) Rousseau, B.; Van Alsenoy, C.; Keuleers, R.; Desseyn, H. O. Solids Modeled by Ab Initio Crystal Field Methods. Part 17. Study of the Structure and Vibrational Spectrum of Urea in the Gas Phase and in Its P4̅21m Crystal Phase. J. Phys. Chem. A 1998, 102, 6540−6548. (2) Dillen, J.; Verhoeven, P. The End of a 30-Year-Old Controversy? A Computational Study of the B−N Stretching Frequency of BH3− NH3 in the Solid State. J. Phys. Chem. A 2003, 107, 2570−2577. (3) Deev, V.; Collins, M. A. Approximate Ab Initio Energies by Systematic Molecular Fragmentation. J. Chem. Phys. 2005, 122, 154102. (4) Netzloff, H. M.; Collins, M. A. Ab Initio Energies of Nonconducting Crystals by Systematic Fragmentation. J. Chem. Phys. 2007, 127, 134113. (5) Collins, M. A. Ab Initio Lattice Dynamics of Nonconducting Crystals by Systematic Fragmentation. J. Chem. Phys. 2011, 134, 164110. (6) Hirata, S. Fast Electron-correlation Methods for Molecular Crystals: An Application to the α, β1, and β2 Modifications of Solid Formic Acid. J. Chem. Phys. 2008, 129, 204104. (7) Sode, O.; Hirata, S. Second-order Many-body Perturbation Study of Solid Hydrogen Fluoride under Pressure. Phys. Chem. Chem. Phys. 2012, 14, 7765−7779. (8) He, X.; Sode, O.; Xantheas, S. S.; Hirata, S. Second-order Manybody Perturbation Study of Ice Ih. J. Chem. Phys. 2012, 137, 204505. (9) Hirata, S.; Gilliard, K.; He, X.; Li, J.; Sode, O. Ab Initio Molecular Crystal Structures, Spectra, and Phase Diagrams. Acc. Chem. Res. 2014, 47, 2721−2730. (10) Beran, G. J. O.; Nanda, K. Predicting Organic Crystal Lattice Energies with Chemical Accuracy. J. Phys. Chem. Lett. 2010, 1, 3480− 3487. (11) Wen, S.; Nanda, K.; Huang, Y.; Beran, G. J. O. Practical Quantum Mechanics-based Fragment Methods for Predicting Molecular Crystal Properties. Phys. Chem. Chem. Phys. 2012, 14, 7578−7590. (12) Heit, Y.; Beran, G. J. O. Exploiting Space-group Symmetry in Fragment-based Molecular Crystal Calculations. J. Comput. Chem. 2014, 35, 2205−2214. (13) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701−706. (14) Fujita, T.; Nakano, T.; Tanaka, S. Fragment Molecular Orbital Calculations under Periodic Boundary Condition. Chem. Phys. Lett. 2011, 506, 112−116. (15) Rice, B. M.; Szalewicz, K.; Podeszwa, R. Predicting Structure of Molecular Crystals from First Principles. Phys. Rev. Lett. 2008, 101, 115503. (16) Yang, J.; Hu, W.; Usvyat, D.; Matthews, D.; Schuetz, M.; Chan, G. K. Ab Initio Determination of the Crystalline Benzene Lattice Energy to Sub-kilojoule/mol Accuracy. Science 2014, 345, 640−643. (17) Bygrave, P. J.; Allan, N. L.; Manby, F. R. The Embedded Manybody Expansion for Energetics of Molecular Crystals. J. Chem. Phys. 2012, 137, 164102. (18) Zhang, P.; Truhlar, D. G.; Gao, J. Fragment-based Quantum Mechanical Methods for Periodic Systems with Ewald Summation and Mean Image Charge Convention for Long-range Electrostatic Interactions. Phys. Chem. Chem. Phys. 2012, 14, 7821−7829. (19) 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. (20) 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. (21) 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. (22) Hua, W.; Fang, T.; Li, W.; Yu, J.; Li, S. Geometry Optimizations and Vibrational Spectra of Large Molecules from a Generalized

If the translation symmetry is employed, we have E2 ′ = E2E1 ″ = E1

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (Grant nos. 21333004 and 21361140376) and the National Basic Research Program (Grant no. 2011CB808501). Part of the calculations was performed at the High Performance Computing Center of Nanjing University and the Shenzhen Supercomputer Center (SSC, China). 2710


Article

The Journal of Physical Chemistry A Energy-Based Fragmentation Approach. J. Phys. Chem. A 2008, 112, 10864−10872. (23) Wang, K.; Li, W.; Li, S. Generalized Energy-Based Fragmentation CCSD(T)-F12a Method and Application to the Relative Energies of Water Clusters (H2O)20. J. Chem. Theory Comput. 2014, 10, 1546−1553. (24) 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. (25) 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. (26) Toukmaji, A. Y.; Board, J. A. Ewald Summation Techniques in Perspective: A Survey. Comput. Phys. Commun. 1996, 95, 73−92. (27) Person, W. B.; Zerbi, G. Vibrational Intensities in Infrared and Raman Spectroscopy; Elsevier Scientific Publishing Co.: Amsterdam/ Oxford/New York, 1982. (28) 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 D.01; Gaussian, Inc.: Wallingford, CT, 2009. (29) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; et al. CRYSTAL14; University of Torino: Torino, 2014. (30) Pascale, F.; Zicovich-Wilson, C. M.; López Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. The Calculation of the Vibrational Frequencies of Crystalline Compounds and Its Implementation in the CRYSTAL code. J. Comput. Chem. 2004, 25, 888−897. (31) Zicovich-Wilson, C. M.; Pascale, F.; Roetti, C.; Saunders, V. R.; Orlando, R.; Dovesi, R. Calculation of The Vibration Frequencies of αquartz: The Effect of Hamiltonian and Basis Set. J. Comput. Chem. 2004, 25, 1873−1881. (32) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158. (33) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (34) Merrick, J. P.; Moran, D.; Radom, L. An Evaluation of Harmonic Vibrational Frequency Scale Factors. J. Phys. Chem. A 2007, 111, 11683−11700. (35) Keuleers, R.; Desseyn, H. O.; Rousseau, B.; Van Alsenoy, C. Vibrational Analysis of Urea. J. Phys. Chem. A 1999, 103, 4621−4630. (36) Grdadolnik, J.; Maréchal, Y. Urea and Urea−water Solutions: An Infrared Study. J. Mol. Struct. 2002, 615, 177−189. (37) Rousseau, B.; Keuleers, R.; Desseyn, H. O.; Geise, H. J.; Van Alsenoy, C. Solids Modeled by Ab-initio Crystal Field Methods. Effects of Intermolecular Interactions on the Vibrational Spectrum of Urea. Chem. Phys. Lett. 1999, 302, 55−59. (38) Klooster, W. T.; Koetzle, T. F.; Siegbahn, P. E. M.; Richardson, T. B.; Crabtree, R. H. Study of the N−H···H−B Dihydrogen Bond Including the Crystal Structure of BH3NH3 by Neutron Diffraction. J. Am. Chem. Soc. 1999, 121, 6337−6343. (39) Custelcean, R.; Dreger, Z. A. Dihydrogen Bonding under High Pressure: A Raman Study of BH3NH3 Molecular Crystal. J. Phys. Chem. B 2003, 107, 9231−9235. (40) Ziparo, C.; Colognesi, D.; Giannasi, A.; Zoppi, M. Raman Spectra of Ammonia Borane: Low Frequency Lattice Modes. J. Phys. Chem. A 2012, 116, 8827−8832.

2711


Vibrational Spectra of Molecular Crystals with the Generalized Energy

Recommend Documents