10.3 First-Order Properties and Gradients

For the ground state first-order properties (expectation values) are implemented at the SCF, MP2 and CC2 level. Note that for the ground state CCS and CIS are equivalent to SCF. For excited states first-order properties are implemented only at the CCS and CC2 level. Gradients are presently only available for the ground state at the MP2 and the CC2 and for excited states only at the CC2 level.

10.3.1 Ground State Properties, Gradients and Geometries

For CC2, one distinguishes between orbital-relaxed and unrelaxed properties. Both are calculated as first derivatives of the respective energy with respect to an external field corresponding to the calculated property. They differ in the treatment of the SCF orbitals. In the orbital-relaxed case the external field is (formally) already included at the SCF stage and the orbitals are allowed to relax in the external field; in the orbital-unrelaxed case the external field is first applied after the SCF calculation and the orbitals do not respond to the external field. Orbital-unrelaxed CC2 properties are calculated as first derivatives of the real part of the unrelaxed Lagrangian [133]

Lur CC2(t,t,β) =  ⟨HF |H |CC ⟩+    tμ1⟨μ1|ˆH + [ˆH,T2]|HF ⟩      (10.12)
                   ∑       ˆ         ˆ
                 + μ  tμ2⟨μ2|H + [F0 + βV,T2]|HF ⟩
with H = H0 + βV —where V is the (one-electron) operator describing the external field, β the field strength, and H0 and F0 are the Hamiltonian and Fock operators of the unperturbed system—by the expression:
   ur CC2      ( ∂Lur CC2(t,t,β))    ∑    ur
⟨V⟩       =  ℜ   -----∂β-------  =     D pq Vpq ,         (10.13)
               (               0     pq
          =  ℜ  ⟨HF|ˆV |HF ⟩+ ∑  t ⟨μ |ˆV +[V,T ]|HF ⟩       (10.14)
                            μ1  μ1  1        2
               ∑                 )
             +    tμ2⟨μ2|[Vˆ,T2]|HF⟩   ,
where indicates that the real part is taken. Relaxed CC2 properties (and gradients) are calculated from the the full variational density including the contributions from the orbital response to the external perturbation, which are derived from the Lagrangian [11,137]
Lrel CC2(t,t) =  ⟨HF|H |CC ⟩+    tμ1⟨μ1|ˆH + [Hˆ, T2]|HF ⟩        (10.15)
                  ∑       ˆ              ∑
                + μ  tμ2⟨μ2|H + [F,T2]|HF ⟩+  μ κμ0Fμ0 ,
                   2                       0
where F is the Fock operator corresponding to the Hamiltonian of the perturbed system H = H0 + βV . One-electron properties are then obtained as:
⟨V⟩rel CC2 =  ℜ  ⟨HF |ˆV |HF ⟩+ ∑  t ⟨μ |ˆV + [V,T ]|HF ⟩       (10.16)
                            μ   μ1 1        2
               ∑             1     ∑        )
             +    tμ2⟨μ2|[V,T2]|HF⟩+    κ μ0Vμ0   ,
                μ2                  μ0
             ∑    rel
          =   pq D pq Vpq .                               (10.17)
The calculation of one-electron first-order properties requires that in addition to the cluster equations also the linear equations for the Lagrangian multipliers tμ are solved, which requires similar resources (CPU, disk space, and memory) as the calculation of a single excitation energy. For orbital-relaxed properties also a CPHF-like linear equation for the Lagrangian multipliers κμ0 needs to be solved and the two-electron density has to be build, since it is needed to set up the inhomogeneity (right-hand side). The calculation of relaxed properties is therefore somewhat more expensive—the operation count for solving the so-called Z-vector equations is similar to what is needed for an SCF calculation—and requires also more disk space to keep intermediates for the two-electron density—about O(2V + 2N)Nx + Nx2 in addition to what is needed for the solution of the cluster equations. For ground states, orbital-relaxed first-order properties are standard in the literature.

