RESEARCH ARTICLE
Advances in molecular quantum chemistry contained in the Q-Chem 4 program package
Yihan Shaoa, Zhengting Gana, Evgeny Epifanovskya,b,c, Andrew T.B. Gilbertd, Michael Wormite, Joerg Kussmannf, Adrian W. Langeg, Andrew Behnc, Jia Dengd, Xintian Fengb, Debashree Ghoshb,1, Matthew Goldeyc, Paul R. Hornc, Leif
D. Jacobsong, Ilya Kalimanh, Rustam Z. Khaliullinc, Tomasz Ku´sb, Arie Landaub,2, Jie Liui,g, Emil I. Proynova,3, Young Min Rheec,4, Ryan M. Richardg, Mary A. Rohrdanzg,5, Ryan P. Steelej, Eric J. Sundstromc, H. Lee Woodcock IIIad, Paul M. Zimmermanc,k, Dmitry Zuevb, Ben Albrechtl, Ethan Alguires, Brian Austinc, Gregory J. O. Beranm, Yves A. Bernardb,
Eric Berquistl, Kai Brandhorstc,6, Ksenia B. Bravayab,7, Shawn T. Browna,8, David Casanovac,n, Chun-Min Changa, Yunqing Chenk, Siu Hung Chiena, Kristina D. Clossera, Deborah L. Crittendend,9, Michael Diedenhofeno, Robert
A. DiStasio Jr.c, Hainam Dop, Anthony D. Dutoiq, Richard G. Edgarr, Shervin Fatehis,j, Laszlo Fusti-Molnara,10, An Ghyselst,11, Anna Golubeva-Zadorozhnayab, Joseph Gomesc, Magnus W.D. Hanson-Heinep, Philipp H.P. Harbache, Andreas W. Hauserc, Edward G. Hohensteinu, Zachary C. Holdeng, Thomas-C. Jagaub, Hyunjun Jiv, Benjamin Kadukw,
Kirill Khistyaevb, Jaehoon Kimv, Jihan Kimc,12, Rollin A. Kingx, Phil Klunzingery, Dmytro Kosenkovh,13, Tim Kowalczykw,14, Caroline M. Krautere, Ka Un Laog, Ad`ele D. Laurentb,15, Keith V. Lawlerc,16, Sergey
V. Levchenkob,17, Ching Yeh Lind, Fenglai Liual, Ester Livshitsz, Rohini C. Lochanc, Arne Luenserf, Prashant Manoharb,18, Samuel F. Manzerc, Shan-Ping Maoaa, Narbe Mardirossianc, Aleksandr V. Marenichab, Simon A. Maurerf, Nicholas J. Mayhallc, Eric Neuscammanc, C. Melania Oanab, Roberto Olivares-Amayar,ac, Darragh P. O’Neilld, John A. Parkhillc,19, Trilisa M. Perrinek,20, Roberto Peveratiab,c, Alexander Prociukk, Dirk R. Rehne, Edina Rostab,21, Nicholas J. Russa, Shaama M. Sharadac, Sandeep Sharmaac, David W. Smallc, Alexander Sodtc,t, Tamar Steinz,5, David St¨uckc, Yu-Chuan Suaa, Alex
J.W. Thomc,22, Takashi Tsuchimochiw, Vitalii Vanovschib, Leslie Vogtr, Oleg Vydrovw, Tao Wangb, Mark A. Watsonac,r, Jan Wenzele, Alec Whitec, Christopher F. Williamsg, Jun Yangac, Sina Yeganehw, Shane R. Yostc,w, Zhi-Qiang Youae,g, Igor
Ying Zhangaf, Xing Zhangg, Yan Zhaoab, Bernard R. Brookst, Garnet K.L. Chanac, Daniel M. Chipmanag, Christopher J. Cramerab, William A. Goddard IIIah, Mark S. Gordonai, Warren J. Hehrey, Andreas Klamto, Henry F. Schaefer IIIaj, Michael W. Schmidtai, C. David Sherrillu, Donald G. Truhlarab, Arieh Warshelb, Xin Xuaf, Al´an Aspuru-Guzikr, Roi Baerz,
Alexis T. Bellc, Nicholas A. Besleyp, Jeng-Da Chaiaa, Andreas Dreuwe, Barry D. Dunietzak, Thomas R. Furlanial, Steven R. Gwaltneyam, Chao-Ping Hsuae, Yousung Jungv, Jing Konga,3, Daniel S. Lambrechtl, WanZhen Liangi, Christian Ochsenfeldf, Vitaly A. Rassolovan, Lyudmila V. Slipchenkoh, Joseph E. Subotniks, Troy Van Voorhisw, John
M. Herbertg, Anna I. Krylovb, Peter M.W. Gilldand Martin Head-Gordonc,∗
aQ-Chem Inc, Pleasanton, CA, USA;bDepartment of Chemistry, University of Southern California, Los Angeles, CA, USA;cDepartment of Chemistry, University of California, Berkeley, CA, USA;dResearch School of Chemistry, Australian National University, Canberra,
Australia;eInterdisciplinary Center for Scientific Computing, Ruprecht-Karls University, Heidelberg, Germany;fDepartment of Chemistry, University of Munich (LMU), M¨unchen, Germany;gDepartment of Chemistry and Biochemistry, The Ohio State University,
Columbus, OH, USA;hDepartment of Chemistry, Purdue University, West Lafayette, IN, USA;iDepartment of Chemistry, Xiamen University, Xiamen, China;jDepartment of Chemistry, The University of Utah, Salt Lake City, UT, USA;kDepartment of Chemistry, University of Michigan, Ann Arbor, MI, USA;lDepartment of Chemistry, University of Pittsburgh, Pittsburgh, PA, USA;mDepartment of
Chemistry, University of California, Riverside, CA, USA;nIKERBASQUE - Basque Foundation for Science & Donostia International Physics Center (DIPC) & Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), Donostia, Spain;oCOSMOlogic GmbH & Co.
KG, Leverkusen, Germany;pSchool of Chemistry, The University of Nottingham, Nottingham, United Kingdom;qDepartment of Chemistry, University of the Pacific, Stockton, CA, USA;rDepartment of Chemistry and Chemical Biology, Harvard University, Cambridge, MA, USA;sDepartment of Chemistry, University of Pennsylvania, Philadelphia, PA, USA;tComputational Biophysics Section, Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD, USA;uSchool of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, GA, USA;vGraduate School of EEWS, KAIST, Daejeon, Republic of Korea;wDepartment of Chemistry, Massachusetts Institute of Technology, Cambridge, MA, USA;xDepartment of Chemistry, Bethel University, St. Paul, MN, USA;yWavefunction Inc., Irvine, CA, USA;zThe Chaim Weizmann Institute of Chemistry, and the Fritz Haber Research Center for Molecular Dynamics, The Hebrew University of Jerusalem, Jerusalem, Israel;aaDepartment of Physics, National Taiwan University, Taipei, Taiwan;abDepartment of Chemistry, University of Minnesota, Minneapolis, MN, USA;acDepartment of Chemistry, Princeton University, Princeton, NJ, USA;adDepartment of Chemistry, University of South Florida, Tampa, FL, USA;aeInstitute of Chemistry, Academia Sinica, Taipei, Taiwan;afDepartment of Chemistry, Fudan University, Shanghai, China;agRadiation Laboratory, University of Notre Dame, Notre Dame, IN, USA;ahMaterials and Process Simulation Center, California Institute of Technology, Pasadena, CA, USA;
aiDepartment of Chemistry, Iowa State University, Ames, IA, USA;ajDepartment of Chemistry, University of Georgia, Athens, GA, USA;
akDepartment of Chemistry and Biochemistry, Kent State University, Kent, OH, USA;alCenter for Computational Research, SUNY at Buffalo, Buffalo, NY, USA;amDepartment of Chemistry, Mississippi State University, Mississippi State, MS, USA;anDepartment of
Chemistry and Biochemistry, University of South Carolina, Columbia, SC, USA (Received 29 May 2014; accepted 1 August 2014)
⃝C 2014 Taylor & Francis
A summary of the technical advances that are incorporated in the fourth major release of the Q-CHEMquantum chemistry program is provided, covering approximately the last seven years. These include developments in density functional theory methods and algorithms, nuclear magnetic resonance (NMR) property evaluation, coupled cluster and perturbation theories, methods for electronically excited and open-shell species, tools for treating extended environments, algorithms for walking on potential surfaces, analysis tools, energy and electron transfer modelling, parallel computing capabilities, and graphical user interfaces. In addition, a selection of example case studies that illustrate these capabilities is given. These include extensive benchmarks of the comparative accuracy of modern density functionals for bonded and non-bonded interactions, tests of attenuated second order Møller–Plesset (MP2) methods for intermolecular interactions, a variety of parallel performance benchmarks, and tests of the accuracy of implicit solvation models. Some specific chemical examples include calculations on the strongly correlated Cr2dimer, exploring zeolite-catalysed ethane dehydrogenation, energy decomposition analysis of a charged ter-molecular complex arising from glycerol photoionisation, and natural transition orbitals for a Frenkel exciton state in a nine-unit model of a self-assembling nanotube.
Keywords: quantum chemistry; software; electronic structure theory; density functional theory; electron correlation;
computational modelling; Q-CHEM
Introduction
Quantum chemistry is a vigorous branch of theoreti- cal chemistry, which is concerned with the development of practical theory, algorithms, and software, based on approximations to the fundamental principles of quantum mechanics (QM). While the electronic Schr¨odinger equa- tion offers an in-principle exact description of the behaviour of electrons in molecules, subject to neglect of relativistic effects and nuclear motion, it is intractable to solve for real-
∗Corresponding author. Email: mhg@cchem.berkeley.edu
1Current address: Physical and Materials Chemistry Division, CSIR-National Chemical Laboratory, Pune 411008, India
2Current address: Faculty of Chemistry, Technion - Israel Institute of Technology, Technion City, Haifa 3200008, Israel
3Current address: Department of Chemistry, Middle Tennessee State University, 1301 East Main Street, Murfreesboro, TN 37132, USA
4Current address: Department of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang, 790-784, Republic of Korea
5Current address: Department of Chemistry, Rice University, Houston, TX 77251, USA
6Current address: Institut f¨ur Anorganische und Analytische Chemie, Technische Universit¨at Braunschweig, 38106 Braunschweig, Germany
7Current address: Department of Chemistry, Boston University, Boston, MA 02215, USA
8Current address: Pittsburgh Supercomputing Center, 300 S. Craig Street, Pittsburgh, PA 15213, USA
9Current address: Department of Chemistry, University of Canterbury, Christchurch, New Zealand
10Current address: OpenEye Scientific Software, 9 Bisbee Court, Suite D, Santa Fe, NM 87508, USA
11Current address: Center for Molecular Modeling, QCMM Alliance Ghent-Brussels, Ghent University, Technologiepark 903, B-9052 Zwijnaarde, Belgium
12Current address: Department of Chemical and Biomolecular Engineering, KAIST, Daejeon 305-701, Republic of Korea
13Current address: Department of Chemistry and Physics, Monmouth University, West Long Branch, New Jersey, 07764, USA
14Current address: Department of Chemistry, Western Washington University, 516 High St, Bellingham, WA 98225., USA
15Current address: University of Nantes, CEISAM UMR 6230 2 rue de la Houssini`ere, 44322 Nantes Cedex 1, France
16Current address: Department of Chemistry, University of Nevada, Las Vegas, NV 89154, USA
17Current address: Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin, Germany
18Current address: Birla Institute of Technology and Science, Pilani 333031, Rajasthan, India
19Current address: Department of Chemistry, and Biochemistry, University of Notre Dame, Notre Dame, IN 46556, USA
20Current address: Department of Chemistry and Biochemistry, Ohio Northern University, Ada, OH 45810, USA
21Current address: Department of Chemistry, King’s College London, London SE1 1DB, United Kingdom
22Current address: Department of Chemistry, Cambridge University, Lensfield Road, Cambridge CB2 1EW, United Kingdom
istic systems without approximations. Several fundamental approaches to developing such approximations have been followed. The predominant methods for present-day appli- cations to larger molecules are based on the framework of density functional theory (DFT). For smaller molecules, accuracy that is higher can be achieved by the use of wave function theory approaches such as perturbation theory, and coupled cluster (CC) theories. The optimal model for a given problem depends on the accuracy required, the
computational resources available, and the size of the sys- tem under consideration. In general, useful electronic struc- ture methods trade off accuracy against computational fea- sibility over a very wide range.
All of the approximate methods of quantum chemistry provide models by which the electronic potential energy of a molecule, E(R), can be evaluated as a function of the clamped nuclear positions,R. Walking on the potential en- ergy surface downwards to local minima leads to stable molecular structures, whose relative energies may be eval- uated to predict reaction energies, and thus the thermody- namics of chemical transformations. Walking downhill in all directions but one (the reaction coordinate), and walk- ing uphill in that direction leads to the first-order saddle points that separate reactants from products, and often play a major role in determining the kinetics of chemical re- actions. Multi-step reaction mechanisms can in principle be identified this way, with the aid of appropriate surface- walking algorithms. Molecular properties, many of which can be used for spectroscopic characterisation, may also be evaluated from quantum chemical models as derivatives of the energy with respect to applied perturbations, such as electric fields or magnetic fields.
Putting together a useful range of quantum chemical models that offer different trade-offs between achievable accuracy and computational effort for a range of molec- ular sizes is a non-trivial matter. Things are further com- plicated by the need to evaluate a range of responses of the energy to such key perturbations as moving the atoms, and applying fields. Therefore, the realisation of electronic structure simulations through useful software has evolved over the past five decades into team science of increasingly large scale. Early efforts such as Gaussian 70 represented essentially the work of a single group (Sir John Pople’s group). Today, there are roughly a dozen or so leading elec- tronic structure codes in chemistry, all of which represent the end result of delocalised collaborations amongst many groups. In addition to Q-CHEM, and its collaborator, Spar- tan (www.wavefun.com), leading commercial programs are represented by the ADF program [1], the Gaussian pro- gram [2], Jaguar [3], MolCAS [4], the Molpro package [5], and the TURBOMOLE program [6,7]. In addition there is a range of non-commercial programs which also repre- sent the result of substantial collaborations. These include ACES III [8], CFOUR [9], Dalton [10], GAMESS US [11]
and UK [12], NWChem [13], and Psi [14]. Many other related codes exist in the condensed matter physics com- munity, where periodic rather than molecular systems are typically the primary focus.
Some 21 years ago, in late 1992, Peter Gill, then a postdoctoral researcher with John Pople, began writing the first lines of a then-new quantum chemistry program, called Q-CHEM, over his Christmas vacation. This paper marks the fourth major release of the resulting software, which now is over 3 million lines of code, and contains a very
wide range of functionality for calculating the structure and properties of molecules using methods based on the principles of QM. The technical developments prior to 2000 were summarised in a first major review on Q-CHEMversion 2 [15], whose author list also illustrates the rapid growth in the number of contributors, which included not only members of the early founders’ groups, but also many new groups including most famously the 1998 Nobel Laureate, Sir John Pople [16]. Subsequent advances between 2000 and 2006 were contained in Q-CHEMversion 3.0, and were also documented in a review [17].
A very recent overview of Q-CHEM[18] provides some further details of the historical development and evolu- tion of the package, as well as a high-level summary of its capabilities. Today Q-CHEMserves the needs of a very large number of users (over 50,000 including both direct users, and the very large number of users who access its capabilities as the back-end of the widely used Spar- tan modelling package). Q-CHEM also serves the needs of one of the larger development communities in quan- tum chemistry, currently consisting of over 200 devel- opers spread across a large number of research groups, primarily in academia. For the developers of the code, Q-CHEM is an open team-ware project, where the source code is provided freely, and distributed and updated through a central code repository. The rights of other developers and the company itself are protected through a straightfor- ward non-disclosure agreement that places no restrictions on a developer’s ability to publish research describing new theory or algorithms. The activity of the developer com- munity is the key driver behind technical advances in the Q-CHEM software, so that this is very much a symbiotic relationship.
This paper summarises the fourth major release of Q-CHEM, and seeks to accomplish three principal pur- poses. The first purpose is to review a selection of the technical advances that have occurred in quantum chemistry over the past seven years or so which are incorporated into Q-CHEM 4. The review is, by neces- sity, relatively non-technical, with a focus on the physi- cal content of the methods and algorithms. We provide brief overviews of the strengths and weaknesses of a large and diverse selection of new methods from the perspec- tive of utility in chemical applications. Complete citations are given to the original literature for readers who are still hungry for further detail. The second purpose is to provide some example case studies of the new methods, partic- ularly those that are not widely used as yet. Such stud- ies provide some specific illustrations of the utility of the methods described in this work for particular chemical applications.
The third purpose of the paper is to serve as the litera- ture citation for release 4 of the Q-CHEMquantum chemistry software package. This purpose is useful because it leads via the technical review to full literature citations for the
key algorithms contained in the program, at the moment in time when this version is current. By contrast, websites are continually updated (and archival material is contin- ually removed), and an author list is an ‘empty citation’
that does not give the researcher any direct path to further information. The author list of this paper comprises the sci- entists who have contributed to Q-CHEMin either release 3 or release 4. Authors of earlier versions who have not sub- sequently contributed may be seen in the reviews describing release 2 [15] and release 3 [17].
The remainder of this paper addresses the challenge of reviewing the science and presenting selected examples under the following organisation. The main areas where methodological advances have been made within our code are reviewed together with a variety of example calcula- tions that illustrate accuracy and/or computational perfor- mance, and a selection of chemical case studies. The se- quence of topics begins with general purpose electronic structure methods. We treat DFT, which is the most widely used family of electronic structure methods, in Section2, and discuss recently added functionals for improved ac- curacy, and algorithmic improvements. Developments in wave function-based methods are reviewed in Section 3.
They have the great strength of systematic improvability, particularly at the level of CC theory, where object-oriented design is vital to facilitate development and implementation of new methods. Support for parallel and graphical process- ing unit (GPU) computing environments is summarised in Section4, with some example timings.
The standard (and many non-standard) electronic struc- ture methods can be used in a great many ways, start- ing with recent developments in moving around on the resulting potential energy surfaces, which are discussed in Section5. We turn next to the problem of treating ex- tended environments in Section6, which is important for modelling molecules in solution, large clusters, or active sites abstracted from complex systems such as proteins or heterogeneous solids. Energy and electron transfer capa- bilities are discussed in Section 7, followed by methods for chemical analysis of (at least some classes of) cal- culations in Section 8. As appropriate for a review of a
‘back-end’ code, as Q-CHEMfundamentally is, we then fin- ish with a short discussion of the available ‘front-ends’ that provide input to the back-end, and visualise the resulting output.
Density functional theory Functionals
Kohn–Sham DFT (KS-DFT)[19] and its extensions provide a foundation for the development of model functionals, but no prescription for how such development should be accomplished. Accordingly, this is an area of great activity.
The full range of density functionals supported in Q-CHEM
is too large to comfortably list here, and includes functionals ranging from vintage to brand new. Furthermore, assessing the strengths and weaknesses of different functionals is a major ongoing effort that involves the entire community of both developers and applications specialists. An overview paper cannot summarise this effort, although we can provide a few leading references to comparative studies [20–23] and the main issues [24,25].
We shall discuss some of the new functionals added to Q-CHEMover recent years by considering the current main directions for improved physical content over the predomi- nant density functional for chemistry, which during the last decade was certainly the global hybrid, B3LYP [26].
(1) Meta-generalised gradient approximation (GGA) and hybrid meta-GGA functionals: including the kinetic energy density, τ , gives flexibility beyond global hybrids, and therefore is used in many modern density functionals, including M06-L, M06, M06-2X, and M06-HF [27], as well as the recently introduced M11 and M11-L functionals [28,29]. These functionals often yield improved ac- curacy for thermochemistry (TC) and non-covalent (NC) interactions relative to functionals that do not depend on τ . M06-L and M11-L have the con- siderable computational advantage of not requiring exact exchange, although there is reduced accu- racy, particularly for reaction barrier heights. Other meta-GGA functionals include the constraint-based TPSS [30], as well as its one-parameter hybrid cousin, TPSSh [30].
(2) Range-separated hybrid functionals: self- interaction error (SIE), where an electron artificially sees a fraction of itself, is a well-known defect of standard density functionals, and causes artefacts that include spurious delocalisation of unpaired electrons [31,32], and charge-transfer excited states that can be drastically too low [33,34]. While very difficult to remove fully, SIE can be significantly reduced by including 100%
exact (wave function) exchange at large electron–
electron distances, and a much smaller fraction at short distances, where DFT exchange functionals are effective. Examples of functionals of this type include the LC-ωPBE family of methods [35–38], the ωB97 functionals [39–42], M11 [28], as well as tuned functionals of the BNL type [43,44], where the range-separation parameter can be chosen for the problem at hand based on physical criteria [45,46]. Another option is to include 100%
Hartree–Fock (HF) exchange at all inter-electronic distances, as in the M06-HF functional [47]. SIE at short inter-electronic distances also affects
TDDFT predictions of core excitation energies and near-edge X-ray absorption fine structure spectra.
These short-range SIE errors can be substantially reduced by short-range corrected functionals [48], which are available in Q-CHEM.
(3) Non-covalent interactions: NC interactions, partic- ularly van der Waals forces, involve non-local cor- relation effects that are very difficult to treat within a standard correlation functional. Thus, the char- acteristic R−6 long-range interaction potential is absent in traditional local and semi-local density functionals [49]. While it is quite often possible to still obtain an accurate result at the van der Waals minimum without including the long-range R−6be- haviour (e.g., with M06-2X or M06-L), several vi- able methods have emerged that recover the correct long-range behaviour [50]:
• A vast range of dispersion-corrected functionals that include damped C6/R6 atom–atom poten- tials, based on either the Grimme-D2 [51] or -D3 [52] parameterisations are available. Computa- tionally virtually free, but not actually density functionals at all, these methods represent the simplest possible treatment of dispersion.
• Becke’s exchange dipole model [53,54] is a novel and accurate method for treating disper- sion which has been implemented in Q-CHEMin an efficient and numerically stable form [55].
• van der Waals density functionals, which numer- ically integrate a nonlocal correlation functional that depends simultaneously on ρ(r) and ρ(r′), are a soundly based approach. Examples include vdW-DF-04 [56], vdW-DF-10 [57], VV09 [58], and VV10 [59]. The VV10 form is also used in the very recently developed ωB97X-V functional [42], a 10-parameter semi-empirical functional that is a further evolution of the ωB97 family that reduces the number of empirical parameters, while improving physical content.
(4) Double hybrid functionals: based on G¨orling–Levy perturbation theory as well as semi-empirical con- siderations, double hybrid functionals (also some- times called ‘doubly hybrid’ functionals) include second-order perturbation theory (PT2) corrections to a KS reference. They can yield improved ac- curacy for both bonded and non-bonded interac- tions, albeit with increased computational cost and a need for larger basis sets. Q-CHEMcontains nu- merous double hybrids, including B2PLYP [60]
and B2PLYP-D [61], XYG3 [62], ωB97X-2 [41], and PBE0-2 [63]. XYGJ-OS is an opposite spin double hybrid [64], which scales as only O(M4),
for which the analytical gradient is also available [65].
(5) Becke post-self-consistent field (SCF) function- als: the semi-empirical B05 post-SCF functional [66,67] uses the Becke–Roussel exchange model [68] to compute a local analogue to the exact exchange hole. The extent of delocalisation of the exact exchange hole is used as a parame- ter for capturing both same-spin and opposite- spin non-dynamical correlation within a single de- terminant framework. Coupled with the modified Bc88 correlation functional (BR94) [69,70] to cap- ture dynamical correlation, the performance of the six-parameter B05 functional parallels the perfor- mance of existing hybrid meta-GGA functionals for atomisation energies and barrier heights. A self-consistent version of the B05 functional has been efficiently implemented [71] into Q-CHEM4, including a resolution-of-the-identity (RI) version that greatly reduces the cost of computing the exact exchange energy density [72].
To provide just a glimpse of the comparative perfor- mance of some of the standard functionals available in Q-CHEM,Table 1shows the root mean square (RMS) errors associated with a variety of density functionals on some es- tablished test sets for bonded and non-bonded interactions.
Of the 10 data-sets, the first 4 (TAE, Alk19, DBH24, and G21IP) correspond to TC datapoints (203 total), while the latter 6 (HW30, S22, S66, A24, X40, and DS14) correspond to NC interaction datapoints (196 total). Computational de- tails and specific data-set information can be found in [42];
comparisons of GGA functionals trained on the same data with different choices of non-local exchange and correla- tion have also been presented recently [73].
For the bonded interactions, it is clear that exact ex- change is very useful, as the two best functionals are the range-separated hybrid GGA ωB97X-V functional and the global hybrid meta-GGA M06-2X functional, with RMS errors around 3 kcal/mol. In comparison, the best local functional for TC (M06-L) has an RMS error that is nearly double that of the best hybrid functional. For non-bonded interactions, the range-separated hybrid GGA ωB97X-V functional outperforms the next best density functional by almost a factor of 2. However, local meta-GGA functionals like M06-L can compete with well-established dispersion- corrected hybrid functionals such as ωB97X-D. On the popular S22 data-set, the three best density functionals are ωB97X-V, ωB97X-D, and M06-L, with RMS errors of 0.23, 0.41, and 0.43 kcal/mol, respectively. For a more comprehensive assessment of these density functionals on a data-set of over 2400 datapoints, the reader is referred to Table 6in [42].
Table 1. RMS errors in kcal/mol for 17 density functionals available in Q-CHEMon 10 data-sets comprising 399 datapoints. The first four data-sets contain bonded interactions, representative of thermochemistry (TC), evaluated using the aug-cc-pVQZ basis. The latter six data-sets contain non-covalent (NC) interactions, which are evaluated using the aug-cc-pVTZ basis, except for the X40 data, which used the def2-TZVPPD basis set. The TC data-sets are (1) non-multi-reference total atomisation energies (TAE) [74], (2) atomisation energies of n = 1–8 alkanes (Alk19) [75], (3) diverse barrier heights (DBH24) [76–78], (4) adiabatic ionisation potentials (G21IP) [79,80], (5) weak hydrocarbon–water interactions (HW30) [81], (6) hydrogen-bonded and dispersion-bonded complexes (S22) [82,83], (7) interaction energies of relevant biomolecular structures (S66) [84,85], (8) small non-covalent complexes (A24) [86], (9) non-covalent interactions of halogenated molecules (X40) [87], and (10) interactions of complexes containing divalent sulphur (DS14) [88]. The final two columns give the overall RMS errors for the four TC data-sets and six NC data-sets.
TAE Alk19 DBH24 G21IP HW30 S22 S66 A24 X40 DS14 TC NC
# 124 19 24 36 30 22 66 24 40 14 203 196
PBE-D2 1 16.94 26.21 10.37 4.81 0.71 0.70 0.63 0.59 0.74 0.57 16.01 0.67
PBE-D3 2 16.85 20.93 10.27 4.81 0.48 0.60 0.46 0.41 0.59 0.47 15.20 0.50
B3LYP-D2 4 5.28 0.64 5.28 4.86 0.53 0.74 0.62 0.39 0.47 0.28 4.95 0.55
B3LYP-D3 5 5.23 5.50 5.23 4.86 0.35 0.50 0.43 0.23 0.34 0.23 5.19 0.38
B3LYP-NL 4 5.92 14.74 6.01 4.35 0.21 0.68 0.46 0.27 0.43 0.24 7.03 0.43
B97-D2 11 4.06 9.28 4.36 3.48 0.35 0.60 0.36 0.26 0.43 0.25 4.74 0.39
B97-D 9 5.18 10.48 7.18 4.47 0.40 0.54 0.52 0.32 0.59 0.37 6.03 0.49
VV10 2 12.46 5.85 9.86 5.43 0.43 0.63 0.52 0.41 0.63 0.52 10.71 0.53
LC-VV10 3 5.30 19.04 3.02 5.23 0.30 0.51 0.31 0.15 0.41 0.12 7.55 0.34
ωB97X 14 3.50 2.84 2.33 3.79 0.45 0.95 0.50 0.40 0.54 0.36 3.38 0.55
ωB97X-D 15 3.65 2.90 2.07 3.82 0.35 0.41 0.52 0.15 0.49 0.18 3.47 0.43
ωB97X-V 10 3.34 0.71 1.81 3.57 0.20 0.23 0.18 0.09 0.21 0.05 3.09 0.18
M06-L 34 5.54 8.11 5.38 5.60 0.35 0.43 0.36 0.23 0.48 0.25 5.82 0.38
M06 33 3.94 4.63 2.97 3.78 0.33 0.77 0.53 0.25 0.57 0.34 3.88 0.51
M06-2X 29 3.24 5.27 1.12 3.49 0.46 0.47 0.29 0.28 0.28 0.20 3.37 0.34
M11-L 44 6.62 29.35 3.54 4.54 0.48 0.91 0.81 0.46 1.23 0.59 10.61 0.84
M11 40 4.37 3.94 1.48 4.64 0.38 0.58 0.41 0.27 0.54 0.30 4.15 0.44
2.2. Electric and magnetic molecular properties The calculation of molecular properties provides an im- portant link to experiment and many linear-scaling meth- ods have been developed over recent years that allow to compute molecular systems with more than 1000 atoms (see, e.g., [89–92]). The Q-CHEMprogram package allows the computation of a wide range of properties at the HF and KS-DFT levels. Apart from determining geometries, vibra- tional spectra, and electronic excitations, the new version of Q-CHEMoffers several new and improved efficient linear- scaling methods to evaluate different electric and magnetic response properties for large systems. These range from the calculation of NMR chemical shieldings using density matrix-based coupled-perturbed SCF theory [93,94] (also for large basis sets with up to g functions in the new version) to electric response properties [95,96] by an implementa- tion of the density matrix-based time-dependent SCF al- gorithm [97] that allows for calculating static and dynamic polarisabilities and first hyperpolarisabilities. Here, an over- all asymptotic linear-scaling behaviour can be reached by employing O(N) integral evaluations based on CFMM [98]
and LinK [99,100] in combination with efficient sparse al- gebra routines [93].
Furthermore, the combination of linear-scaling QM methods for calculating molecular properties with simple
molecular mechanics (MM) schemes (QM/MM) has proven to be a very valuable tool for studying complex molecular systems. The linear-scaling methods allow to systemati- cally converge the property with the chosen QM sphere, and convergence for QM/MM schemes is typically clearly faster than in pure QM calculations, since in complex sys- tems long-range electrostatics are accounted for (see, e.g., [92,101,102]).
As another new feature, the calculation of indirect nu- clear spin–spin coupling constants (J-coupling) [103,104]
is introduced. The implementation uses the LinK scheme [99,100] for the construction of exchange-type matrices for non-metallic systems. A fully density matrix-based algo- rithm is currently in development [105]. Basis functions with angular momenta up to g are supported. Predictions of good accuracy for J-couplings can be obtained [106,107]
especially when using specialised basis sets [108–111], sev- eral of which have been added to the basis set library.
The J-based configurational analysis [112] is a robust technique for the structural elucidation of even large organic molecules [113]. In the past, analyses of J-couplings have mostly utilised the Karplus equations [114–116], which re- late3J-coupling constants to the dihedral angle between atoms. Moving from empirical equations to predictions from first principles, for example with DFT, is desirable
not only because the predictions are expected to be more reliable, but also because it expands their applicability.
Long-range J-couplings [117], or even couplings through hydrogen bonds [118], neither of which are tractable with current empirical methods, can naturally be studied by com- putational chemistry if an adequate level of electronic struc- ture theory is selected.
2.3. Algorithm developments
A great variety of algorithmic improvements for DFT calculations have been incorporated into Q-CHEM. Those to do with parallel computing are covered separately in Section4.
Resolution-of-the-identity methods: for large atomic or- bital (AO) basis sets in particular, RI methods offer substan- tially improved performance relative to exact evaluation of two-electron integrals with nearly negligible loss of ac- curacy, if appropriately optimised auxiliary basis sets are employed. The Karlsruhe group associated with TurboMole has been amongst the leading developers of such basis sets.
The standard RI method may be further enhanced by per- forming local fitting of a density or function pair element.
This is the basis of the atomic-RI method (ARI), which has been developed for both Coulomb (J) matrix [119] and exchange (K) matrix evaluation [120]. In ARI, only nearby auxiliary functions K(r) are employed to fit the target func- tion. This reduces the asymptotic scaling of the matrix- inversion step as well as that of many intermediate steps in the digestion of RI integrals. Briefly, atom-centred auxiliary functions on nearby atoms are only used if they are within the outer radius (R1) of the fitting region. Between R1and the inner radius (R0), the amplitude of interacting auxiliary functions is smoothed by a function that goes from zero to one and has continuous derivatives. To optimise efficiency, the van der Waals radius of the atom is included in the cut- off so that smaller atoms are dropped from the fitting radius sooner. Energies and gradients are available.
Multi-resolution exchange-correlation (MrXC) quadra- ture: MrXC [121–123] can accelerate the numerical quadra- ture associated with computation of the XC energy and the XC matrix needed in the SCF procedure. It is an algorithm for seamlessly combining the standard atom-centred grid of quantum chemistry, with a cubic grid of uniform spacing by placing the calculation of the smooth part of the density and XC matrix onto the uniform grid. The computation as- sociated with the smooth fraction of the electron density is the major bottleneck of the XC part of a DFT calculation and can be done at a much faster rate on the cubic grid due to its low resolution. Fast Fourier transform and B-spline interpolation are employed for the accurate transformation between the two types of grids such that the final results remain the same as they would be on the atom-centred grid alone. By this means, a speedup of several times for the cal- culations of the XC matrix is achieved. The smooth part of
the calculation with MrXC can also be combined with the Fourier transform Coulomb method [124] to achieve even higher efficiency, particularly for calculations using large basis sets and diffuse functions.
TDDFT gradients and Hessians: a recent implementa- tion of TDDFT analytical gradients (also in the Tamm–
Dancoff approximation [125]) is available [126], with parallel capabilities. A more distinctive capability is the availability of an implementation of TDDFT analytical frequencies (greatly extending an existing analytical con- figuration interaction with single substitutions [CIS]
frequency code [127]), both in the Tamm–Dancoff approxi- mation [128] and for full TDDFT [129]. In addition, analyt- ical TDDFT gradients and frequencies have been extended [130] to include the smooth polarisable continuum models for solvation that are discussed in Section 6.1. Compared to numerical differentiation, analytical second derivatives of the excitation energy yield higher precision and need much less computer time, but require much more memory.
The memory usage is mainly dominated by the geometric derivatives of the MO coefficients and transition ampli- tudes, which is dealt with by solving the coupled-perturbed equations in segments. To ensure high precision, a fine grid for numerical integration should be used, since up to fourth-order functional derivatives with respect to the den- sity variables as well as their derivatives with respect to the nuclear coordinates are needed.
Dual basis methods: an effective method for reducing the computational cost of SCF calculations is to perform the SCF calculation in a small (‘primary’) basis set, and subsequently correct that result in a larger (‘secondary’) basis set, using perturbation theory. Q-CHEMcontains two related approaches for performing dual basis calculations.
Both energies and analytical gradients are available for the dual basis approach of Head-Gordon and co-workers [131,132], who have also developed dual basis pairings for the Dunning cc-pVXZ [133] and aug-cc-pVXZ [134] basis sets. From a good reference, the HF perturbation theory ap- proach [135,136] provides more accurate corrections even at first order, and is also applicable to DFT calculations. Not only jumps in basis set, but also in the choice of quadrature grid, and density functional itself are possible via the ‘triple jumping’ approach [137].
Metadynamics and non-orthogonal configuration inter- action (NOCI): SCF metadynamics [138] allows one to find alternative minima within the SCF framework (either HF or KS-DFT). Alternative minima are obtained by applying a bias in density matrix space at the locations of previ- ously found minima and using standard convergence algo- rithms on this modified potential energy surface. It is then possible to perform NOCI [139,140] using the resulting non-orthogonal determinants as a basis. One then builds and diagonalises the Hamiltonian in this representation.
Q-CHEM supports the use of general and complex HF orbitals for this purpose. While calculation of the
1 2
[kcal/mol]
PBE PBE-D2
-1 0
e Energy Error [kca
PBE-D2 B3LYP B3LYP-D2 ωB97X-D ωB97X-V
-3 -2
Rela"ve En ωB97X-V
M06-L M06-2X M11 Isomer
15 20
mol] PBE
0 5 10
rgy Error [kcal/m PBE
PBE-D2 B3LYP B3LYP-D2 ωB97X-D
-15 -10 -5 0
Binding Energy
B3LYP D2 ωB97X-D ωB97X-V M06-L M06-2X -20
-15
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
Bin
Isomer
M06-2X M11
Figure 1. Relative and binding energy errors for 10 isomers of F−(H2O)10with respect to RI-CCSD(T)/CBS benchmark values.
Hamiltonian is more complicated than in orthogonal CI, it has been shown for some systems that the number of deter- minants to obtain qualitatively accurate results for ground and excited states of challenging systems such as polyenes is rather small (less than 100 or so) [140].
2.4. Case study: relative and binding energies of 10 F−(H2O)10isomers
In a recent paper [141], Herbert and co-workers discov- ered that halide–water clusters present a challenge for density functionals such as LC-VV10, ωB97X-D, and M06-2X. In particular, the binding energies of 10 isomers of F−(H2O)10 proved to be the most notorious case. Us- ing the same geometries and reference values, nine den- sity functionals available in Q-CHEMwere benchmarked on these 10 isomers, and the results, calculated in the aug-cc- pVTZ basis set with a (99,590) grid, are provided inTable 2 and Figure 1. Table 2 confirms that a majority of these
Figure 2. Timings and parallel speedup of an ADC(3) calcula- tion of the benzene molecule using aug-cc-pVTZ basis set.
density functionals are unable to accurately predict the binding energies of these isomers. However, the newly de- veloped ωB97X-V functional performs at least two times better than the next best functional, which is (surprisingly) PBE. In order to identify if the functionals in question are underbinding or overbinding the clusters, it is useful to considerFigure 1. Besides PBE and B3LYP, all of the density functionals overbind the isomers, with B3LYP-D2, M06-2X, and PBE-D2 overbinding more severely than the rest. ωB97X-D, M06-L, and M11 overbind considerably as well, though by approximately 5 kcal/mol instead of more than 10 kcal/mol. For the relative energies of the clusters, the Minnesota functionals perform poorly, with errors larger than 1 kcal/mol. The two best functionals are ωB97X-V and ωB97X-D, with RMS errors of 0.40 and 0.45 kcal/mol, respectively. The parallel performance of ADC(3) calculations on the benzene molecule is illustrated in Figure 2.
3. Wave function methods 3.1. Perturbative methods
Second-order Møller–Plesset perturbation theory [142,143]
is widely used as the simplest and most computation- ally inexpensive wave function treatment of dynamic cor- relation. Q-CHEM’s workhorse implementation based on RI algorithms for the energy and gradient [144,145] is
Table 2. Relative and binding energy RMS errors in kcal/mol for 10 isomers of F−(H2O)10with respect to RI-CCSD(T)/CBS benchmark values.
RMSD PBE PBE-D2 B3LYP B3LYP-D2 ωB97X-D ωB97X-V M06-L M06-2X M11
Relative energy 0.68 0.63 0.68 0.54 0.45 0.40 1.10 0.97 1.35
Binding energy 4.11 14.19 14.42 11.13 5.76 1.94 6.17 11.30 5.74
Table 3. CPU times on a single CPU core and scaling behaviour for conventional RI-MP2 as well as RI-CDD-MP2 calculations on model DNA systems in a def2-SVP basis. The index of the DNA systems denotes the number of A–T base pairs. For full details and additional performance data can be found in [148].
# Basis RI-MP2 RI-CDD-MP2
System Functions Time (h) Scaling Time (h) Scaling
DNA1 625 0.16 – 0.23 –
DNA2 1332 6.36 4.87 4.75 4.02
DNA4 2746 231.63 4.97 53.22 3.34
DNA8 5574 – – 449.53 3.01
highly efficient for small- and medium-size molecules.
For greater efficiency in larger basis sets, energies [133]
and gradients [146] are also available for the dual basis RI-MP2 method. For larger molecules, an efficient cubic- scaling MP2 method has been implemented in Q-CHEM. The method is grounded on the atomic orbital-based MP2 formulation and uses a Cholesky decomposition of pseudo- density matrices (CDD) [147,148] in combination with in- tegral screening procedures using QQR integral estimates [148–150]. Using the RI approach and efficient sparse ma- trix algebra, the RI-CDD-MP2 method shows a fairly small prefactor for a reduced-scaling method. Due to the asymp- totically cubic scaling of the computational cost of the RI-CDD-MP2 method with the size of the molecule, the approach is faster for larger systems than the conventional fifth-order scaling RI-MP2 method. The crossover between RI-CDD-MP2 and conventional RI-MP2 is found already for systems as small as, e.g., two DNA base pairs as shown by the timings inTable 3.
While MP2 greatly improves on the mean-field refer- ence in many cases, it also has some well-known weak- nesses. These include a need for large basis sets, overes- timation of intermolecular interactions, and susceptibility
to spin-contamination. Q-CHEM contains a variety of re- cently developed methods that partially lift some of these limitations. For ground state treatment of intermolecular interactions, Q-CHEMcontains newly developed attenuated MP2 methods, which offer remarkable improvements in accuracy for small- and medium-sized basis sets. Attenu- ated MP2 [151–153] is available with the Dunning aug-cc- pVDZ (small) and aug-cc-pVTZ (medium) basis sets. This approach works by cancelling the overestimation of inter- molecular interactions by attenuation of the long-range part of the correlation energy. A summary of the RMS errors (kcal/mol) obtained for a series of inter- and intramolec- ular non-bonded interactions is given inTable 4. Consis- tent improvement relative to unattenuated MP2 is found for databases of hydrogen-bonded, dispersion, and mixed interactions (divalent sulphur, A24, S22, S66, and L7).
Relative conformational energies for sulphate–water clus- ters, alkane conformers (ACONF), cysteine conformers (CYCONF), sugar conformers (SCONF), and dipeptide and tripeptide conformers (P76) are in good agreement with benchmarks.
As a further evolution of spin-component scaled MP2 methods for systems susceptible to spin-contamination,
Table 4. Root mean squared errors (in kcal/mol) for databases of non-bonded interactions, grouped by intermolecular or intramolecular interactions. Only equilibrium geometries were examined from the divalent sulphur database [154]. Complete basis set estimates (CBS) for MP2 were taken from references for the divalent sulphur, SW49 [155–157], ACONF [158], CYCONF [159], and SCONF [160] databases.
MP2/CBS results for the S22 [82,83], S66 [84,85], L7 [161], and P76 [162] databases were obtained from the Benchmark Energy and Geometry DataBase (BEGDB) [163]. MP2/CBS results for the A24 databases [164] were generated for this work.
MP2/aDZ MP2(terfc, aDZ) MP2/aTZ MP2(terfc, aTZ) MP2/CBS
Divalent sulphur 1.25 0.28 0.80 0.16 0.41
A24 0.52 0.26 0.31 0.18 0.21
S22 3.91 0.61 2.5 0.48 1.39
S66 2.66 0.43 1.53 0.25 0.73
L7 24.14 1.10 14.00 1.87 8.78
SW49(bind) 1.23 1.03 0.84 0.36 0.34
SW49(rel) 0.49 0.34 0.34 0.39 0.10
ACONF 0.31 0.29 0.24 0.08 0.11
CYCONF 0.20 0.28 0.30 0.21 0.25
SCONF 0.28 0.52 0.22 0.12 0.21
P76 1.06 0.33 0.59 0.31 0.42
the orbital optimised opposite spin (O2) method is avail- able (energies [165] and gradients [166]). Relative to full OO-MP2, which exhibits systematic overbinding, O2 yields higher accuracy by virtue of its single semi-empirical scal- ing parameter, and also lower computational cost (O(M4) vs. O(M5)) by virtue of containing no same-spin contribu- tion [167]. In addition to cleaning up spin-contamination problems [165], O2 also avoids the n-representability and force discontinuity issues of MP2 [168]. At higher compu- tational cost than O2 (same cost as OO-MP2), the recently introduced δ-OO-MP2 method [169] simultaneously solves the problems of overestimating correlation effects and di- vergences from vanishing denominators in small-gap sys- tems by using a regularisation, or level shift, parameter.
Beyond ground states, the corresponding second-order correction to single excitation CI for excited states [170], CIS(D), offers similar advantages to MP2 and suffers from very similar limitations. Within Q-CHEM, more computa- tionally efficient excited state scaled opposite-spin meth- ods are available, that, like SOS-MP2 and O2, scale as O(M4). SOS-CIS(D), a non-degenerate method [171], is available for excited state energies. SOS-CIS(D0), a quasi- degenerate approach, which is therefore more robust at the cost of some additional computation, has both energies [173] and analytic gradients [174,175] available.
3.2. Coupled cluster methods
Q-CHEM4 features a wide variety of computational meth- ods for the ground and excited states based on CC theory [176,177]. These methods are amongst the most versatile and accurate electronic structure approaches. The equation- of-motion (EOM) approach extends single-reference CC methods to various multi-configurational wave functions.
Q-CHEM 4 includes EOM-CC methods for electronically excited states (EOM-EE), ionised/electron-attached ones (EOM-IP/EA), as well as doubly ionised (EOM-DIP) [178]
and spin-flip (EOM-SF and EOM-2SF) [179,180] exten- sions that enable robust and reliable treatment of bond- breaking, diradicals/triradicals, and other selected multi- configurational wave functions. Gradient and properties calculations (including interstate properties) are available for most CC/EOM-CC methods.
Q-CHEM4 offers an efficient multi-core parallel imple- mentation of these methods based on a general purpose tensor library [181]. The library provides a convenient ten- sor expressions C++ interface that aids new developments.
In order to reduce computational requirement for the CC methods and improve parallel performance, we exploited two reduced-rank approaches based on RI and Cholesky de- composition (CD) of two-electron repulsion integrals [172].
The equations were rewritten to eliminate the storage of the largest four-dimensional intermediates leading to a signif- icant reduction in disk storage requirements, reduced I/O penalties, and, as a result, improved parallel performance.
Table 5. Timings of one CD-CCSD iteration (in hours) for (mU)2-H2O (test 4 in [172]) using 1E-2 threshold for Cholesky decomposition. This calculation takes 12 CCSD iterations to converge.
# Basis Memory Wall
Method Basis functions limit time
CD-CCSD 6-31+G(d,p) 489 100 GB 5.1
CD-CCSD/FNO 6-31+G(d,p) 489 100 GB 1.4
cc-pVTZ/FNO cc-pVTZ 882 300 GB 12.2
For medium-size examples, RI/CD calculations are ap- proximately 40%–50% faster compared with the canonical implementation. More significant speedups (two- to five- fold) are obtained in larger basis sets, e.g., cc-pVTZ.
Even more considerable speedups (six- to seven-fold) are achieved by combining RI/CD with the frozen natural orbitals approach [182]. Importantly, with Q-CHEM, one can perform CC/EOM-CC calculations for relatively large sys- tems (up to ∼1000 basis functions) on mainstream single- node servers. Detailed performance benchmarks are avail- able in [172,181].Table 5shows selected timings obtained on a single 16-core Xeon-Dell node for dimethyl-uracyl dimer solvated by one water molecule ((mU)2-H2O, C1 symmetry, 158 electrons). Frozen natural orbital (FNO) cal- culations inTable 5used an occupation threshold of 99.5%.
As an example, for the 6-31+G(d,p) basis, this corresponds to 292 active virtual orbitals and 118 frozen virtuals. Using FNO leads to errors in IEs that are less than 0.02 eV rela- tive to the full calculation, which is typical for this threshold [182].
While conventional CC and EOM methods allow one to tackle electronic structure ranging from well-behaved closed-shell molecules to various open-shell and electron- ically excited species [177], metastable electronic states, so-called resonances, present a difficult case for theory.
By using complex scaling and complex absorbing potential techniques, we extended these powerful methods to describe autoionising states, such as transient anions, highly excited electronic states, and core-ionised species [183–185]. In addition, users can employ stabilisation techniques using charged sphere and scaled atomic charges options [186].
Various improvements of iterative diagonalisation algo- rithms enable access to high-lying interior eigenvalues.
3.3. Algebraic diagrammatic construction (ADC) methods
Algebraic diagrammatic construction (ADC) methods con- stitute a series of methods for the calculation of excited states which derive from the perturbation expansion of the polarisation propagator [187,188]. Each method dif- fers in the approximation to the Hamiltonian matrix for which the eigenvalue problem has to be solved, as well
as in the way state properties and transition properties are computed from the eigenvectors. In first order, the ADC(1) eigenvalue problem is identical to CIS, but additional terms enter in the computation of transition moments. The second- order approximation, ADC(2) provides excitation energies (and transition properties) comparable to those obtained with CIS(D) [170] and CC2 [189,190]. An extension to the second-order scheme ADC(2)-x adds additional terms to the Hamiltonian matrix which put more weight on doubly excited configurations. As a result, the excitation energies are shifted to lower energies, in particular if the respective states possess strong double excitation character. Accord- ingly, the comparison of ADC(2)-x and ADC(2) results can yield useful insights about the importance of double exci- tations in the spectrum [191]. With the third-order method, ADC(3) the accuracy of the excitation energies improves further, getting close to the results obtained by CC3 calcu- lations, though at computational costs which are an order of magnitude smaller.
In Q-CHEM, ADC methods up to third order are available for the computation of excited states [192]. They have been implemented based on the same general purpose tensor library [181] as the CC methods, offering shared-memory parallelisation and a LaTeX style programming interface for new equations. The implementation allows for the calcula- tion of excitation energies and transition properties from the ground state, as usual. In addition, excited state properties and transition properties between excited states can be com- puted on request. From those, two-photon absorption cross- sections can be deduced via sum-over-states expressions.
Alternatively, the two-photon absorption cross-section can be obtained by inversion of the ADC matrix [193]. For visualisation of the excited states, transition densities or at- tachment and detachment densities may be exported as grid data for later display by standard visualisation tools.
Furthermore, Q-CHEM features spin-opposite scaled (SOS) ADC variants for both second-order schemes ADC(2) and ADC(2)-x [194]. They follow the idea of SOS-MP2 to reduce computational costs and improve the resulting energies. Therefore, same-spin contributions in the ADC matrix are neglected, while opposite-spin con- tributions are scaled using appropriate semi-empirical pa- rameters. SOS-ADC(2) requires two scaling parameters cos and cc, while for SOS-ADC(2)-x another parameter cx is needed. The parameter cos = 1.3 is inherited from SOS-MP2 for the scaling of the T2 amplitudes, while the parameters ccand cxare used to scale the ph/2p2h block and the off-diagonal part of the 2p2h/2p2h block of the ADC matrix, respectively. For SOS-ADC(2) the optimal value of cc= 1.17 was determined by fitting against the Thiel bench- mark set [195]. A similar fit for SOS-ADC(2)-x yielded cc= 1.0 and cx= 0.9 as optimal values [194]. With these parameters, a mean absolute error of 0.14 eV in the excita- tion energies is achieved by SOS-ADC(2) for the Thiel benchmark set. For SOS-ADC(2)-x, the mean absolute
error for predominantly single excitations becomes 0.17 eV, while for states with large double excitation character, it is 0.21 eV.
Another set of ADC variants in Q-CHEMuses the core–
valence separation (CVS) approximation [196] to calculate core excitations. In general, the calculation of core-excited states is quite difficult, since with standard implementations the valence excited states need to be calculated before any core-excited state can be obtained. The CVS approxima- tion solves this problem by decoupling core and valence excitations in the ADC Hamiltonian. Thereby, it makes use of the fact that the interactions between core and valence excitations are negligible due to the strong localisation of the core orbitals and the large energy separation between core and valence orbitals. As result, core-excited states can be computed independently from the core-excitation part of the ADC Hamiltonian, which significantly reduces the computational costs compared to the calculation of valence- excited states. CVS variants for ADC(1), ADC(2), and ADC(2)-x are available with CVS-ADC(2)-x showing ex- cellent agreement with experimental data [197].
3.4. Density matrix renormalisation group
Q-CHEM now includes an interface to the density matrix renormalisation group (DMRG) code of Sharma and Chan (‘Block’) [198–201]. The DMRG [198–211] is a varia- tional wave function method based on a class of wave func- tions known as matrix product states (MPS) [208,212]. The DMRG allows for unique kinds of quantum chemical cal- culations to be performed. The accuracy of the DMRG can be continuously tuned based on a single parameter, the number of renormalised states, denoted M. In typical cal- culations, M ranges from 1000 to 10,000. By increasing M, it is possible to push DMRG calculations to yield highly accurate energies (e.g., to within 10–100 microHartrees of the exact result) for systems much larger than can be treated with full configuration interaction (FCI) [213–215]. While convergence with M is system dependent, as a guide, for a modest number of electrons (10–20), high accuracy can be achieved for more than 100 orbitals, while for larger number of electrons (e.g., up to 40), high accuracy can be achieved for about 40 orbitals, using M up to 10,000.
Furthermore, because the MPS is not built on an excita- tion expansion around a Slater determinant, DMRG calcu- lations are well suited to describe strong or multi-reference correlation, as found in transition metals or excited states [208]. Here, the DMRG is often used to replace a com- plete active space calculation [205,216,217]. Active spaces with up to 40 orbitals can be treated reliably, and the DMRG has been applied to bioinorganic complexes with as many as four transition metal ions, such as the Mn4Ca cluster of pho- tosystem II [218], and [4Fe-4S] clusters. Finally, the MPS mathematically represents one-dimensional chain-like cor- relations very efficiently. The DMRG can thus be used with
great effect in treating correlations in π-systems of many conjugated molecules [216,219,220].
The version of Block included with Q-CHEMcan be run in an entirely black-box fashion; orbital ordering, one of the more unusual inputs into a DMRG calculation, can be determined automatically using a graph-theoretical algo- rithm, or a genetic algorithm optimisation. The user need only specify the final number of states (M) desired. The DMRG module is also completely parallelised; larger cal- culations as described above should be run on 10–100 cores.
3.5. Active space spin-flip methods
This is a family of methods capable of treating strong cor- relations via an active space at lower computational cost than complete active space SCF (CASSCF) type methods, and with greater ease of use. A molecule with strong corre- lations requires multi-configurational wave functions to be even qualitatively correct. For two strongly correlated elec- trons, the first such approach is the spin-flip extended single configuration interaction (SF-XCIS) model [221]. For gen- eral numbers of strongly correlated electrons (though com- putational cost increases exponentially with the number of spin-flips), two implementations [222,223] of the restricted active space spin-flip (RAS-SF) model [224,225] are avail- able. These methods start from a high spin restricted open- shell HF determinant, where the strongly correlated elec- trons are initially all high spin. The target low-spin strongly correlated states are accessed by flipping the spins of half the high spin levels (which define the active space), and per- forming a full CI calculation in the active space, augmented by single excitations into and out of the active space. These additional ‘particle’ (p) and ‘hole’ (h) excitations provide state-specific relaxation of the orbitals.
The efficient implementation of the RAS-SF (i.e., spin-flipping in the RAS-CI formalism) has enabled de- tailed electronic structure studies of singlet fission in dimers plus environment models of pentacene and tetracene crystals [226–229]. Relevant electronic states include delocalised excitonic states, with a variable admixture of charge-resonance configurations, interacting with a dark multi-exciton state of a doubly excited character. Impor- tantly, the RAS-2SF method allows one to treat all electronic states within the same computational framework in dimers and even trimers of relevant compounds (tetracene, pen- tacene, hexacene, etc). RAS-SF calculations enabled inves- tigations of the effect of morphology on the state couplings.
Very recently, the efficiency of the methods has been increased by treating the excitations into and out of the ac- tive space perturbatively, to define the SF-CAS(h,p) method [230]. The basic idea is that if the states of interest are pre- dominantly described by active space configurations, then the small state-specific relaxations that are accounted for by the particle and hole excitations in RAS-SF can be accu- rately approximated by perturbation theory at much lower
Figure 3. H2O potential energy surfaces for double dissocia- tion in the aug-cc-pVTZ basis set. SF-CAS is the dotted line, SF-CAS(h,p)1 is the dashed line, and SF-CAS(S)1 is the solid line. The upper red line is the HF binding energy (no correlation), and the lower red line is the CCSD(T) binding energy (nearly complete treatment of correlation effects at equilibrium vs. disso- ciation).
computational cost. The perturbative framework also per- mits treatment of extended single excitations that go from the hole space to the particle space (i.e., ‘hole–particle’
excitations), as the active space is rearranged, defining the SF-CAS(S) method [231]. SF-CAS(S) is a method that con- tains physics that goes beyond RAS-SF, and therefore be- gins to account for effects that we would normally identify as being associated with dynamic correlation.
As an example of the performance of the newest SF-CAS methods, one may consider the simultaneous bond dissociation of H2O. In the aug-cc-pVTZ basis set, CCSD(T) provides a binding energy of 231 kcal/mol, with 75 kcal/mol attributed to dynamical electron correlation.
In Figure 3, the computed bond dissociation curves for the SF-CAS, SF-CAS(h,p)1, and SF-CAS(S)1methods are compared. In general, a single spin-flip takes care of a single bond dissociation. In this example, two bonds are broken, requiring two spin-flips from the quintet restricted open shell Hartree–Fock (ROHF) reference.
While all the methods smoothly dissociate to correct products without spin-contamination, the binding energy for SF-CAS is significantly less than even that of unre- stricted HF (UHF), reflecting the biasing of quintet orbitals against the bound singlet state. The perturbative treatment
