Features
In this section the available features of the Mrcc code are summarized. We also specify what type of reference states (orbitals) can be used, and if a particular feature requires one of the interfaces or is available with Mrcc in standalone mode. We also give the corresponding references which describe the underlying methodological developments.
Singlepoint energy calculations
Available methods
 1.
conventional and densityfitting (resolutionofidentity) Hartree–Fock SCF (Ref. 129): restricted HF (RHF), unrestricted HF (UHF), and restricted openshell HF (ROHF)
 2.
conventional and densityfitting (resolutionofidentity) multiconfigurational SCF (MCSCF)
 3.
conventional and densityfitting (resolutionofidentity) Kohn–Sham (KS) density functional theory (DFT) (Ref. 67): restricted KS (RKS), unrestricted KS (UKS), restricted openshell KS (ROKS); local density approximation (LDA), generalized gradient approximation (GGA), metaGGA (depending on kineticenergy density and/or the Laplacian of the density), hybrid, rangeseparated hybrid (RSH), and doublehybrid (DH) functionals (for the available functionals see the description of keyword dft); dispersion corrections
 4.
 5.
densityfitting (resolutionofidentity) MP2, spincomponent scaled MP2 (SCSMP2), and scaled oppositespin MP2 (SOSMP2); currently only for RHF and UHF references (Ref. 67)
 6.
densityfitting (resolutionofidentity) randomphase approximation (RPA, also known as ringCCD, rCCD), direct RPA (dRPA, also known as direct ringCCD, drCCD), secondorder screened exchange (SOSEX), renormalized secondorder perturbation theory (rPT2), and approximate RPA with exchange (RPAX2); currently the dRPA, SOSEX, rPT2, and RPAX2 methods are available for RHF/RKS and UHF/UKS references, while the RPA method is only implemented for RHF/RKS (Refs. 67 and 68)
 7.
 8.
densityfitting (resolutionofidentity) secondorder coupledcluster singles and doubles (CC2), configuration interaction singles with perturbative correction for double excitations [CIS(D)], iterative doubles correction to configuration interaction singles [CIS(D${}_{\mathrm{\infty}}$)], and secondorder algebraic diagrammatic construction [ADC(2)] approaches; spincomponent scaled CC2, CIS(D), CIS(D${}_{\mathrm{\infty}}$), and ADC(2) [SCSCC2, SCSCIS(D), SCSCIS(D${}_{\mathrm{\infty}}$), and SCSADC(2)]; scaled oppositespin CC2, CIS(D), CIS(D${}_{\mathrm{\infty}}$), and ADC(2) [SOSCC2, SOSCIS(D), SOSCIS(D${}_{\mathrm{\infty}}$), and SOSADC(2)] (Refs. 87, 88, 84, and 86)
 9.
arbitrary singlereference coupledcluster methods (Ref. 65): CCSD, CCSDT, CCSDTQ, CCSDTQP, …, CC($n$)
 10.
arbitrary singlereference configurationinteraction methods (Ref. 65): CIS, CISD, CISDT, CISDTQ, CISDTQP, …, CI($n$), …, full CI
 11.
arbitrary perturbative coupledcluster models (Refs. 59, 8, and 61):
 •
CCSD[T], CCSDT[Q], CCSDTQ[P], …, CC($n1$)[$n$]
 •
CCSDT[Q]/A, CCSDTQ[P]/A, …, CC($n1$)[$n$]/A
 •
CCSDT[Q]/B, CCSDTQ[P]/B, …, CC($n1$)[$n$]/B
 •
CCSD(T), CCSDT(Q), CCSDTQ(P), …, CC($n1$)($n$)
 •
CCSDT(Q)/A, CCSDTQ(P)/A, …, CC($n1$)($n$)/A
 •
CCSDT(Q)/B, CCSDTQ(P)/B, …, CC($n1$)($n$)/B
 •
CCSD(T)${}_{\mathrm{\Lambda}}$, CCSDT(Q)${}_{\mathrm{\Lambda}}$, CCSDTQ(P)${}_{\mathrm{\Lambda}}$, …, CC($n1$)($n$)${}_{\mathrm{\Lambda}}$
 •