The calculation of the gradient implies the calculation of the same variational densities as needed for relaxed one-electron properties and the solution of the same equations. The construction of the gradient contributions from the densities and derivative integrals takes about the same CPU time as 3–4 SCF iterations and only minor extra disk space. For details of the implementation of CC2 relaxed first-order properties and gradients and a discussion of applicability and trends of CC2 ground-state equilibrium geometries see ref.  [11]. The following is in example input for a MP2 and CC2 single point calculation of first-order properties and gradients:

  static relaxed operators=diplen,qudlen  

A different input is required for geometry optimizations: in this case the model for which the geometry should be optimized must be specified in the data group $ricc2 by the keyword geoopt:

  geoopt model=cc2

For CC2 calculations, the single-substitution part of the Lagrangian multipliers tμ are saved in the file CCL0--1--1---0 and can be kept for a restart (for MP2 and CCS, the single-substitution part tμ vanishes).

For MP2 only relaxed first-order properties and gradients are implemented (unrelaxed MP2 properties are defined differently than in CC response theory and are not implemented). For MP2, only the CPHF-like Z-vector equations for κμ0 need to be solved, no equations have to be solved for the Lagrangian multipliers tμ. CPU time and disk space requirements are thus somewhat smaller than for CC2 properties or gradients.

For SCF/CIS/CCS it is recommended to use the modules grad and rdgrad for the calculation of, ground state gradients and first-order properties.

10.3.2 Excited State Properties, Gradients and Geometries

Also for excited states presently unrelaxed and relaxed first-order properties are available in the ricc2 program. These are implemented for CCS and CC2. Note, that in the unrelaxed case CIS and CCS are not equivalent for excited-states first-order properties and no first-order properties are implemented for CIS in the ricc2 program.

Orbital-unrelaxed first-order properties

The unrelaxed first-order properties are calculated from the variational excited states Lagrangian [138], which for the calculation of unrelaxed properties is composed of the unrelaxed ground state Lagrangian, Eq. (10.12), and the expression for the excitation energy:

 ur CC2, ex   (ex)                    ∑   
L       (E, E,t,t   ,β)  =  ⟨HF|H |CC ⟩+    E μAμν(t,β)Eν       (10.18)
                             ∑          μν
                           +    t(eμx1) ⟨μ1|ˆH + [Hˆ, T2]|HF ⟩
                           + ∑  t(ex)⟨μ2|ˆH + [F0 + βˆV,T2]|HF ⟩
                             μ2  μ2
where it is assumed that the left and right eigenvectors are normalized such that μνEμμ|τνEν = 1 and H = H0 + βV . The first-order properties are calculated as first derivatives of Lur CC2, ex(E,E,t,t(ex)) with respect to the field strength β and are evaluated via a density formalism:
   ur,ex       ( ∂Lur,ex(E, E,t,t(ex),β))    ∑    ur,ex
⟨V ⟩     =  R   --------∂β----------  =    D pq Vpq ,      (10.19)
                                     0   pq
(Again indicates that the real part is taken.) The unrelaxed excited-state properties obtained thereby are related in the same way to the total energy of the excited states as the unrelaxed ground-state properties to the energy of the ground state and the differences between excited- and ground-state unrelaxed properties are identical to those identified from the second residues of the quadratic response function. For a detailed description of the theory see refs.  [137,138]; the algorithms for the RI-CC2 implementation are described in refs.  [10,135]. ref.  [135] also contains a discussion of the basis set effects and the errors introduced by the RI approximation.

The calculation of excited-state first-order properties thus requires the calculation of both the right (Eμ) and left (Eμ) eigenvectors and of the excited state Lagrangian multipliers tμ(ex). The disk space and CPU requirements for solving the equations for Eμ and tμ(ex) are about the same as those for the calculation of the excitation energies. For the construction of the density matrices in addition some files with O(nrootN2) size are written, where nroot is the number of excited states.

