WWW.C-CHEM.ORG
FULL PAPER
Acceleration of Self-Consistent Field Convergence in Ab Initio Molecular Dynamics Simulation with Multiconfigurational Wave Function Masaki Okoshi,[a] and Hiromi Nakai*[a,b,c,d] The Lagrange interpolation of molecular orbital (LIMO) method, which reduces the number of self-consistent field iterations in ab initio molecular dynamics simulations with the Hartree–Fock method and the Kohn–Sham density functional theories, is extended to the theory of multiconfigurational wave functions. We examine two types of treatments for the active orbitals that are partially occupied. The first treatment, as denoted by LIMO(C), is a simple application of the conventional LIMO method to the union of the inactive core and the
active orbitals. The second, as denoted by LIMO(S), separately treats the inactive core and the active orbitals. Numerical tests to compare the two treatments clarify that LIMO(S) is superior to LIMO(C). Further applications of LIMO(S) to various systems C 2014 Wiley demonstrate its effectiveness and robustness. V Periodicals, Inc.
Introduction
coefficient matrices of a few previous MD steps. The technique, called LIMO, drastically reduced the SCF iterations, with observed reductions ranging from 50 to 70%. They further extended the LIMO method to geometry optimization and Monte Carlo simulations, which have no information on the simulation (or physical) time.[19] Instead, a least-squares (LS) technique utilizing geometry information was applied to MO propagation (LSMO), and the performance of this technique was comparable to the LIMO technique in AIMD simulations. The single-configurational (SC) theories, including HF and DFT, and the single-reference (SR) theories, such as the Mïller– Plesset perturbation and the coupled-cluster theories, successfully reproduce potential energies and forces acting on atoms in most of ground-state systems. However, especially in quasi-
The ab initio molecular dynamics (AIMD) method, which utilizes ab initio potential energies and forces, has become a useful tool to investigate chemical reaction dynamics.[1] Its major drawback, however, is the high computational cost of the ab initio calculations performed at each molecular dynamics (MD) step. One of the computational bottlenecks in Hartree–Fock (HF) and Kohn–Sham (KS) density functional theory (DFT) calculations is the convergence of the electronic wave function onto the Born–Oppenheimer potential energy surface, which requires the iterative self-consistent field (SCF) procedure. Utilizations of fairly reasonable approximations, such as the linear-scaling technique[2,3] and the dual basis set method,[4] for instance, are capable of reducing computational efforts without the collapse of simulations. Alternatively, several techniques that propagate the electronic structure instead of performing the SCF procedure have been extensively investigated.[5–13] As the SCF convergence in the standard HF/ KS DFT calculations has already been well optimized through methods such as direct inversion in the iterative subspace[14] and its expansions,[15] it might be difficult to achieve further improvements. However, it is possible to obtain information regarding the electronic structure in the AIMD simulation, and this information could potentially be used to further optimize the convergence calculations. Pulay and Fogarasi[16] proposed the Fock matrix dynamics (FMD) method, which extrapolates the Fock matrix forward in time. In their numerical tests, the FMD technique accelerated the SCF convergence by about 50%. Herbert and Head-Gordon[17] further assessed the performances of FMD, in the viewpoints of the efficiencies and accuracies in energy drift. Atsumi and Nakai[18] developed a molecular orbital (MO) propagation technique, which proposes an efficient MO initial guess by the Lagrange interpolation (LI) of MO
DOI: 10.1002/jcc.23617
[a] M. Okoshi, H. Nakai Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan E-mail: [email protected] [b] H. Nakai Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan [c] H. Nakai CREST, Japan Science and Technology Agency, Saitama 332-0012, Japan [d] H. Nakai ESICB, Kyoto University, Kyotodaigaku-Katsura, Kyoto 615-8520, Japan Contract grant sponsor: Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST; Theoretical Design of Materials with Innovative Functions Based on Relativistic Electronic Theory); Contract grant sponsor: Elements Strategy Initiative to Form Core Research Center, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan [Elements Strategy Initiative for Catalysts & Batteries (ESICB)]; Contract grant sponsor: Computational Materials Science Initiative (CMSI), MEXT, Japan [Strategic Programs for Innovative Research (SPIRE)]; Contract grant sponsor: Grants for Excellent Graduate Schools, MEXT, Japan (Practical Chemical Wisdom). C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
1
FULL PAPER
WWW.C-CHEM.ORG
degenerate systems, multiconfigurational (MC) or multireference (MR) type wave function theories should be utilized: for instance, an ozone molecule, multibond formation/cleavage systems, and nonadiabatic regions, in which two or more electronic states are (pseudo-) crossing each other. While the MC/ MR wave function theories generally need much more computational efforts than the SC/SR wave function theories, AIMD with the MC/MR wave function theories have been extensively performed, particularly for the systems, in which photochemical reactions with nonradiative transitions take place.[20–26] In the present article, we extended the LIMO method, which was originally developed for the SC wave function, to the MC wave function theories, that is, the complete active space SCF (CASSCF) and restricted active space SCF (RASSCF) methods. Theoretical considerations were initially made and numerical tests were performed.
Theoretical Aspects LIMO for SC wave functions This subsection briefly reviews the original LIMO method developed for the SC wave function.[18] In AIMD, the conðnÞ verged MO coefficient matrix of the last (nth) step, Cconv , is often used as the initial guess for the (n 1 1)th step because the geometric and electronic structures of the neighboring MD steps are expected to be similar to each other. We refer to this conventional technique as PREMO, to refer to the fact that only the MO coefficient matrix of the previous step is used for the initial MO guess. The LIMO technique predicts better initial MO guesses than PREMO, by using the converged MO coefficient matrices in the last step as well as in several previous steps. The main concept of LIMO is that the MO coefficient matrices are expected to change smoothly because the orbital centers move continuously in the AIMD simulations. Thus, an effective initial MO guess can be predicted as the linear combination of the several converged MO coefficient matrices. The predicted initial ðn11Þ MO guess for the (n 1 1)th step, Cprd , is given by the LI technique:
ðn21Þ!n
TpSCF
21 t ðn21Þ n ðn21Þ ðn21Þ 5 t Cconv Cconv Cconv Cconv
(3)
and ðn2lÞ!n
TpSCF
5
l21 Y
ðn2l1mÞ!ðn2l1m11Þ
TpSCF
(4)
m50
where the subscript “pSCF” refers to the pseudo-SCF effect that is present in the (n-l)th and nth converged MO coefficient matrices. In the case of the PREMO method, the transformaðn21Þ!n tion matrix for the neighboring steps, TpSCF , represents the SCF effect in the nth MD step. Note that the pSCF effect is considered for the occupied orbitals only; that is, the transformation matrix T has the dimensions of Nocc 3 Nocc, and the MO coefficient matrix C and its transpose tC have the dimensions of Nbasis 3 Nocc and Nocc 3 Nbasis, respectively, where Nbasis and Nocc are the numbers of the basis functions and of the occupied orbitals. Mathematical proofs have shown that the orbital-invariant LIMO method for all occupied and unoccupied orbitals is completely meaningless, since the effective ðn11Þ MO initial guess, Cprd , becomes identical to the converged [18] MO coefficient matrix of the last step, CðnÞ conv . The effectiveness of the MO initial guess predicted by LIMO depends on the degree of the Lagrange polynomial (k). Numerical assessments have shown the existence of the optimal degree kopt, which depends mainly on the time step of AIMD, Dt.[18] While it is unfeasible to determine kopt prior to performing SCF, applying the optimal degree at the nth step to the (n 1 1)th step results in acceptable performance. Therefore, the optimal degree kopt is determined to minimize the following deviation sum: k X X X ðn2lÞ!n n ðn2lÞ dk 5 Lk;l Cconv ðl; jÞTpSCF ðj; iÞ Cconv ðl; iÞ2 l;i j l50
(5)
where l runs over the basis functions, whereas i and j run over the occupied MOs.
Extension of LIMO to the MC wave function ðn11Þ
Cprd 5
k X
ðn2lÞ!n
Lk;l Cðn2lÞ conv TpSCF
(1)
l50
where fCðn2lÞ conv g (l 5 0, 1, 2, . . ., k) are the converged MO coefficient matrices at the (n-l)th step, and fLk;l g are the lth Lagrange coefficients for the kth degree of polynomial given by Lk;l 5
k Y t ðn11Þ 2tðn2mÞ m50
tðn2lÞ 2t ðn2mÞ
(2)
When adopting a constant time step in AIMD, the Lagrange coefficients become simple integers. The transformation ðn2lÞ!n matrix, fTpSCF g, leads to the orbital-invariant formula, by considering the orbital crossing and mixing effects within the ðn2lÞ!n occupied orbitals. The explicit formula of TpSCF is given by: 2
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
As described above, pSCF can only be considered for the occupied orbitals; therefore, it is indispensable to separate the occupied and the unoccupied orbital spaces. In the case of the SC wave function, the occupied/unoccupied orbitals are well defined. However, in the MC wave function case, there remains ambiguity to define the occupied orbitals, since the active orbitals are partially occupied. Therefore, the main problem is how to treat the pSCF effect. In the present study, we examined two types of effective MO initial guesses of dimensions Nbasis 3 Nocc, where Nocc refers to the summation of the number of the inactive core and the active orbitals. One method to treat the pSCF effect is the conventional LIMO scheme for the occupied orbitals, that is, the union of the inactive core and the active orbitals. In this case, all occupied orbitals are treated using eqs. (1)–(5) in the similar procedure as the SC case, as if they are doubly occupied orbitals. WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 1. Initial geometries of (a) HCNH1 and (b) benzene.
Another method is to treat the inactive core and the active orbitals separately and independently: ðn11Þ
Cprd; core 5
k X
ðn2lÞ!n
Lk;l Cðn2lÞ conv; core TpSCF; core
(6)
l50 ðn11Þ
Cprd; active 5
k X
ðn2lÞ
ðn2lÞ!n
Lk;l Cconv; active TpSCF; active
(7)
l50
where the subscripts “core” and “active” refer to the inactive core and the active orbitals, respectively. In this method, the matrices Ccore, Tcore, Cactive, and Tactive have the dimensions of Nbasis 3 Ncore, Ncore 3 Ncore, Nbasis 3 Nactive, and Nactive 3 Nactive, respectively, where Ncore and Nactive refer to the numbers of the inactive core and the active orbitals. The total predicted MO coefficient matrix is formed as the direct sum of Ccore and Cactive: ðn11Þ
ðn11Þ
ðn11Þ
Cprd 5Cprd; core 丣Cprd; active
(8)
ðn11Þ
FULL PAPER
on the atoms were evaluated using the CASSCF method with Pople’s double-zeta-plus-polarization type basis set, 6-31G**.[28] The used active spaces were 10 electrons within 10 orbitals and six electrons within six orbitals as denoted by 10e/10o and 6e/6o for HCNH1 and benzene, respectively. The active space included full valence electrons/orbitals in the case of HCNH1, while all p electrons/orbitals were used as the active space for benzene. The degrees of the Lagrange polynomial k were determined automatically using eq. (5), with a maximum prediction degree kmax of three. The Velocity Verlet method with a time step (Dt) of 0.1 fs was adopted for the numerical integration of the Newtonian equation of motion of nuclei under the microcanonical (NVE) ensemble. The AIMD simulations were performed for 300 steps (30.0 fs) after 50 steps (5.0 fs) of equilibration, for each system. Randomly strained structures from the equilibrium with no velocities, which are shown in Figures 1a and 1b for HCNH1 and benzene, respectively, were used as the initial conditions. The distortion energies for HCNH1 and benzene were 1.05 and 7.44 eV, whereas their zero-point vibrational energies (ZVEs) were 0.75 and 2.92 eV, respectively. ZVEs were calculated by summing up the normal mode frequencies, calculated for the equilibrium structures at the same levels of theory as the AIMD simulations; that is, CASSCF(10e/10o)/6-31G** and CASSCF(6e/6o)/6-31G** for HCNH1 and benzene, respectively. Figures 2a and 2b show the energy variations in the AIMD simulations of HCNH1 and benzene, respectively. The dotted, dashed, and solid lines correspond to the results for the kinetic, potential, and total energies of the system. The potential and total energies are relative values to those at t 5 25.0 fs, that is, the initial conditions. The simulations at t 5 25.0–0.0 fs correspond to the equilibration, whereas those
The initial MO guess, Cprd , is orthonormalized prior to the MCSCF procedure through Schmidt’s procedure in all cases, including PREMO, the conventional LIMO method, and the LIMO method for the separated orbitals. For convenience, the conventional scheme is denoted by LIMO(C) in this article, whereas the separated procedure given in eq. (8) is denoted by LIMO(S). The optimal degree of the Lagrange polynomial is determined in a manner similar to the conventional LIMO through use of eq. (5) for the total MO coefficient matrices.
Numerical Assessments Numerical performances of the LIMO technique were investigated further. AIMD simulations, as well as MO prediction procedures, are typically performed using an original MD code combined with quantum chemistry program packages that use Gaussian basis functions. In the present AIMD simulations, the Gaussian 03 suite of programs[27] was utilized to perform ab initio MCSCF calculations. Numerical tests for prediction methods The performances of the LIMO(C) and LIMO(S) prediction schemes were compared by AIMD simulations for HCNH1 and benzene molecules in the S0 ground state. The forces acting
Figure 2. Energy variations in the AIMD simulations of (a) HCNH1 and (b) benzene. The dotted, dashed, and solid lines correspond to the kinetic, potential, and total energies, respectively.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 3. Changes in (a) the number of SCF iterations NSCF and (b) the energy difference DE over time for HCNH1. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
at t 5 0.0–30.0 were considered as the production runs. For HCNH1 (Fig. 2a), the kinetic energy (dotted line) rapidly increased to about 1.0 eV at 25.0–0.0 fs, whereas the potential energy (dashed line) decreased as much. This meant that a part of distortion energy changed into the kinetic energy with the geometric relaxation. After that, the kinetic energy as well as the potential energy oscillated with a period of about 7.0– 9.0 fs. The total energy (solid line) was approximately zero, with the maximum absolute value of 0.53 meV. As for benzene (Fig. 2b), the similar tendency was observed with the large kinetic energy of about 5.0 eV. The maximum absolute value for the total energy was 2.8 meV. Figures 3a and 3b show the number of SCF iterations (NSCF) and the energy difference DE, respectively, over time for the AIMD simulations of HCNH1. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. The energy difference DE was used to assess the performance of each prediction method: ðnÞ
ðnÞ DE5jE½Cprd; k 2E½Cconv j
(9)
where the two terms on the right-hand side of the equation indicate the potential energies evaluated using the MO coeffiðnÞ cient matrices of Cprd; k and CðnÞ conv , respectively. For the PREMO method (dotted line), NSCF slightly oscillated between 4 and 6 cycles, whereas DE ranged from 1.3 3 1026 to 1.4 3 1025 hartrees. In the case of LIMO(C) (dashed line), DE was observed to be lower compared to the results from the PREMO method, with DE ranging from 2.5 3 1029 to 4.5 3 1027 hartrees. The values for NSCF were also relatively low, ranging from 2 to 4 cycles. In both PREMO and LIMO(C), oscil4
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
lations of DE with a period of 7.0–9.0 fs were observed. These correspond to the oscillation of the kinetic energy of the system, as shown in Figure 2a. On the contrary, unphysical fluctuations with high amplitudes, ranging from 2.9 3 10211 to 2.7 3 1027 hartrees, were observed in the case of LIMO(S). Despite the existence of the unphysical fluctuations, DE for LIMO(S) was generally lower than that of LIMO(C) with exceptions at several timings of t at 3.4, 3.9, 10.8, 25.4, 25.5, 25.9–26.3, 26.5–26.8, 27.4, and 27.5 fs, and the NSCF value was comparable to that of LIMO(C) at 2–4 cycles. The unphysical fluctuation of LIMO(S) might reflect the difficulties in MO prediction of the active orbitals. In the case of LIMO for SR wave function,[18] DE oscillates around 10212– 10211 hartrees. For the core orbitals, LIMO(S) is expected to predict effective initial guesses in such accuracies, that is, DE of 10212–10211 hartrees, as the prediction scheme is the same as LIMO for SR wave function. As for the active orbitals, however, LIMO(S) as well as LIMO(C) cannot consider the changes in the occupation numbers, namely, the effects of configuration interaction procedure, which is specific in MC wave function theories. The neglect of the change in the occupation numbers would result in the lower accuracy of predicted initial MO guesses. In the case of LIMO(C), the pSCF procedure mixes the active and core orbitals. Therefore, the core orbitals also suffered from the defects in the active orbitals. As the core orbitals contributed to the electronic energy dominantly, LIMO(C) showed much larger DE values than LIMO(S), while the unphysical fluctuations were suppressed. Figures 4a and 4b show the changes in NSCF and DE, respectively, over time for the AIMD simulations of benzene. The
Figure 4. Changes in (a) the number of SCF iterations NSCF and (b) the energy difference DE over time for benzene. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 1. Mean values of the number of SCF iterations NSCF and the energy differences DE (in hartrees) for the AIMD simulations of HCNH1 and benzene. HCNH1
PREMO LIMO(C) LIMO(S)
Benzene
NSCF
DE
NSCF
DE
5.0 3.0 (61%) 2.6 (53%)
6.7 3 1026 1.5 3 1027 5.1 3 1028
9.1 7.9 (87%) 4.5 (50%)
4.0 3 1025 6.9 3 1026 6.0 3 1027
dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods. For PREMO, NSCF oscillated between 8 and 10 cycles, whereas DE ranged from 2.7 3 1025 to 5.3 3 1025 hartrees. In the case of LIMO(C), DE as well as NSCF were relatively lower than the PREMO results, with observed values of 5.1 3 1026–9.5 3 1026 hartrees and 7–9 cycles, respectively. For the LIMO(S) method, DE was observed to be even lower at 6.0 3 10210–1.3 3 1026 hartrees, with an intense unphysical oscillation observed at about 15–25 fs, as in the case of HCNH1. The NSCF values also decreased to 2–6 cycles, reflecting the reduction of DE. Table 1 summarizes the performances of the PREMO, LIMO(C), and LIMO(S) methods discussed in Figures 3 and 4: namely, the mean values of the number of SCF iterations and the energy differences, NSCF and DE, are shown. Approximately two orders of magnitude of reduction in DE was achieved by LIMO(S) compared to PREMO, while DE decreased by one order of magnitude in the LIMO(C) method. Reflecting the reduction in DE, about 50% of reduction in NSCF was achieved by the LIMO(S) method, while NSCF decreased by 15–40% in the case of LIMO(C). Consequently, LIMO(S) showed a superior performance to LIMO(C), because of the avoidance of the degrade in the core orbital space, as noted above.
Numerical tests for prediction degrees of Lagrange polynomial and time steps The dependencies of the performance of LIMO(S) on the prediction degrees of the Lagrange polynomial and the time steps were examined. The prediction degrees were fixed to k 5 0–5, while the dynamic optimization procedures in eq. (5) were used with maximum degrees kmax of three and five. By definition, the simulation with the Lagrange polynomial degree k 5 0 corresponded to PREMO. Fixed time steps (Dt) of 0.1 and 0.25 fs were used. The AIMD simulations were performed for 300 steps after 50 steps of equilibration. The used systems were HCNH1 and benzene, with the initial conditions and the computational levels for ab initio calculations, which were the same as the simulations shown in Figures 3 and 4, respectively. Table 2 summarizes mean values of the number of SCF iterations (NSCF ), averaged over 300 steps. For the simulations with the fixed Lagrange polynomial degrees ranging from k 5 0 to 5, k 5 2 produced the fewest number of iterations for all systems, though the reason for this is unknown. However, utilization of a high prediction degree of Lagrange polynomial such as k 5 5 would be expected to produce worse results than PREMO. For example, the NSCF value for the simulation of
HCNH1 with a time step of 0.1 fs was 7.5 cycles at k 5 5, whereas it was 5.0 cycles for PREMO where k is defined as 0. For the simulations with dynamically optimized prediction degrees using eq. (5), a dependency on the maximum prediction degree kmax was observed, as the observed results for kmax 5 3 and 5 were different from each other. This meant that eq. (5) did not work completely. However, NSCF values were less than PREMO in both cases of kmax 5 3 and 5. In the case of kmax 5 3, NSCF values were comparable to the simulations with the fixed prediction degree of k 5 2, which gave the best performance in the present assessments.
Numerical tests for various systems The performances of LIMO(S) were examined for the various systems and electronic states, including the excited states. The used systems were HCNH1 in the singlet ground (S0) and the singlet first excited (S1) states, nitroxide radical (NOHH) in the doublet ground (D0), and the doublet first excited (D1) states, benzene in the S0 and the S1 states, azobenzene in the S1 state, 4-dimethylamino benzonitrile (DMABN) in the S1 state, and ozone in the S0 and the singlet fourth excited (S4) states. HCNH1 in the excited state is known to be a precursor of interstellar matters.[22] NOHH is the simplest analogue of the nitroxide radicals that are stable under an air atmosphere at the room temperature. Benzene is one of the most simple p conjugated systems, which have doubly degenerated HOMOs and LUMOs. Photo-isomerization of azobenzene[23] attracts interests for the importance in photochromic materials. DMABN, in which the electron-donating dimethylamino group and the electron-withdrawing cyano group connect through a benzene ring, has been extensively studied for its dual fluorescence due to low lying p-p* type locally excited state and intramolecular charge transfer state. Photochemical reaction of Table 2. Dependencies of mean values of the number of SCF iterations NSCF for the AIMD simulations of HCNH1 and benzene on the time steps and the Lagrange prediction degrees. HCNH1
System Time step/fs Prediction degree
Fixed
Auto
0 1 2 3 4 5 3 5
Benzene
0.10
0.25
0.10
0.25
5.0 2.9 2.5 3.8 5.6 7.5 2.6 3.5
6.1 4.3 3.8 3.9 5.6 7.6 3.8 4.5
9.1 4.5 4.4 6.9 10.5 14.0 4.5 6.4
11.3 6.5 6.1 6.9 10.4 14.1 6.2 8.4
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
5
FULL PAPER
WWW.C-CHEM.ORG
Table 3. Mean values of the number of SCF iterations (NSCF ) and CPU times (in s) for SCF (TSCF ), force evaluation (Tforce ), propagations of nucleus and MO (Tprop ), and total calculations (Ttotal ), for the simulations of various systems. PREMO System 1
HCNH HCNH1 NOHH NOHH Benzene Benzene Azobenzene DMABN Ozone Ozone
LIMO(S)
State
Active space
NSCF
TSCF
Tforce
Tprop
Ttotal
S0 S1 D0 D1 S0 S1 S1 S1 S0 S4
10e/10o 10e/10o 13e/10o 13e/10o 6e/6o 6e/6o 6e/4o 6e/5o 16e/15o 12e/9o
5.0 13.3 6.4 9.9 9.1 16.1 16.0 22.9 15.0 12.3
31.8 108.1 38.7 75.1 46.5 76.4 467.4 410.5 65.5 34.7
0.8 0.8 0.8 0.8 16.1 16.0 154.2 106.1 1.7 0.5
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
32.6 108.9 39.5 75.9 62.5 92.4 621.6 516.6 67.6 35.2
ozone plays an essential role in the problem of ozone depletion. The computational levels were CASSCF with the active spaces of 10e/10o, 13e/10o, 6e/6o, 6e/4o, and 6e/5o, for HCNH1, NOHH, benzene, azobenzene, and DMABN, respectively. For ozone molecule, CASSCF calculation with the active space of 12e/9o was performed for the S4 excited state, while the RASSCF method with the active space of 16e/15o was applied in the case of the S0 ground state. The active space was divided into seven, three, and five orbitals in the RAS-1, RAS-2, and RAS-3 spaces, respectively. The electrons up to two were excited from the RAS-1 and/or to RAS-3 spaces. The CASSCF/RASSCF method was used for the systems in the ground states, whereas the state-averaged CASSCF (SACASSCF) method was used for the systems in the excited states. In the case of the excited-state calculations, two states (S0/D0 and S1/D1) were averaged with the same weights for the S1/D1 excited-state calculations, whereas five states (S0, S1, S2, S3, and S4) were averaged with the same weights for the S4 excited state. The 6-31G** Gaussian basis set was used. The degrees of the Lagrange polynomial were determined automatically, using eq. (5) with a maximum prediction degree (kmax) of three. All AIMD simulations were performed under the NVE ensemble with a constant time step of Dt 5 0.1 fs, with slightly distorted initial structures with no velocities. Table 3 summarizes the mean values of the number of SCF iterations and the central processing unit (CPU) times for the AIMD simulations by PREMO and LIMO(S), averaged over 300 steps after 50 steps of equilibration, with the exception of the ozone molecule in the S4 excited state. In this case, photochemical reactions were observed and thus no equilibration procedures were performed to investigate the bond cleavage processes. The AIMD simulations of ozone in the S4 excited state were aborted at the 189th and 212th steps (t 5 18.9 and 21.2 fs) for PREMO and LIMO(S), respectively, due to failures in the SCF convergence. Therefore, the mean values were averaged over 189 and 212 steps for PREMO and LIMO(S), respectively. The CPU times were measured using a single CPU core on an Intel Xeon E5450 3.0 GHz quad-core processor. The mean values of the number of SCF iterations were observed to depend on the system type for the PREMO and LIMO(S) methods. In particular, the excited-state systems with 6
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
NSCF 2.6 4.3 4.7 8.1 4.5 9.7 8.7 9.2 12.3 4.1
(53%) (32%) (73%) (82%) (50%) (60%) (54%) (40%) (82%) (34%)
TSCF 23.4 65.2 31.7 65.9 25.4 47.3 290.3 172.8 57.5 19.8
(74%) (60%) (82%) (88%) (55%) (62%) (62%) (42%) (88%) (57%)
Tforce
Tprop
0.8 0.8 0.8 0.8 16.0 16.0 153.9 106.1 1.7 0.5
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Ttotal 24.2 66.0 32.5 66.7 41.5 63.3 444.3 278.9 59.6 20.3
(74%) (61%) (82%) (88%) (66%) (69%) (71%) (54%) (88%) (58%)
SA wave functions generally required more SCF iterations to achieve convergence compared to the ground-state systems. Despite an observed reduction in NSCF for the LIMO methods in all systems, the performance as indicated by the LIMO/ PREMO ratio ranged from 30 to 80%, which meant that the effectiveness of LIMO depended on the type of system. The mean CPU times for the SCF procedures (TSCF ) were observed to decrease, reflecting the decrease in NSCF , while the LIMO/PREMO ratios for TSCF were about 40–90% and were generally larger than those observed for NSCF (ca. 30–80%). This difference was mainly caused by the extra SCF iteration and the difference in the convergence of iterative diagonalization calculations using Davidson’s algorithm. The extra iteration was needed in each MD step to perform the postprocessing procedures, including the population analysis, while Davidson’s algorithm was used to perform configuration interaction calculations in each MCSCF iteration for active space that included more than six active orbitals. The mean CPU times for force evaluation calculations of PREMO and LIMO (Tforce ) were intrinsically identical to each other. The mean CPU times for propagation of the positions and velocities of nuclei and MO initial guess predictions (Tprop ) were negligibly small for all calculations. Consequently, total CPU times per an AIMD step were reduced by about 10–45%. A photochemical reaction was observed in the simulations of ozone in the S4 excited state. Figure 5 shows the changes of the three OAO bond lengths, the potential energies of the five lowest electronic states, the number of SCF iterations NSCF, and the energy difference DE, over time. In Figure 5a, O1AO2, O1AO3, and O2AO3 are described by dotted, dashed, and solid lines, respectively. At the initial time (t 5 0.0 fs), the O1AO2, ˚, O1AO3, and O2AO3 bond lengths were 1.26, 1.47, and 2.36 A respectively. Based on the reported covalent bond radius of ˚ for oxygen atoms (or bond lengths of 1.28–1.36 0.64–0.68 A [29] ˚ A), the O1AO2 and O1AO3 bonds corresponded to covalent bonds with less than a 10% difference in bond length, while O2 and O3 atoms did not appear to covalently bind to each other. As the AIMD simulation proceeded, both the O1AO3 and O2AO3 bond lengths monotonically increased to 2.74 and ˚ , respectively, and indicated that bond cleavage of the 2.91 A O1AO3 bond occurred in the simulation. In contrast, the ˚ , and O1AO2 bond slightly oscillated between 1.23 and 1.31 A WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 5c shows the change in the energy difference DE, over time. The results for PREMO and LIMO(S) are shown in dotted and solid lines, respectively. As the simulation with PREMO was aborted, the dotted lines are ended at t 5 18.9 fs. For both PREMO and LIMO, DE showed high values of 6.5 3 1022 hartrees at t 5 0.0 fs, because of the absence of initial MO guesses in the first MD step. Except for the spikes around t 5 6.0 and 9.1 fs, DE oscillated around 1026–1024 and 10210– 1027 hartrees for PREMO and LIMO(S), respectively. Around t 5 6.0 and 9.1 fs, the MO predictions completely failed in both PREMO and LIMO(S) cases, with the observed DE values up to about 1022 hartrees. Indeed this was due to the crossings with the higher electronic states, which were not considered to average. However, LIMO(S) as well as PREMO did not show corrupting behaviors in the case of crossings of averaged states, that is, t 5 4.3 and 7.8 fs. The present result showed that LIMO(S) was able to accelerate the SCF convergence in AIMD simulations, involving chemical reactions in the excited state in addition to the ground state.
Conclusions
Figure 5. Changes in (a) OAO bond lengths, (b) potential energies of the five lowest electronic states, and (c) the energy difference DE over time for the AIMD simulations of an ozone molecule in the S4 excited state. In Figure 5a, the dotted, dashed, and solid lines correspond to the O1AO2, O1AO3, and O2AO3 bond lengths, respectively. In Figure 5b, the solid, dashed, dotted, chain, and chain double-dashed lines correspond to the S0, S1, S2, S3, and S4 states, respectively. In Figure 5c, the dotted and solid lines correspond to the results for the PREMO and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
consequently, a photochemical reaction of the ozone molecule was observed. In the photochemical reaction, (pseudo-) crossings of electronic states were observed. Figure 5b shows the potential energies of S0, S1, S2, S3, and S4 states by the solid, dashed, dotted, chain, and chain double-dashed lines, respectively. At the initial time (t 5 0,0 fs), the state energies were 2224.462, 2224.417, 2224.377, 2224.322, and 2224.307 hartrees, for S0–S4 states, respectively. The S0-state energy slightly oscillated with the maximum and the minimum values of 2224.449 and 2224.471 hartrees, respectively, while the S1-state energy monotonically decreased. At t 5 4.3 fs, a crossing between the S3 and S4 states was observed. Also, a crossing between the S2, S3, and S4 states was observed at t 5 7.8 fs. At t 5 6.0 and 9.1 fs, the S4-state energy showed discontinuous behavior, implying the crossings with higher electronic states, such as the S5 and S6 states, which were not considered to average in the present SA-CASSCF calculations. At the dissociated state, that is, t 5 20.0 fs, all states were degenerated with the state energies of 2224.4712224.468 hartrees.
In the present study, we proposed an extension of the LIMO method, which enhances the SCF convergence in AIMD, to the MC wave function theory using the CASSCF and RASSCF methods. While the basic concepts of LIMO could be applied to the MC wave function theories fairly easily, uncertainties existed due to the pSCF effect, that is, the orbital crossing/mixing effect. To take into account the pSCF effect, it is important to separate the occupied and unoccupied orbital spaces. In the case of MC wave functions, the pSCF effects were incorporated by dividing the orbital space into three subspaces consisting of the doubly occupied, active, and unoccupied orbitals. Numerical tests confirmed that the extended LIMO method successfully enhanced the SCF convergence in the AIMD simulations for both groundand excited-state systems with MC wave functions.
Acknowledgments Some of the present calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). Keywords: ab initio molecular dynamics simulation acceleration technique self-consistent field convergence multiconfigurational wave function theory Lagrange interpolation technique
How to cite this article: M. Okoshi, H. Nakai. J. Comput. Chem. 2014, DOI: 10.1002/jcc.23617
[1] B. Kirchner, P. J. di Dio, J. Hutter, Top. Curr. Chem. 2012, 307, 109. [2] Y. Komeiji, Y. Mochizuki, T. Nakano, D. G. Fedorov, J. Mol. Struct. THEOCHEM. 2009, 898, 2. [3] T. Fujiwara, Y. Mochizuki, Y. Komeiji, Y. Okiyama, H. Mori, T. Nakano, E. Miyoshi, Chem. Phys. Lett. 2010, 490, 41.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
7
FULL PAPER
WWW.C-CHEM.ORG
[4] R. P. Steele, M. Head-Gordon, J. C. Tully, J. Phys. Chem. A 2010, 114, 11853. [5] R. Car, M. Parrinello, Phys. Rev. Lett. 1985, 55, 2471. [6] H. B. Schlegel, J. M. Millam, S. S. Iyenger, G. A. Voth, A. D. Daniels, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2001, 114, 9758. [7] S. S. Iyenger, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2001, 115, 10291. [8] H. B. Schlegel, S. S. Iyenger, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2002, 117, 8694. [9] J. L. Alonso, X. Andrade, P. Echenique, F. Falceto, D. Prada-Gracia, A. Rubio, Phys. Rev. Lett. 2008, 101, 096403. [10] J. Jakowski, K. Morokuma, J. Chem. Phys. 2009, 130, 224106. [11] A. M. N. Niklasson, Phys. Rev. Lett. 2008, 100, 123004. [12] G. Zheng, A. M. N. Niklasson, M. Karplus, J. Chem. Phys. 2011, 135, 044122. [13] P. Souvatzis, A. M. N. Niklasson, J. Chem. Phys. 2014, 140, 044117. [14] P. Pulay, J. Comput. Chem. 1982, 3, 556. [15] K. N. Kudin, G. E. Scuseria, E. Cance`s, J. Chem. Phys. 2002, 116, 8255. [16] P. Pulay, G. Fogarasi, Chem. Phys. Lett. 2004, 386, 272. [17] J. M. Herbert, M. Head-Gordon, Phys. Chem. Chem. Phys. 2005, 7, 3269. [18] T. Atsumi, H. Nakai, J. Chem. Phys. 2008, 128, 094101. [19] T. Atsumi, H. Nakai, Chem. Phys. Lett. 2010, 490, 102. [20] A. Sanchez-Galvez, P. Hunt, M. A. Robb, M. Olivucci, T. Vreven, H. B. Schlegel, J. Am. Chem. Soc. 2000, 122, 2911. [21] P. A. Hunt, M. A. Robb, J. Am. Chem. Soc. 2005, 127, 5720. [22] T. Taketsugu, A. Tajima, K. Ishii, T. Hirano, Astrophys. J. 2004, 608, 323. [23] Y. Ootani, K. Satoh, A. Nakaya, T. Noro, T. Taketsugu, J. Chem. Phys. 2009, 131, 194306. [24] A. Nakayama, Y. Harabuchi, S. Yamazaki, T. Taketsugu, Phys. Chem. Chem. Phys. 2013, 15, 12322.
8
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
[25] B. Sellner, M. Barbatti, T. M€ uller, W. Domcke, H. Lischka, Mol. Phys. 2013, 111, 2439. [26] J. Cao, L.-H. Liu, W.-H. Fang, Z.-Z. Xie, Y. Zhang, J. Chem. Phys. 2013, 138, 134306. [27] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, J. A. Pople, Gaussian 03, Revision D.02; Gaussian Inc.: Wallingford, CT, 2004. [28] P. C. Hariharan, J. A. Pople, Theor. Chim. Acta 1973, 28, 213. [29] B. Cordero, V. G omez, A. E. Platero-Prats, M. Rev es, J. Echeverrıa, E. Cremades, F. Barragan, S. Alvarez, Dalton Trans. 2008, 2832.
Received: 25 February 2014 Revised: 7 April 2014 Accepted: 10 April 2014 Published online on 00 Month 2014
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
Acceleration of Self-Consistent Field Convergence in Ab Initio Molecular Dynamics Simulation with Multiconfigurational Wave Function Masaki Okoshi,[a] and Hiromi Nakai*[a,b,c,d] The Lagrange interpolation of molecular orbital (LIMO) method, which reduces the number of self-consistent field iterations in ab initio molecular dynamics simulations with the Hartree–Fock method and the Kohn–Sham density functional theories, is extended to the theory of multiconfigurational wave functions. We examine two types of treatments for the active orbitals that are partially occupied. The first treatment, as denoted by LIMO(C), is a simple application of the conventional LIMO method to the union of the inactive core and the
active orbitals. The second, as denoted by LIMO(S), separately treats the inactive core and the active orbitals. Numerical tests to compare the two treatments clarify that LIMO(S) is superior to LIMO(C). Further applications of LIMO(S) to various systems C 2014 Wiley demonstrate its effectiveness and robustness. V Periodicals, Inc.
Introduction
coefficient matrices of a few previous MD steps. The technique, called LIMO, drastically reduced the SCF iterations, with observed reductions ranging from 50 to 70%. They further extended the LIMO method to geometry optimization and Monte Carlo simulations, which have no information on the simulation (or physical) time.[19] Instead, a least-squares (LS) technique utilizing geometry information was applied to MO propagation (LSMO), and the performance of this technique was comparable to the LIMO technique in AIMD simulations. The single-configurational (SC) theories, including HF and DFT, and the single-reference (SR) theories, such as the Mïller– Plesset perturbation and the coupled-cluster theories, successfully reproduce potential energies and forces acting on atoms in most of ground-state systems. However, especially in quasi-
The ab initio molecular dynamics (AIMD) method, which utilizes ab initio potential energies and forces, has become a useful tool to investigate chemical reaction dynamics.[1] Its major drawback, however, is the high computational cost of the ab initio calculations performed at each molecular dynamics (MD) step. One of the computational bottlenecks in Hartree–Fock (HF) and Kohn–Sham (KS) density functional theory (DFT) calculations is the convergence of the electronic wave function onto the Born–Oppenheimer potential energy surface, which requires the iterative self-consistent field (SCF) procedure. Utilizations of fairly reasonable approximations, such as the linear-scaling technique[2,3] and the dual basis set method,[4] for instance, are capable of reducing computational efforts without the collapse of simulations. Alternatively, several techniques that propagate the electronic structure instead of performing the SCF procedure have been extensively investigated.[5–13] As the SCF convergence in the standard HF/ KS DFT calculations has already been well optimized through methods such as direct inversion in the iterative subspace[14] and its expansions,[15] it might be difficult to achieve further improvements. However, it is possible to obtain information regarding the electronic structure in the AIMD simulation, and this information could potentially be used to further optimize the convergence calculations. Pulay and Fogarasi[16] proposed the Fock matrix dynamics (FMD) method, which extrapolates the Fock matrix forward in time. In their numerical tests, the FMD technique accelerated the SCF convergence by about 50%. Herbert and Head-Gordon[17] further assessed the performances of FMD, in the viewpoints of the efficiencies and accuracies in energy drift. Atsumi and Nakai[18] developed a molecular orbital (MO) propagation technique, which proposes an efficient MO initial guess by the Lagrange interpolation (LI) of MO
DOI: 10.1002/jcc.23617
[a] M. Okoshi, H. Nakai Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan E-mail: [email protected] [b] H. Nakai Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan [c] H. Nakai CREST, Japan Science and Technology Agency, Saitama 332-0012, Japan [d] H. Nakai ESICB, Kyoto University, Kyotodaigaku-Katsura, Kyoto 615-8520, Japan Contract grant sponsor: Core Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST; Theoretical Design of Materials with Innovative Functions Based on Relativistic Electronic Theory); Contract grant sponsor: Elements Strategy Initiative to Form Core Research Center, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan [Elements Strategy Initiative for Catalysts & Batteries (ESICB)]; Contract grant sponsor: Computational Materials Science Initiative (CMSI), MEXT, Japan [Strategic Programs for Innovative Research (SPIRE)]; Contract grant sponsor: Grants for Excellent Graduate Schools, MEXT, Japan (Practical Chemical Wisdom). C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
1
FULL PAPER
WWW.C-CHEM.ORG
degenerate systems, multiconfigurational (MC) or multireference (MR) type wave function theories should be utilized: for instance, an ozone molecule, multibond formation/cleavage systems, and nonadiabatic regions, in which two or more electronic states are (pseudo-) crossing each other. While the MC/ MR wave function theories generally need much more computational efforts than the SC/SR wave function theories, AIMD with the MC/MR wave function theories have been extensively performed, particularly for the systems, in which photochemical reactions with nonradiative transitions take place.[20–26] In the present article, we extended the LIMO method, which was originally developed for the SC wave function, to the MC wave function theories, that is, the complete active space SCF (CASSCF) and restricted active space SCF (RASSCF) methods. Theoretical considerations were initially made and numerical tests were performed.
Theoretical Aspects LIMO for SC wave functions This subsection briefly reviews the original LIMO method developed for the SC wave function.[18] In AIMD, the conðnÞ verged MO coefficient matrix of the last (nth) step, Cconv , is often used as the initial guess for the (n 1 1)th step because the geometric and electronic structures of the neighboring MD steps are expected to be similar to each other. We refer to this conventional technique as PREMO, to refer to the fact that only the MO coefficient matrix of the previous step is used for the initial MO guess. The LIMO technique predicts better initial MO guesses than PREMO, by using the converged MO coefficient matrices in the last step as well as in several previous steps. The main concept of LIMO is that the MO coefficient matrices are expected to change smoothly because the orbital centers move continuously in the AIMD simulations. Thus, an effective initial MO guess can be predicted as the linear combination of the several converged MO coefficient matrices. The predicted initial ðn11Þ MO guess for the (n 1 1)th step, Cprd , is given by the LI technique:
ðn21Þ!n
TpSCF
21 t ðn21Þ n ðn21Þ ðn21Þ 5 t Cconv Cconv Cconv Cconv
(3)
and ðn2lÞ!n
TpSCF
5
l21 Y
ðn2l1mÞ!ðn2l1m11Þ
TpSCF
(4)
m50
where the subscript “pSCF” refers to the pseudo-SCF effect that is present in the (n-l)th and nth converged MO coefficient matrices. In the case of the PREMO method, the transformaðn21Þ!n tion matrix for the neighboring steps, TpSCF , represents the SCF effect in the nth MD step. Note that the pSCF effect is considered for the occupied orbitals only; that is, the transformation matrix T has the dimensions of Nocc 3 Nocc, and the MO coefficient matrix C and its transpose tC have the dimensions of Nbasis 3 Nocc and Nocc 3 Nbasis, respectively, where Nbasis and Nocc are the numbers of the basis functions and of the occupied orbitals. Mathematical proofs have shown that the orbital-invariant LIMO method for all occupied and unoccupied orbitals is completely meaningless, since the effective ðn11Þ MO initial guess, Cprd , becomes identical to the converged [18] MO coefficient matrix of the last step, CðnÞ conv . The effectiveness of the MO initial guess predicted by LIMO depends on the degree of the Lagrange polynomial (k). Numerical assessments have shown the existence of the optimal degree kopt, which depends mainly on the time step of AIMD, Dt.[18] While it is unfeasible to determine kopt prior to performing SCF, applying the optimal degree at the nth step to the (n 1 1)th step results in acceptable performance. Therefore, the optimal degree kopt is determined to minimize the following deviation sum: k X X X ðn2lÞ!n n ðn2lÞ dk 5 Lk;l Cconv ðl; jÞTpSCF ðj; iÞ Cconv ðl; iÞ2 l;i j l50
(5)
where l runs over the basis functions, whereas i and j run over the occupied MOs.
Extension of LIMO to the MC wave function ðn11Þ
Cprd 5
k X
ðn2lÞ!n
Lk;l Cðn2lÞ conv TpSCF
(1)
l50
where fCðn2lÞ conv g (l 5 0, 1, 2, . . ., k) are the converged MO coefficient matrices at the (n-l)th step, and fLk;l g are the lth Lagrange coefficients for the kth degree of polynomial given by Lk;l 5
k Y t ðn11Þ 2tðn2mÞ m50
tðn2lÞ 2t ðn2mÞ
(2)
When adopting a constant time step in AIMD, the Lagrange coefficients become simple integers. The transformation ðn2lÞ!n matrix, fTpSCF g, leads to the orbital-invariant formula, by considering the orbital crossing and mixing effects within the ðn2lÞ!n occupied orbitals. The explicit formula of TpSCF is given by: 2
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
As described above, pSCF can only be considered for the occupied orbitals; therefore, it is indispensable to separate the occupied and the unoccupied orbital spaces. In the case of the SC wave function, the occupied/unoccupied orbitals are well defined. However, in the MC wave function case, there remains ambiguity to define the occupied orbitals, since the active orbitals are partially occupied. Therefore, the main problem is how to treat the pSCF effect. In the present study, we examined two types of effective MO initial guesses of dimensions Nbasis 3 Nocc, where Nocc refers to the summation of the number of the inactive core and the active orbitals. One method to treat the pSCF effect is the conventional LIMO scheme for the occupied orbitals, that is, the union of the inactive core and the active orbitals. In this case, all occupied orbitals are treated using eqs. (1)–(5) in the similar procedure as the SC case, as if they are doubly occupied orbitals. WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 1. Initial geometries of (a) HCNH1 and (b) benzene.
Another method is to treat the inactive core and the active orbitals separately and independently: ðn11Þ
Cprd; core 5
k X
ðn2lÞ!n
Lk;l Cðn2lÞ conv; core TpSCF; core
(6)
l50 ðn11Þ
Cprd; active 5
k X
ðn2lÞ
ðn2lÞ!n
Lk;l Cconv; active TpSCF; active
(7)
l50
where the subscripts “core” and “active” refer to the inactive core and the active orbitals, respectively. In this method, the matrices Ccore, Tcore, Cactive, and Tactive have the dimensions of Nbasis 3 Ncore, Ncore 3 Ncore, Nbasis 3 Nactive, and Nactive 3 Nactive, respectively, where Ncore and Nactive refer to the numbers of the inactive core and the active orbitals. The total predicted MO coefficient matrix is formed as the direct sum of Ccore and Cactive: ðn11Þ
ðn11Þ
ðn11Þ
Cprd 5Cprd; core 丣Cprd; active
(8)
ðn11Þ
FULL PAPER
on the atoms were evaluated using the CASSCF method with Pople’s double-zeta-plus-polarization type basis set, 6-31G**.[28] The used active spaces were 10 electrons within 10 orbitals and six electrons within six orbitals as denoted by 10e/10o and 6e/6o for HCNH1 and benzene, respectively. The active space included full valence electrons/orbitals in the case of HCNH1, while all p electrons/orbitals were used as the active space for benzene. The degrees of the Lagrange polynomial k were determined automatically using eq. (5), with a maximum prediction degree kmax of three. The Velocity Verlet method with a time step (Dt) of 0.1 fs was adopted for the numerical integration of the Newtonian equation of motion of nuclei under the microcanonical (NVE) ensemble. The AIMD simulations were performed for 300 steps (30.0 fs) after 50 steps (5.0 fs) of equilibration, for each system. Randomly strained structures from the equilibrium with no velocities, which are shown in Figures 1a and 1b for HCNH1 and benzene, respectively, were used as the initial conditions. The distortion energies for HCNH1 and benzene were 1.05 and 7.44 eV, whereas their zero-point vibrational energies (ZVEs) were 0.75 and 2.92 eV, respectively. ZVEs were calculated by summing up the normal mode frequencies, calculated for the equilibrium structures at the same levels of theory as the AIMD simulations; that is, CASSCF(10e/10o)/6-31G** and CASSCF(6e/6o)/6-31G** for HCNH1 and benzene, respectively. Figures 2a and 2b show the energy variations in the AIMD simulations of HCNH1 and benzene, respectively. The dotted, dashed, and solid lines correspond to the results for the kinetic, potential, and total energies of the system. The potential and total energies are relative values to those at t 5 25.0 fs, that is, the initial conditions. The simulations at t 5 25.0–0.0 fs correspond to the equilibration, whereas those
The initial MO guess, Cprd , is orthonormalized prior to the MCSCF procedure through Schmidt’s procedure in all cases, including PREMO, the conventional LIMO method, and the LIMO method for the separated orbitals. For convenience, the conventional scheme is denoted by LIMO(C) in this article, whereas the separated procedure given in eq. (8) is denoted by LIMO(S). The optimal degree of the Lagrange polynomial is determined in a manner similar to the conventional LIMO through use of eq. (5) for the total MO coefficient matrices.
Numerical Assessments Numerical performances of the LIMO technique were investigated further. AIMD simulations, as well as MO prediction procedures, are typically performed using an original MD code combined with quantum chemistry program packages that use Gaussian basis functions. In the present AIMD simulations, the Gaussian 03 suite of programs[27] was utilized to perform ab initio MCSCF calculations. Numerical tests for prediction methods The performances of the LIMO(C) and LIMO(S) prediction schemes were compared by AIMD simulations for HCNH1 and benzene molecules in the S0 ground state. The forces acting
Figure 2. Energy variations in the AIMD simulations of (a) HCNH1 and (b) benzene. The dotted, dashed, and solid lines correspond to the kinetic, potential, and total energies, respectively.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 3. Changes in (a) the number of SCF iterations NSCF and (b) the energy difference DE over time for HCNH1. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
at t 5 0.0–30.0 were considered as the production runs. For HCNH1 (Fig. 2a), the kinetic energy (dotted line) rapidly increased to about 1.0 eV at 25.0–0.0 fs, whereas the potential energy (dashed line) decreased as much. This meant that a part of distortion energy changed into the kinetic energy with the geometric relaxation. After that, the kinetic energy as well as the potential energy oscillated with a period of about 7.0– 9.0 fs. The total energy (solid line) was approximately zero, with the maximum absolute value of 0.53 meV. As for benzene (Fig. 2b), the similar tendency was observed with the large kinetic energy of about 5.0 eV. The maximum absolute value for the total energy was 2.8 meV. Figures 3a and 3b show the number of SCF iterations (NSCF) and the energy difference DE, respectively, over time for the AIMD simulations of HCNH1. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. The energy difference DE was used to assess the performance of each prediction method: ðnÞ
ðnÞ DE5jE½Cprd; k 2E½Cconv j
(9)
where the two terms on the right-hand side of the equation indicate the potential energies evaluated using the MO coeffiðnÞ cient matrices of Cprd; k and CðnÞ conv , respectively. For the PREMO method (dotted line), NSCF slightly oscillated between 4 and 6 cycles, whereas DE ranged from 1.3 3 1026 to 1.4 3 1025 hartrees. In the case of LIMO(C) (dashed line), DE was observed to be lower compared to the results from the PREMO method, with DE ranging from 2.5 3 1029 to 4.5 3 1027 hartrees. The values for NSCF were also relatively low, ranging from 2 to 4 cycles. In both PREMO and LIMO(C), oscil4
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
lations of DE with a period of 7.0–9.0 fs were observed. These correspond to the oscillation of the kinetic energy of the system, as shown in Figure 2a. On the contrary, unphysical fluctuations with high amplitudes, ranging from 2.9 3 10211 to 2.7 3 1027 hartrees, were observed in the case of LIMO(S). Despite the existence of the unphysical fluctuations, DE for LIMO(S) was generally lower than that of LIMO(C) with exceptions at several timings of t at 3.4, 3.9, 10.8, 25.4, 25.5, 25.9–26.3, 26.5–26.8, 27.4, and 27.5 fs, and the NSCF value was comparable to that of LIMO(C) at 2–4 cycles. The unphysical fluctuation of LIMO(S) might reflect the difficulties in MO prediction of the active orbitals. In the case of LIMO for SR wave function,[18] DE oscillates around 10212– 10211 hartrees. For the core orbitals, LIMO(S) is expected to predict effective initial guesses in such accuracies, that is, DE of 10212–10211 hartrees, as the prediction scheme is the same as LIMO for SR wave function. As for the active orbitals, however, LIMO(S) as well as LIMO(C) cannot consider the changes in the occupation numbers, namely, the effects of configuration interaction procedure, which is specific in MC wave function theories. The neglect of the change in the occupation numbers would result in the lower accuracy of predicted initial MO guesses. In the case of LIMO(C), the pSCF procedure mixes the active and core orbitals. Therefore, the core orbitals also suffered from the defects in the active orbitals. As the core orbitals contributed to the electronic energy dominantly, LIMO(C) showed much larger DE values than LIMO(S), while the unphysical fluctuations were suppressed. Figures 4a and 4b show the changes in NSCF and DE, respectively, over time for the AIMD simulations of benzene. The
Figure 4. Changes in (a) the number of SCF iterations NSCF and (b) the energy difference DE over time for benzene. The dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Table 1. Mean values of the number of SCF iterations NSCF and the energy differences DE (in hartrees) for the AIMD simulations of HCNH1 and benzene. HCNH1
PREMO LIMO(C) LIMO(S)
Benzene
NSCF
DE
NSCF
DE
5.0 3.0 (61%) 2.6 (53%)
6.7 3 1026 1.5 3 1027 5.1 3 1028
9.1 7.9 (87%) 4.5 (50%)
4.0 3 1025 6.9 3 1026 6.0 3 1027
dotted, dashed, and solid lines correspond to the results for the PREMO, LIMO(C), and LIMO(S) methods. For PREMO, NSCF oscillated between 8 and 10 cycles, whereas DE ranged from 2.7 3 1025 to 5.3 3 1025 hartrees. In the case of LIMO(C), DE as well as NSCF were relatively lower than the PREMO results, with observed values of 5.1 3 1026–9.5 3 1026 hartrees and 7–9 cycles, respectively. For the LIMO(S) method, DE was observed to be even lower at 6.0 3 10210–1.3 3 1026 hartrees, with an intense unphysical oscillation observed at about 15–25 fs, as in the case of HCNH1. The NSCF values also decreased to 2–6 cycles, reflecting the reduction of DE. Table 1 summarizes the performances of the PREMO, LIMO(C), and LIMO(S) methods discussed in Figures 3 and 4: namely, the mean values of the number of SCF iterations and the energy differences, NSCF and DE, are shown. Approximately two orders of magnitude of reduction in DE was achieved by LIMO(S) compared to PREMO, while DE decreased by one order of magnitude in the LIMO(C) method. Reflecting the reduction in DE, about 50% of reduction in NSCF was achieved by the LIMO(S) method, while NSCF decreased by 15–40% in the case of LIMO(C). Consequently, LIMO(S) showed a superior performance to LIMO(C), because of the avoidance of the degrade in the core orbital space, as noted above.
Numerical tests for prediction degrees of Lagrange polynomial and time steps The dependencies of the performance of LIMO(S) on the prediction degrees of the Lagrange polynomial and the time steps were examined. The prediction degrees were fixed to k 5 0–5, while the dynamic optimization procedures in eq. (5) were used with maximum degrees kmax of three and five. By definition, the simulation with the Lagrange polynomial degree k 5 0 corresponded to PREMO. Fixed time steps (Dt) of 0.1 and 0.25 fs were used. The AIMD simulations were performed for 300 steps after 50 steps of equilibration. The used systems were HCNH1 and benzene, with the initial conditions and the computational levels for ab initio calculations, which were the same as the simulations shown in Figures 3 and 4, respectively. Table 2 summarizes mean values of the number of SCF iterations (NSCF ), averaged over 300 steps. For the simulations with the fixed Lagrange polynomial degrees ranging from k 5 0 to 5, k 5 2 produced the fewest number of iterations for all systems, though the reason for this is unknown. However, utilization of a high prediction degree of Lagrange polynomial such as k 5 5 would be expected to produce worse results than PREMO. For example, the NSCF value for the simulation of
HCNH1 with a time step of 0.1 fs was 7.5 cycles at k 5 5, whereas it was 5.0 cycles for PREMO where k is defined as 0. For the simulations with dynamically optimized prediction degrees using eq. (5), a dependency on the maximum prediction degree kmax was observed, as the observed results for kmax 5 3 and 5 were different from each other. This meant that eq. (5) did not work completely. However, NSCF values were less than PREMO in both cases of kmax 5 3 and 5. In the case of kmax 5 3, NSCF values were comparable to the simulations with the fixed prediction degree of k 5 2, which gave the best performance in the present assessments.
Numerical tests for various systems The performances of LIMO(S) were examined for the various systems and electronic states, including the excited states. The used systems were HCNH1 in the singlet ground (S0) and the singlet first excited (S1) states, nitroxide radical (NOHH) in the doublet ground (D0), and the doublet first excited (D1) states, benzene in the S0 and the S1 states, azobenzene in the S1 state, 4-dimethylamino benzonitrile (DMABN) in the S1 state, and ozone in the S0 and the singlet fourth excited (S4) states. HCNH1 in the excited state is known to be a precursor of interstellar matters.[22] NOHH is the simplest analogue of the nitroxide radicals that are stable under an air atmosphere at the room temperature. Benzene is one of the most simple p conjugated systems, which have doubly degenerated HOMOs and LUMOs. Photo-isomerization of azobenzene[23] attracts interests for the importance in photochromic materials. DMABN, in which the electron-donating dimethylamino group and the electron-withdrawing cyano group connect through a benzene ring, has been extensively studied for its dual fluorescence due to low lying p-p* type locally excited state and intramolecular charge transfer state. Photochemical reaction of Table 2. Dependencies of mean values of the number of SCF iterations NSCF for the AIMD simulations of HCNH1 and benzene on the time steps and the Lagrange prediction degrees. HCNH1
System Time step/fs Prediction degree
Fixed
Auto
0 1 2 3 4 5 3 5
Benzene
0.10
0.25
0.10
0.25
5.0 2.9 2.5 3.8 5.6 7.5 2.6 3.5
6.1 4.3 3.8 3.9 5.6 7.6 3.8 4.5
9.1 4.5 4.4 6.9 10.5 14.0 4.5 6.4
11.3 6.5 6.1 6.9 10.4 14.1 6.2 8.4
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
5
FULL PAPER
WWW.C-CHEM.ORG
Table 3. Mean values of the number of SCF iterations (NSCF ) and CPU times (in s) for SCF (TSCF ), force evaluation (Tforce ), propagations of nucleus and MO (Tprop ), and total calculations (Ttotal ), for the simulations of various systems. PREMO System 1
HCNH HCNH1 NOHH NOHH Benzene Benzene Azobenzene DMABN Ozone Ozone
LIMO(S)
State
Active space
NSCF
TSCF
Tforce
Tprop
Ttotal
S0 S1 D0 D1 S0 S1 S1 S1 S0 S4
10e/10o 10e/10o 13e/10o 13e/10o 6e/6o 6e/6o 6e/4o 6e/5o 16e/15o 12e/9o
5.0 13.3 6.4 9.9 9.1 16.1 16.0 22.9 15.0 12.3
31.8 108.1 38.7 75.1 46.5 76.4 467.4 410.5 65.5 34.7
0.8 0.8 0.8 0.8 16.1 16.0 154.2 106.1 1.7 0.5
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
32.6 108.9 39.5 75.9 62.5 92.4 621.6 516.6 67.6 35.2
ozone plays an essential role in the problem of ozone depletion. The computational levels were CASSCF with the active spaces of 10e/10o, 13e/10o, 6e/6o, 6e/4o, and 6e/5o, for HCNH1, NOHH, benzene, azobenzene, and DMABN, respectively. For ozone molecule, CASSCF calculation with the active space of 12e/9o was performed for the S4 excited state, while the RASSCF method with the active space of 16e/15o was applied in the case of the S0 ground state. The active space was divided into seven, three, and five orbitals in the RAS-1, RAS-2, and RAS-3 spaces, respectively. The electrons up to two were excited from the RAS-1 and/or to RAS-3 spaces. The CASSCF/RASSCF method was used for the systems in the ground states, whereas the state-averaged CASSCF (SACASSCF) method was used for the systems in the excited states. In the case of the excited-state calculations, two states (S0/D0 and S1/D1) were averaged with the same weights for the S1/D1 excited-state calculations, whereas five states (S0, S1, S2, S3, and S4) were averaged with the same weights for the S4 excited state. The 6-31G** Gaussian basis set was used. The degrees of the Lagrange polynomial were determined automatically, using eq. (5) with a maximum prediction degree (kmax) of three. All AIMD simulations were performed under the NVE ensemble with a constant time step of Dt 5 0.1 fs, with slightly distorted initial structures with no velocities. Table 3 summarizes the mean values of the number of SCF iterations and the central processing unit (CPU) times for the AIMD simulations by PREMO and LIMO(S), averaged over 300 steps after 50 steps of equilibration, with the exception of the ozone molecule in the S4 excited state. In this case, photochemical reactions were observed and thus no equilibration procedures were performed to investigate the bond cleavage processes. The AIMD simulations of ozone in the S4 excited state were aborted at the 189th and 212th steps (t 5 18.9 and 21.2 fs) for PREMO and LIMO(S), respectively, due to failures in the SCF convergence. Therefore, the mean values were averaged over 189 and 212 steps for PREMO and LIMO(S), respectively. The CPU times were measured using a single CPU core on an Intel Xeon E5450 3.0 GHz quad-core processor. The mean values of the number of SCF iterations were observed to depend on the system type for the PREMO and LIMO(S) methods. In particular, the excited-state systems with 6
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
NSCF 2.6 4.3 4.7 8.1 4.5 9.7 8.7 9.2 12.3 4.1
(53%) (32%) (73%) (82%) (50%) (60%) (54%) (40%) (82%) (34%)
TSCF 23.4 65.2 31.7 65.9 25.4 47.3 290.3 172.8 57.5 19.8
(74%) (60%) (82%) (88%) (55%) (62%) (62%) (42%) (88%) (57%)
Tforce
Tprop
0.8 0.8 0.8 0.8 16.0 16.0 153.9 106.1 1.7 0.5
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Ttotal 24.2 66.0 32.5 66.7 41.5 63.3 444.3 278.9 59.6 20.3
(74%) (61%) (82%) (88%) (66%) (69%) (71%) (54%) (88%) (58%)
SA wave functions generally required more SCF iterations to achieve convergence compared to the ground-state systems. Despite an observed reduction in NSCF for the LIMO methods in all systems, the performance as indicated by the LIMO/ PREMO ratio ranged from 30 to 80%, which meant that the effectiveness of LIMO depended on the type of system. The mean CPU times for the SCF procedures (TSCF ) were observed to decrease, reflecting the decrease in NSCF , while the LIMO/PREMO ratios for TSCF were about 40–90% and were generally larger than those observed for NSCF (ca. 30–80%). This difference was mainly caused by the extra SCF iteration and the difference in the convergence of iterative diagonalization calculations using Davidson’s algorithm. The extra iteration was needed in each MD step to perform the postprocessing procedures, including the population analysis, while Davidson’s algorithm was used to perform configuration interaction calculations in each MCSCF iteration for active space that included more than six active orbitals. The mean CPU times for force evaluation calculations of PREMO and LIMO (Tforce ) were intrinsically identical to each other. The mean CPU times for propagation of the positions and velocities of nuclei and MO initial guess predictions (Tprop ) were negligibly small for all calculations. Consequently, total CPU times per an AIMD step were reduced by about 10–45%. A photochemical reaction was observed in the simulations of ozone in the S4 excited state. Figure 5 shows the changes of the three OAO bond lengths, the potential energies of the five lowest electronic states, the number of SCF iterations NSCF, and the energy difference DE, over time. In Figure 5a, O1AO2, O1AO3, and O2AO3 are described by dotted, dashed, and solid lines, respectively. At the initial time (t 5 0.0 fs), the O1AO2, ˚, O1AO3, and O2AO3 bond lengths were 1.26, 1.47, and 2.36 A respectively. Based on the reported covalent bond radius of ˚ for oxygen atoms (or bond lengths of 1.28–1.36 0.64–0.68 A [29] ˚ A), the O1AO2 and O1AO3 bonds corresponded to covalent bonds with less than a 10% difference in bond length, while O2 and O3 atoms did not appear to covalently bind to each other. As the AIMD simulation proceeded, both the O1AO3 and O2AO3 bond lengths monotonically increased to 2.74 and ˚ , respectively, and indicated that bond cleavage of the 2.91 A O1AO3 bond occurred in the simulation. In contrast, the ˚ , and O1AO2 bond slightly oscillated between 1.23 and 1.31 A WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 5c shows the change in the energy difference DE, over time. The results for PREMO and LIMO(S) are shown in dotted and solid lines, respectively. As the simulation with PREMO was aborted, the dotted lines are ended at t 5 18.9 fs. For both PREMO and LIMO, DE showed high values of 6.5 3 1022 hartrees at t 5 0.0 fs, because of the absence of initial MO guesses in the first MD step. Except for the spikes around t 5 6.0 and 9.1 fs, DE oscillated around 1026–1024 and 10210– 1027 hartrees for PREMO and LIMO(S), respectively. Around t 5 6.0 and 9.1 fs, the MO predictions completely failed in both PREMO and LIMO(S) cases, with the observed DE values up to about 1022 hartrees. Indeed this was due to the crossings with the higher electronic states, which were not considered to average. However, LIMO(S) as well as PREMO did not show corrupting behaviors in the case of crossings of averaged states, that is, t 5 4.3 and 7.8 fs. The present result showed that LIMO(S) was able to accelerate the SCF convergence in AIMD simulations, involving chemical reactions in the excited state in addition to the ground state.
Conclusions
Figure 5. Changes in (a) OAO bond lengths, (b) potential energies of the five lowest electronic states, and (c) the energy difference DE over time for the AIMD simulations of an ozone molecule in the S4 excited state. In Figure 5a, the dotted, dashed, and solid lines correspond to the O1AO2, O1AO3, and O2AO3 bond lengths, respectively. In Figure 5b, the solid, dashed, dotted, chain, and chain double-dashed lines correspond to the S0, S1, S2, S3, and S4 states, respectively. In Figure 5c, the dotted and solid lines correspond to the results for the PREMO and LIMO(S) methods, respectively. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
consequently, a photochemical reaction of the ozone molecule was observed. In the photochemical reaction, (pseudo-) crossings of electronic states were observed. Figure 5b shows the potential energies of S0, S1, S2, S3, and S4 states by the solid, dashed, dotted, chain, and chain double-dashed lines, respectively. At the initial time (t 5 0,0 fs), the state energies were 2224.462, 2224.417, 2224.377, 2224.322, and 2224.307 hartrees, for S0–S4 states, respectively. The S0-state energy slightly oscillated with the maximum and the minimum values of 2224.449 and 2224.471 hartrees, respectively, while the S1-state energy monotonically decreased. At t 5 4.3 fs, a crossing between the S3 and S4 states was observed. Also, a crossing between the S2, S3, and S4 states was observed at t 5 7.8 fs. At t 5 6.0 and 9.1 fs, the S4-state energy showed discontinuous behavior, implying the crossings with higher electronic states, such as the S5 and S6 states, which were not considered to average in the present SA-CASSCF calculations. At the dissociated state, that is, t 5 20.0 fs, all states were degenerated with the state energies of 2224.4712224.468 hartrees.
In the present study, we proposed an extension of the LIMO method, which enhances the SCF convergence in AIMD, to the MC wave function theory using the CASSCF and RASSCF methods. While the basic concepts of LIMO could be applied to the MC wave function theories fairly easily, uncertainties existed due to the pSCF effect, that is, the orbital crossing/mixing effect. To take into account the pSCF effect, it is important to separate the occupied and unoccupied orbital spaces. In the case of MC wave functions, the pSCF effects were incorporated by dividing the orbital space into three subspaces consisting of the doubly occupied, active, and unoccupied orbitals. Numerical tests confirmed that the extended LIMO method successfully enhanced the SCF convergence in the AIMD simulations for both groundand excited-state systems with MC wave functions.
Acknowledgments Some of the present calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). Keywords: ab initio molecular dynamics simulation acceleration technique self-consistent field convergence multiconfigurational wave function theory Lagrange interpolation technique
How to cite this article: M. Okoshi, H. Nakai. J. Comput. Chem. 2014, DOI: 10.1002/jcc.23617
[1] B. Kirchner, P. J. di Dio, J. Hutter, Top. Curr. Chem. 2012, 307, 109. [2] Y. Komeiji, Y. Mochizuki, T. Nakano, D. G. Fedorov, J. Mol. Struct. THEOCHEM. 2009, 898, 2. [3] T. Fujiwara, Y. Mochizuki, Y. Komeiji, Y. Okiyama, H. Mori, T. Nakano, E. Miyoshi, Chem. Phys. Lett. 2010, 490, 41.
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
7
FULL PAPER
WWW.C-CHEM.ORG
[4] R. P. Steele, M. Head-Gordon, J. C. Tully, J. Phys. Chem. A 2010, 114, 11853. [5] R. Car, M. Parrinello, Phys. Rev. Lett. 1985, 55, 2471. [6] H. B. Schlegel, J. M. Millam, S. S. Iyenger, G. A. Voth, A. D. Daniels, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2001, 114, 9758. [7] S. S. Iyenger, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2001, 115, 10291. [8] H. B. Schlegel, S. S. Iyenger, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frisch, J. Chem. Phys. 2002, 117, 8694. [9] J. L. Alonso, X. Andrade, P. Echenique, F. Falceto, D. Prada-Gracia, A. Rubio, Phys. Rev. Lett. 2008, 101, 096403. [10] J. Jakowski, K. Morokuma, J. Chem. Phys. 2009, 130, 224106. [11] A. M. N. Niklasson, Phys. Rev. Lett. 2008, 100, 123004. [12] G. Zheng, A. M. N. Niklasson, M. Karplus, J. Chem. Phys. 2011, 135, 044122. [13] P. Souvatzis, A. M. N. Niklasson, J. Chem. Phys. 2014, 140, 044117. [14] P. Pulay, J. Comput. Chem. 1982, 3, 556. [15] K. N. Kudin, G. E. Scuseria, E. Cance`s, J. Chem. Phys. 2002, 116, 8255. [16] P. Pulay, G. Fogarasi, Chem. Phys. Lett. 2004, 386, 272. [17] J. M. Herbert, M. Head-Gordon, Phys. Chem. Chem. Phys. 2005, 7, 3269. [18] T. Atsumi, H. Nakai, J. Chem. Phys. 2008, 128, 094101. [19] T. Atsumi, H. Nakai, Chem. Phys. Lett. 2010, 490, 102. [20] A. Sanchez-Galvez, P. Hunt, M. A. Robb, M. Olivucci, T. Vreven, H. B. Schlegel, J. Am. Chem. Soc. 2000, 122, 2911. [21] P. A. Hunt, M. A. Robb, J. Am. Chem. Soc. 2005, 127, 5720. [22] T. Taketsugu, A. Tajima, K. Ishii, T. Hirano, Astrophys. J. 2004, 608, 323. [23] Y. Ootani, K. Satoh, A. Nakaya, T. Noro, T. Taketsugu, J. Chem. Phys. 2009, 131, 194306. [24] A. Nakayama, Y. Harabuchi, S. Yamazaki, T. Taketsugu, Phys. Chem. Chem. Phys. 2013, 15, 12322.
8
Journal of Computational Chemistry 2014, DOI: 10.1002/jcc.23617
[25] B. Sellner, M. Barbatti, T. M€ uller, W. Domcke, H. Lischka, Mol. Phys. 2013, 111, 2439. [26] J. Cao, L.-H. Liu, W.-H. Fang, Z.-Z. Xie, Y. Zhang, J. Chem. Phys. 2013, 138, 134306. [27] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, J. A. Pople, Gaussian 03, Revision D.02; Gaussian Inc.: Wallingford, CT, 2004. [28] P. C. Hariharan, J. A. Pople, Theor. Chim. Acta 1973, 28, 213. [29] B. Cordero, V. G omez, A. E. Platero-Prats, M. Rev es, J. Echeverrıa, E. Cremades, F. Barragan, S. Alvarez, Dalton Trans. 2008, 2832.
Received: 25 February 2014 Revised: 7 April 2014 Accepted: 10 April 2014 Published online on 00 Month 2014
WWW.CHEMISTRYVIEWS.COM