CCSDT1a, CCSDTQ1a, CCSDTQP1a, …, CC($n$)1a
 •
CCSDT1b, CCSDTQ1b, CCSDTQP1b, …, CC($n$)1b
 •
CC2, CC3, CC4, CC5, …, CC$n$
 •
CCSDT3, CCSDTQ3, CCSDTQP3, …, CC($n$)3
 •
 12.
multireference CI approaches (Ref. 66)
 13.
multireference CC approaches using a stateselective ansatz (Ref. 66)
 14.
arbitrary singlereference linearresponse (equationofmotion, EOM) CC methods (Ref. 58): LRCCSD (EOMCCSD), LRCCSDT (EOMCCSDT), LRCCSDTQ (EOMCCSDTQ), LRCCSDTQP (EOMCCSDTQP), …, LRCC($n$) [EOMCC($n$)]
 15.
linearresponse (equationofmotion) MRCC schemes (Ref. 58)
 16.
DFT/WFT embedding (Ref. 46): DFTinDFT, WFTinDFT, WFTinWFT, and WFTinWFTinDFT embedding (where WFT stands for wave function theory); currently embedding is only available for closedshell systems using densityfitting; LDA, GGA, metaGGA, or hybrid functionals can be used as the DFT method; any WFT method implemented in Mrcc can by used in WFTinDFT type embedding calculations; WFTinWFT multilevel methods (via keyword corembed) is only available for local correlation methods.
 17.
The CI and CC approaches listed above are available with RHF, UHF, standard/semicanonical ROHF, MCSCF orbitals. Densityfitting with CI or highorder CC is currently enabled only for RHF references. For the perturbative CC approaches with ROHF reference determinant, for theoretical reasons, semicanonical orbitals are used (see Ref. 61).
Features available via interfaces
 1.
The CI and CC approaches listed above are also available with the following interfaces and references.
 RHF:

Cfour, Columbus, and Molpro
 ROHF, standard orbitals:

Cfour, Columbus, and Molpro
 ROHF, semicanonical orbitals:

Cfour
 UHF:

Cfour, and Molpro
 MCSCF:

Columbus and Molpro
 2.
QM/MM calculations can be performed by the Amber interface using any method implemented in Mrcc for the QM region.
Notes
 1.
Singlepoint calculations are also possible with several types of relativistic Hamiltonians and reference functions, see Sect. 6.7 for more details.
 2.
 3.
 4.
CC$n$ calculations with ROHF orbitals are not possible for theoretical reasons, see Ref. 61 for explanation.
 5.
Singlepoint CI and CC calculations are, in principle, possible with RKS, UKS, ROKS orbitals.
 6.
Solvation effects can be modeled by the polarizable continuum model (PCM) [138] via an interface to the PCMSolver library [16, 17, 1]. The PCM treatment is selfconsistent for HF and KS SCF calculations, while, in the postSCF steps, the potential of the solvent optimized at the SCF level is frozen. See the descriptions of keywords pcm*.
Geometry optimizations and firstorder properties
Available methods
Geometry optimizations and firstorder property calculations can be performed using analytic gradient techniques with the following methods.
 1.
conventional and DF (RI) HFSCF (Ref. 129): RHF and UHF
 2.
conventional and DF (RI) DFT (Ref. 67): RKS and UKS with LDA, GGA, metaGGA (depending only on kineticenergy density), and hybrid functionals as well as dispersion corrections
 3.
double hybrid density functional methods (Ref. 67), such as B2PLYP, B2PLYPD3, B2GPPLYP, etc. (current limitations: only MP2 correlation, closed shell RKS, no spincomponent scaling, no metaGGA functionals, only DH functionals for which the DFT contribution to the energy is stationary with respect to the variation of the MO coefficients)
 4.
DFMP2 (RIMP2), currently only for RHF references (Ref. 67)
 5.
 6.
 7.
 8.
 9.
 10.
 11.
Currently only unrestricted geometry optimizations are possible, and electric dipole, quadrupole, and octapole moments as well as the electric field at the atomic centers can be evaluated. In addition, Mulliken, Löwdin, and intrinsic atomic orbital (IAO) atomic charges, and Mayer bond orders can be computed using the SCF wave functions. Analytic gradients for the CI and CC methods listed above are available with RHF, UHF, and standard ROHF orbitals without density fitting.
The following keywords are available to control the optimization process
 optalg