The single-substitution parts of the excited-states Lagrangian multipliers tμ(ex) are saved in files named CCNL0-s--m-xxx.

For the calculation of first-order properties for excited states, the keyword exprop must be added with appropriate options to the data group $excitations; else the input is same as for the calculation of excitation energies:

  fop unrelaxed_only operators=diplen,qudlen  
  irrep=a1 nexc=2  
  exprop states=all operators=diplen,qudlen

Orbital-relaxed first-order properties and gradients

To obtain orbital-relaxed first-order properties or analytic derivatives (gradients) the Lagrange functional for the excited state in Eq. (10.18) is—analogously to the treatment of ground states—augmented by the equations for the SCF orbitals and the perturbations is also included in the Fock operator:

Lrel CC2, ex(E,E, t,t(ex),β) = ⟨HF|H |CC ⟩+    EμAμν(t,β)Eν           (10.20)
                            + ∑  t(ex)⟨μ1|ˆH + [Hˆ, T2]|HF ⟩
                              μ1  μ1
                              ∑   (ex)                  ∑   (ex)
                            +    tμ2 ⟨μ2|ˆH + [F,T2]|HF ⟩+    κμ0 Fμ0 .
                              μ2                        μ0
Compared to unrelaxed properties, the calculation of relaxed properties needs in addition for each excited state the solution of a CPHF equations for the Lagrangian multipliers κμ0(ex), for which the computational costs are similar to those of a Hartree-Fock calculation.

Orbital-relaxed properties are requested by adding the flag relaxed to the input line for the exprop option. The following is an example for a CC2 single point calculation for orbital-relaxed excited state properties:

  irrep=a1 nexc=2  
  exprop states=all relaxed operators=diplen,qudlen

Note that during the calculation of orbital-relaxed excited-state properties the corresponding unrelaxed properties are also automatically evaluated at essentially no additional costs. Therefore, the calculation of unrelaxed properties can not be switched off when relaxed properties have been requested.

Again the construction of gradients requires the same variational densities as needed for relaxed one-electron properties and the solution of the same equations. The construction of the gradient contributions from one- and two-electron densities and derivative integrals takes approximately the same time as for ground states gradients (approx. 3–4 SCF iterations) and only minor extra disk space. The implementation of the excited state gradients for the RI-CC2 approach is described in detail in Ref.  [12]. There one can also find some information about the performance of CC2 for structures and vibrational frequencies of excited states.

For the calculation of an excited state gradient with CC2 at a single point (without geometry optimization and if it is not a calculation with NumForce) one can use the input:

  irrep=a1 nexc=2  
  xgrad states=(a1 1-2)

For geometry optimizations or a numerical calculation of the Hessian with NumForce the wavefunction model and the excited state for which the geometry should be optimized have to be specified in the data group $ricc2 with the keyword geoopt:

  geoopt model=cc2 state=(a1 2)  
  irrep=a1 nexc=2

If the geometry optimization should carried out for the lowest excited state (of those for which an excitation energy is requested in $excitation), one can use alternatively state=(s1).

Since the calculation of unrelaxed and relaxed first-order properties can be combined gradient calculations without significant extra costs, a request for excited state gradients will automatically enforce the calculation of the relaxed and unrelaxed dipole moments. If the keyword geoopt is used, the relaxed dipole moment for the specified excited state and wavefunction model will be written to the control file and used in calculations with NumForce for the evaluation of the IR intensities.

10.3.3 Visualization of densities and Density analysis

As most other programs which allow for the calculation of wavefunctions and densities also the ricc2 module is interfaced to wavefunction analysis and visualization toolbox described in chapter 17. From ricc2 module this interface can used in two different ways

If through the geoopt keyword in $ricc2 a unique method and state has been specified for which the density, gradient and properties are evaluated, the density analysis and visualization routines will called by default with the (orbital-relaxed) density for this state and method similar as in dscf, ridft, mpgrad, etc.
The ricc2 program can be called in a special analysis mode which allows to analyse densities and combination (e.g. differences) of densities evaluated in preceeding ricc2 calculations.

