J. Chem. Theory Comput. 2009, 5, 287–294
287
Fully Automated Incremental Evaluation of MP2 and CCSD(T) Energies: Application to Water Clusters Joachim Friedrich* and Michael Dolg Institute for Theoretical Chemistry, UniVersity of Cologne, Greinstrasse 4, 50939 Cologne, Germany Received January 1, 2008
Abstract: A fully automated implementation of the incremental scheme for CCSD energies has been extended to treat MP2 and CCSD(T) energies. It is shown in applications on water clusters that the error of the incremental expansion for the total energy is below 1 kcal/mol already at second or third order. It is demonstrated that the approach saves CPU time, RAM, and disk space. Finally it is shown that the calculations can be run in parallel on up to 50 CPUs, without significant loss of computer time.
1. Introduction It is well-known that a quantitative description of the electronic structure in molecules is usually not possible with the Hartree-Fock (HF) method. One way to achieve higher accuracy in electronic structure calculations is to improve the wave function, which can be routinely done by manybody perturbation theory (MBPT), configuration interaction theory (CI), or coupled-cluster theory (CC). The major drawback of these approaches is their strong dependence of the computational effort on the size of the one-particle basis set. This means that these approaches depend heavily on the system size too, if atom-centered basis functions are used. Since for large systems the canonical HF orbitals are not necessarily the best choice for a PT, CI, or CC expansion of the wave function, many groups use a local orbital basis instead to include electron correlation. This allows one to screen out insignificant contributions to the energy, and therefore the computational cost is reduced.1-19 Conceptually different approaches divide the total system into parts and then perform a perturbation expansion to obtain the total correlation energy.20-23 An approach designed in this way is the incremental scheme of Stoll.24-26 It is based on the Bethe-Goldstone expansion, which was introduced to quantum chemistry by Nesbet27-29 more than 40 years ago. The incremental scheme was successfully applied during the past 15 years to various periodic systems30-34 and mole* To whom correspondence should be addressed. Tel.: (+49) (0)221-470-6886. Fax: (+49) (0)221-470-6896. E-mail:
[email protected].
cules.35-39 Besides the treatment of closed-shell systems, extensions to open-shell cases have been developed too.40-42 Recently we proposed a fully automated implementation of the incremental scheme for CCSD energies,36 implemented an automatic distance screening,37 extended the approach for the usage of symmetry39 and to the RCCSD method for open-shell calculations.42 To account for the local character of the core electrons, we introduced an efficient scheme to treat the core and core-valence correlation.43 In this work we now extend our fully automated implementation to second-order Møller-Plesset perturbation theory (MP2) and to the coupled-cluster ansatz with singles, doubles, and perturbative triples excitations CCSD(T). Furthermore, as recently proposed, we use a second basis set to describe the environment of the orbital domains in a computationally efficient manner.44 We note that the idea of multiple basis sets is not entirely new, since Jurgens-Lutovsky and Almlo¨f made use of a reduced basis for the occupied space in MP2 calculations already in 199145 and Klopper et al. made use of a reduced basis for the treatment of the triples in CCSD(T) calculations in 1997.46 However, we exploit this idea in the framework of the incremental scheme and check how it performs with respect to the overall accuracy and the CPU time requirements.
2. Theory 2.1. Incremental Scheme. In an incremental calculation we divide the total system into small domains consisting of groups of localized occupied orbitals according to the procedure outlined in refs 36-38. Then we calculate the
10.1021/ct800355e CCC: $40.75 2009 American Chemical Society Published on Web 01/15/2009
288 J. Chem. Theory Comput., Vol. 5, No. 2, 2009
Friedrich and Dolg
correlation energies for these domains. To include the nonadditivity corrections, we also calculate correction energies of pairs and triples, etc., of domains, until we reach the desired accuracy. The correlation energy is then computed according to Ecorr )
∑ ∆εi + 2!1 ∑ ∆εij + 3!1 ∑ ∆εijk + ... i
ij
∆εi ) εi
(1)
ijk
∆εij ) εij - ∆εi - ∆εj
where εi is the correlation energy of the subsystem i and εij the correlation energy of the subsystem i and j together. For the convenient treatment of higher order terms we use a notation based on simple set theory.36 Now eq 1 reads
with C′L ) S1/2CL. In the canonical basis, E is diagonal and it contains the orbital energies, but in the local basis it is not diagonal. Now we classify the occupied orbitals into four classes, the frozen core orbitals, which are not correlated in all calculations, the environment orbitals, which are frozen in a specific calculation, the domain orbitals, which have to be correlated in this specific domain and the virtual orbitals, which are unoccupied in the HF reference. Furthermore we do not distinguish between core and environment orbitals, since they are treated equally in the subsequent steps. According to these criteria, we classify the blocks of the E matrix into nine blocks:
(
)
(2)
core domain virtual cd cv core cc ) E ) C′†F˜(CL)C′ (6) dd dv domain dc virtual vc vd vv
D is the set of domains, P(D) stands for the power set of the set of the domains and O denotes the order of the expansion. The summation index in eq 2 runs over all increments up to the order O (see refs 36 and 37 for details). When εX represents the correlation energy of the unified subsystems of the general correlation energy increment, ∆εX is given as
with (frozen core + environment) ) core. In our approach we diagonalize E in the subspace of the domain. In other words, we diagonalize the dd block of E. The unitary matrix U which is necessary to transform the MOs has the form
∑
Ecorr )
∆εX
X
X∈P(X)∧|X|eO
∆εX ) εX -
∑
∆εY
(3)
Y
Y∈P(X)∧|Y|