– to select an algorithm for the optimization
 optmaxit

– maximum number of iterations allowed
 optetol

– convergence criterion for energy change
 optgtol

– convergence criterion for the gradient change
 optstol

– convergence criterion for the stepsize
The optimization will be terminated and regarded as successful when the maximum gradient component becomes less than optgtol and either an energy change from the previous step is less than optetol or the maximum displacement from the previous step is less than optstol. For their detailed description see Sect. 12.
Features available via interfaces
 1.
The implemented analytic gradients for the CI and CC approaches listed above can also be utilized via the Cfour and Columbus interfaces with the following references.
 RHF:

Cfour and Columbus
 ROHF, standard orbitals:

Cfour and Columbus
 UHF:

Cfour
 MCSCF:

Columbus
In addition to geometries, most of the firstorder properties (dipole moments, quadrupole moments, electric field gradients, relativistic contributions, etc.) implemented in Cfour and Columbus can be calculated with Mrcc.
 2.
QM/MM geometry optimizations and MD calculations can be performed by the Amber interface using any method implemented in Mrcc for which analytic gradients are available.
Notes
 1.
Geometry optimizations and firstorder property calculations can also be performed via numerical differentiation for all methods available in Mrcc using the Cfour interface.
 2.
Analytic gradients are are also available with several types of relativistic Hamiltonians and reference functions, see Sect. 6.7 for more details.
 3.
Harmonic frequencies and secondorder properties
Available methods
Harmonic vibrational frequencies, infrared (IR) intensities, and ideal gas thermodynamic properties can be evaluated using numerically differentiated analytic gradients for all the methods listed in Sect. 6.2.
Features available via interfaces
CC and CI harmonic frequency and secondorder property calculations for RHF and UHF references can also be performed using analytic second derivatives (linear response functions) with the aid of the Cfour interface. Analytic Hessians (LR functions) are available for the following approaches.
 1.
 2.
 3.
 4.
In addition to harmonic vibrational frequencies [57], the analytic Hessian code has been tested for NMR chemical shifts [57], static and frequencydependent electric dipole polarizabilities [60], magnetizabilities and rotational $g$tensors [30], electronic $g$tensors [29], spinspin coupling constants, and spin rotation constants. These properties are available via the Cfour interface.
Notes
 1.
Using the Cfour interface harmonic frequency calculations are also possible via numerical differentiation of energies for all implemented methods with RHF, ROHF, and UHF orbitals.
 2.
Using the Cfour or the Columbus interface harmonic frequency calculations are also possible via numerical differentiation of analytic gradients for all implemented methods for which analytic gradients are available (see Sect. 6.2 for a list of these methods). With Cfour the calculation of static polarizabilities is also possible using numerical differentiation.
 3.
NMR chemical shifts can be computed for closedshell molecules using gaugeincluding atomic orbitals and RHF reference function.
Higherorder properties
Features available via interfaces
Thirdorder property calculations can be performed using analytic third derivative techniques (quadratic response functions) invoking the Cfour interface for the following methods with RHF and UHF orbitals.
Notes
 1.
The analytic third derivative code has been tested for static and frequencydependent electricdipole first (general, secondharmonicgeneration, opticalrectification) hyperpolarizabilities [107] and Raman intensities [106]. Please note that the orbital relaxation effects are not considered for the electricfield. These properties are available via the Cfour interface.
 2.
Using the Cfour interface anharmonic force fields and the corresponding spectroscopic properties can be computed using numerical differentiation techniques together with analytic first and/or analytic second derivatives at all computational levels for which these derivatives are available (see Sect. 6.2 and 6.3 for a list of these methods).
Diagonal Born–Oppenheimer corrections
Features available via interfaces
Diagonal Born–Oppenheimer correction (DBOC) calculations can be performed using analytic second derivative techniques via the Cfour interface for the following methods with RHF and UHF references.
 1.
 2.
 3.
 4.