Default density analysis and visualization:  
As in a single calculations with the ricc2 program one–electron densities can be calculated for more than one method and/or electronic state, the interface to the analysis and visualization routines require the specification of a unique level of calculation and a unique state. This is presently done through the geoopt flag which determines the method/state for which results are written to interface files (e.g. control, gradient, or xxx.map).

In ground state calculations ricc2 will pass to the density analysis routines the correlated total (and for UHF based calculations also the spin) density and the canonical SCF orbitals from which the SCF (spin) density is constructed. All options described in chapter 17 are available from within the ricc2 program apart from the evaluation of electrostatic moments, which would interfere with the calculation of expectation values requested through the fop option in $response.

In excited state calculation ricc2 will pass the excited state total (and for UHF based calculation in addition the spin) density. But no ground state densities and/or uncorrelated densities or orbitals. Thus, for excited states the ricc2 program does, in difference to egrad not print out a comparison with the ground state SCF density. Also, all some options which require orbitals (as e.g. the generation and visualization of localized orbitals or some population analysis options) and not available for excited states in ricc2.

As other modules, also ricc2 provides the -proper flag to bypass a re-calculation of the density and gradient to enter immediately the density analysis routines with a previously calculated density. The ricc2 program will then pass the densities found on the interface file for the density analysis routines without further check on the method and state for which they have been evaluated. If both, ground and excited state densities are found on file, both will be passed to the density analysis, thereby providing a shortcut to the -fanal and the $anadens keyword for the analysis of differences between ground and excited state densities.

The general density analysis option:  
In general ricc2 saves by default all relaxed densities generated during a calculation in files named cc1td-<type>-<mult><irrep>-<number>, where cc1td stands for “coupled-cluster one-electron total density”. <type> is one of mp2-gs (MP2 ground state), cc2-gs (CC2 ground state), ccs-xs (CCS excited state), cc2-xs (CC2 excited state), or adc2-xs (ADC(2) excited state) and the other entries specify multiplicity, irreducible representation and the number of the state. Having specified the calculation of relaxed densities—e.g. by requesting relaxed one-electron properties or as a by-product of a gradient calculation—you will end up with two files named like


In case of open shell molecules, additional files with names cc1sd... (for one-electron spin-densities) will be generated.

These files are (currently) in a binary format, similar as the files dens, mdens and edens. Therefore be aware that a transfer between different computer architectures may result in trouble.

The densities on these files can be analysed with the tools and interfaces provided by Moloch (see Section 17.2). This can be done by calling ricc2 with the option -fanal which bypasses the usual wavefunction calculation and triggers the program into an analysis mode for densities. In this mode the program interpretes $anadens and the keywords described in Section 17.2. To plot, for example, the difference density of the two above mentioned total densities you have to add the following lines in your control file

 calc my_favourite_diffden from  
 1d0 cc1td-cc2-xs-3a2-001  
-1d0 cc1td-cc2-gs-1a1-001  

and invoke

  ricc2 -fanal

This will generate the files my_favourite_diffden and my_favourite_diffden.map. The latter can be converted into gOpenMol format as described in Section 17.2.

10.3.4 Fast geometry optimizations with RI-SCF based gradients

If geometry optimizations on MP2 or CC2 level are performed with large basis set, especially with diffuse basis functions, the N4–steps might become the dominant part of the overall timings. In these cases, the integral screening in the Hartree–Fock part often becomes inefficient. The resolution–of–the–identity can be applied here to speed up the calculation of the HF reference wavefunction, as well as the solution of the coupled–perturbed Hartree–Fock (CPHF) equations in the MP2 or CC2 gradient calculation.

An additional auxiliary basis (denoted jkbas) set has to be assigned via the General Options Menu in the define program. In the submenu rijk choose on and select your auxiliary basis set. Then, run the jobex script the additional rijk-flag:

> jobex -level cc2 -rijk