Electronically excited states
Available methods
Excitation energies, firstorder excitedstate properties, and ground to excitedstate transition moments can be computed as well as excitedstate geometry optimizations can be performed using linear response theory and analytic gradients with the following methods.
 1.
 2.
 3.
 4.
Excitation energy and property calculations for the aforementioned methods are available with RHF, UHF, and standard ROHF orbitals. Densityfitting is only possible for RHFbased single point calculations. So far electric and magnetic dipole transition moments, both in the length and the velocity gauge, as well as the corresponding oscillator and rotator strengths have been implemented. For the list of implemented firstorder properties see Sect. 6.2.
Excitation energies can also be computed for closedshell systems using the densityfitting approximation and RHF (RKS) orbitals with the following methods (Refs. 87, 88, 85, and 84).
 1.
CIS, timedependent HF (TDHF), Tamm–Dancoff approximation (TDA), timedependent DFT (TDDFT)
 2.
CIS(D) and ADC(2)based double hybrid TDDFT methods.
 3.
secondorder coupledcluster singles and doubles (CC2) method
 4.
configuration interaction singles with perturbative correction for double excitations [CIS(D)]
 5.
iterative doubles correction to configuration interaction singles [CIS(D${}_{\mathrm{\infty}}$)] method
 6.
secondorder algebraic diagrammatic construction [ADC(2)] approach
 7.
spinscaled versions of the latter approaches: SCSCC2, SCSCIS(D), SCSCIS(D${}_{\mathrm{\infty}}$), SCSADC(2), SOSCC2, SOSCIS(D), SOSCIS(D${}_{\mathrm{\infty}}$), and SOSADC(2)
Analytic gradients are implemented for CIS and TDHF. Ground to excitedstate transition moments are available for TDHF, TDA, TDDFT, double hybrid TDDFT, CIS(D), SCSCIS(D), SOSCIS(D), ADC(2), SCSADC(2), and SOSADC(2). For the CIS, TDHF, TDA, TDDFT, CC2, CIS(D${}_{\mathrm{\infty}}$), ADC(2), SCSCC2, SCSCIS(D${}_{\mathrm{\infty}}$), SCSADC(2), SOSCC2, SOSCIS(D${}_{\mathrm{\infty}}$), and SOSADC(2) methods efficient reducedcost approaches are also implemented, see Sect. 6.8 for details.
Features available via interfaces
Excitation energies, firstorder excitedstate properties, and ground to excitedstate transition moments can also be calculated as well as excitedstate geometry optimizations can also be carried out using the following interfaces and reference states.
 RHF:

Cfour, Columbus, and Molpro (only excitation energy)
 ROHF, standard orbitals:

Cfour, Columbus, and Molpro (only excitation energy)
 UHF:

Cfour and Molpro (only excitation energy)
 MCSCF:

Columbus and Molpro (only excitation energy)
Notes
 1.
Please note that for excitation energies and geometries LRCC methods are equivalent to the corresponding EOMCC models. It is not true for firstorder properties and transition moments.
 2.
With CI methods excited to excitedstate transition moments can also be evaluated.
 3.
Excitedstate harmonic frequencies can be evaluated for the above methods with the help of numerical differentiation of analytical gradients, see Sect. 6.3.
 4.
Excitedstate harmonic frequencies can also be calculated for the above methods via numerical differentiation using the Cfour or Columbus interface.
 5.
Excitedstate harmonic frequencies and secondorder properties can be evaluated for CI methods using analytic second derivatives and the Cfour interface.
Relativistic calculations
Treatment of special relativity in singlepoint energy calculations is possible for all the CC and CI methods listed in Sect. 6.1 using various relativistic Hamiltonians with the following interfaces (Refs. 103 and 63).
 1.
With Molpro relativistic calculations can be performed with Douglas–Kroll–Hess Hamiltonians using RHF, UHF, ROHF, and MCSCF orbitals. The interface also enables the use of effective core potentials (see Molpro’s manual for the specification of the Hamiltonian and effective core potentials).
 2.
With Cfour exact twocomponent (X2C) and spinfree Dirac–Coulomb (SFDC) calculations can be performed. The evaluation of massvelocity and Darwin corrections is also possible using analytic gradients for all the methods and reference functions listed in Sect. 6.2. (See the description of the RELATIVISTIC keyword in the Cfour manual for the specification of the Hamiltonian.)
 3.
Treatment of special relativity in analytic gradient calculations is possible for all the CC and CI methods listed in Sect. 6.2 using various relativistic Hamiltonians with the following interfaces.
 1.
With Cfour analytic gradient calculations can be performed with the exact twocomponent (X2C) treatment.
 2.
Reducedscaling and local correlation calculations
Orbital transformation techniques
The computational expenses of the CC and CI methods listed in Sect. 6.1 can be reduced via orbital transformation techniques (Ref. 127). In this framework, to reduce the computation time the dimension of the properly transformed virtual oneparticle space is truncated. Currently optimized virtual orbitals (OVOs) or MP2 natural orbitals (NOs) can be chosen. This technique is recommended for small to mediumsize molecules. This scaling reduction approach is available using RHF or UHF orbitals. See the description of keywords ovirt, eps, and ovosnorb for more details.
Natural auxiliary functions
The cost of densityfitting methods can be reduced using natural auxiliary functions (NAFs) introduced in Ref. 67. The approach is very efficient for dRPA, but considerable speedups can also be achieved for MP2, CC2, and ADC(2). See the description of keywords naf_cor and naf_scf for more details.
Reducedcost techniques for excited states
The computational expenses of CIS, TDHF, TDA, TDDFT, CC2, CIS(D${}_{\mathrm{\infty}}$), ADC(2), SCSCC2, SCSCIS(D${}_{\mathrm{\infty}}$), SCSADC(2), SOSCC2, SOSCIS(D${}_{\mathrm{\infty}}$), and SOSADC(2) excitedstate calculations can be efficiently reduced using local fitting domains as well as stateaveraged NOs and NAFs (Refs. 87, 88, and 85). The scaling of CC2 and ADC(2) calculations can also be decreased by local correlation approaches (Ref. 84). See the description of keywords redcost_exc and redcost_tddft for more details.
Local correlation methods
The cost of MP2, dRPA, SOSEX as well as singlereference iterative and perturbative coupledcluster calculations can be reduced for large molecules by the local natural orbital CC (LNOCC) approach (Refs. 126, 129, 68, 101, 99, 102, and 100). This method combines ideas from the clusterinmolecule approach of Li and coworkers [76], the incremental approach of Stoll et al. [135], domain and pair approximations introduced first by Pulay et al. (see, e.g., Ref. 121) with frozen natural orbital, natural auxiliary function, and Laplace transform techniques. It is currently available only for closedshell molecules using RHF (RKS) orbitals. See the description of keywords localcc, lnoepso, lnoepsv, domrad, lmp2dens, dendec, nchol, osveps, spairtol, wpairtol, laptol, lccrest, and lcorthr for further details. Extensive benchmarks regarding the accuracy and efficiency of the local correlation methods are also provided in Refs. 68, 101, 102, and 100 and Section II.G of Ref. 62.
Multilevel local correlation methods
Utilizing the above local correlation techniques a multilevel scheme is defined in which the LMOs are classified as active or environment (Ref. 46, 47). The contributions of these LMOs to the total correlation energy are evaluated using different models for the two subsystems, for instance, one can choose a LNOCC model for the active subsystem and LMP2 for the environment. See the description of keyword corembed for further details.
Optimization of basis sets
The optimization of basis set’s exponents and contraction coefficients can be performed with any method for which singlepoint energy calculations are available (see Sect. 6.1). The implementation is presented in Ref. 83. The related keywords are
 basopt

– to turn on/off basis set optimization
 optalg

– to select an algorithm for the optimization
 optmaxit

– maximum number of iterations allowed
 optetol

– convergence criterion for energy change
 optstol

– convergence criterion for parameter (exponent, contraction coefficient) change
For their detailed description see Sect. 12.
For the optimization of basis sets it is important to know the format for the storage of the basis set parameters. In Mrcc the format used by the Cfour package is adapted. The format is communicated by the following example.
actual lines  description  

C:631G  $\hookleftarrow $ Carbon atom:basis name  
Pople’s Gaussian basis set  $\hookleftarrow $ comment line  
$\hookleftarrow $ blank line  
2  $\hookleftarrow $ number of angular momentum types  
0  1  $\hookleftarrow $ 0$\to $s , 1$\to $p  
3  2  $\hookleftarrow $ number of contracted functions  
10  4  $\hookleftarrow $ number of primitives  
$\hookleftarrow $ blank line  
3047.5249  457.36952  …  $\hookleftarrow $ exponents for s functions  
$\hookleftarrow $ blank line  
0.0018347  0.0000000  0.0000000  $\hookleftarrow $ contraction coefficients  
0.0140373  0.0000000  0.0000000  for s functions  
0.0688426  0.0000000  0.0000000  
0.2321844  0.0000000  0.0000000  
0.4679413  0.0000000  0.0000000  
0.3623120  0.0000000  0.0000000  
0.0000000  0.1193324  0.0000000  
0.0000000  0.1608542  0.0000000  
0.0000000  1.1434564  0.0000000  
0.0000000  0.0000000  1.0000000  
$\hookleftarrow $ blank line  
7.8682724  1.8812885  …  $\hookleftarrow $ exponents for p functions  
$\hookleftarrow $ blank line  
0.0689991  0.0000000  $\hookleftarrow $ contraction coefficients  
0.3164240  0.0000000  for p functions  
0.7443083  0.0000000  
0.0000000  1.0000000 
In a basis set optimization process you need two files in the working directory: the appropriate MINP file with the basopt keyword set and a user supplied GENBAS file that contains the basis set information in the above format. You do not need to write the GENBAS file from scratch, you can use the files in the BASIS directory of Mrcc to generate one or you can use the Basis Set Exchange [5, 22, 133, 120] to download a basis in the appropriate form (CFOUR format). Note that you can optimize several basis sets at a time: all the basis sets which are added to the GENBAS file will be optimized.
You can perform unconstrained optimization when all the exponents and contraction coefficients are optimized except the ones which are exactly 0.0 or 1.0. Alternatively, you can run constrained optimizations when particular exponents/coefficients or all exponents and coefficients for a given angular momentum quantum number are kept fixed during the optimization. The parameters to be optimized can be specified in the GENBAS file as follows.
 1.
Unconstrained optimization: no modifications are needed—by default all exponents and contraction coefficients will be optimized except the ones which are exactly 0.0 or 1.0.
 2.
Constrained optimization: by default all the exponents and coefficients will be optimized just as for the unconstrained optimization. To optimize/freeze particular exponents or coefficients special marks should be used:
 •
use the “” mark (without quotes) if you want to keep an exponent or coefficient fixed during the optimization. You should put this mark right after the fixed parameter (no blank space is allowed). If this mark is attached to an angular momentum quantum number, none of the exponents/coefficients of the functions in the given shell will be optimized except the ones which are marked by “++”.
 •
use the “++” mark (without quotes) if you want a parameter to be optimized. Then you should put this mark right after it (no blanks are allowed). You might wonder why this is needed if the default behavior is optimization. Well, this makes life easier. If you want to optimize just a few parameters, it is easier to constrain all parameters first then mark those, which are needed to be optimized (see the example below).
 •
Examples:
 1.

To reoptimize all parameters in the above basis set but the exponents and coefficients of stype functions you should copy the basis set to the GENBAS file and put mark “” after the angular momentum quantum number of 0. The first lines of the GENBAS file:
C:631G Pople’s Gaussian basis set 2 0 1 3 2 10 4 3047.5249 457.36952 …  2.

Both s and ptype functions are fixed but the first sexponent:
C:631G Pople’s Gaussian basis set 2 0 1 3 2 10 4 3047.5249++ 457.36952 …
During the optimization the GENBAS file is continuously updated, and if the optimization terminated successfully, it will contain the optimized values (in this case it is equivalent to the GENBAS.opt file, see below, the only difference is that the file GENBAS.opt may contain the special marks, i.e., “++”, “”). Further files generated in the optimization are:
 •
GENBAS.init – the initial GENBAS file saved
 •
GENBAS.tmp – temporary file, updated after each iteration, can be used to restart conveniently a failed optimization process
 •
GENBAS.opt – this file contains the optimized parameters after a successful optimization.