# The closed shell CCSD program

Bibliography:

C. Hampel, K. Peterson, and H.-J. Werner, Chem. Phys. Lett. **190**, 1 (1992)

All publications resulting from use of this program must acknowledge the above.

The CCSD program is called by the CISD, CCSD, BCCD, or QCI directives. CID or CCD can be done as special cases using the `NOSINGL`

directive. The code also allows to calculate Brueckner orbitals (QCI and CCSD are identical in this case). Normally, no further input is needed if the `CCSD`

card follows the corresponding HF-SCF. Optional `ORBITAL`

, `OCC`

, `CLOSED`

, `CORE`

, `SAVE`

, `START`

, `PRINT`

options work as described for the MRCI program in section the MRCI program. The only special input directives for this code are `BRUECKNER`

and `DIIS`

, as described below.

The following options may be specified on the command line:

Ignore convergence checks.`NOCHECK`

Do calculation integral direct.`DIRECT`

Do not include singly external configurations.`NOSING`

Maximum number of iterations.`MAXIT=`

*value*Denominator shift for update of singles.`SHIFTS=`

*value*Denominator shift for update of doubles.`SHIFTP=`

*value*Convergence threshold for the energy.`THRDEN=`

*value*Convergence threshold for CC amplitudes. This applies to the square sum of the changes of the amplitudes.`THRVAR=`

*value*Save CCSD wavefunction to this record`SAVE=`

*record*Restart CCSD wavefunction from a previously written record`START=`

*record*Specifies the zeroth-order Hamiltonian in closed-shell CCSD calculations with KS orbitals.`HO=`

*operator*`H0=MP`

: Moeller Plesset partitioning (default); the Fock matrix is computed and block diagonalised. The eigenvalues are used in the MP2 and (T) denominators.`H0=KS`

: The KS Fock matrix is used as zeroth-order Hamiltonian, and its eigenvalues are used in the denominators.Forces the use of the KS Fock operator in the subsequent calculation. This is required for double-hybrid and RPA calculations.`KSFOCK`

The convergence thresholds can also be modified using

`THRESH`

,`ENERGY=`

*thrden*,`COEFF=`

*thrvar*

Convergence is reached if the energy change is smaller than *thrden* (default 1.d-6) *and* the square sum of the amplitude changes is smaller than *thrvar* (default (1.d-10). The `THRESH`

card must follow the command for the method (e.g., `CCSD`

) and then overwrites the corresponding global options (see `GTHRESH`

, sec. global Thresholds (GTHRESH)).

The computed energies are stored in variables as explained in section special variables. As well as the energy, the $T_1$ diagnostic (T. J. Lee and P. R. Taylor, Int. J. Quant. Chem. S23 (1989) 199) and the $D_1$ diagnostic (C. L. Janssen and I. M. B. Nielsen, Chem. Phys. Lett. 290 (1998), 423, and T. J. Lee. Chem. Phys. Lett. 372 (2003), 362) are printed and stored for later analysis in the variables T1DIAG and D1DIAG, respectively.

For geometry optimization and computing first-order properties see section expectation values.

## Coupled-cluster, CCSD

The command `CCSD`

performs a closed-shell coupled-cluster calculation. Using the `CCSD(T)`

command, the perturbative contributions of connected triple excitations are also computed.

If the CCSD is not converged, an error exit will occur if triples are requested. This can be avoided using the `NOCHECK`

option:

`CCSD(T),NOCHECK`

In this case the (T) correction will be computed even if the CCSD did not converge. Note: `NOCHECK`

has no effect in geometry optimizations or frequency calculations.

For further information on triples corrections see under RCCSD.

In cases that the Hartree-Fock reference determinant does not sufficiently dominate the wavefunction, the program will stop with a message `UNREASONABLE NORM. CALCULATION STOPPED`

. This can be avoided by setting the option `CC_NORM_MAX`

=*value* on the CCSD command line. The calculation is stopped if the square norm of the wavefunction (in intermediate normalization) exceeds `1+nval*CC_NORM_MAX`

, where nval is the number of correlated valence orbitals. The default value is `CC_NORM_MAX=0.2`

.

CCSD calculations with KS orbitals are possible. Optionally, the `H0`

option can be used to specify the zeroth-order Hamiltonian (see options above). Note that some singles terms are still missing in the closed-shell (T) program in case that non-HF orbitals are used. This will be fixed as soon as possible.

## Quadratic configuration interaction, QCI

`QCI`

or `QCISD`

performs quadratic configuration interaction, QCISD. Using the `QCI(T)`

or `QCISD(T)`

commands, the contributions of connected triples are also computed by perturbation theory. Normally, no further input is needed if the `QCI`

card follows the corresponding HF-SCF. Otherwise, occupancies and orbitals can be specified as in the CI program. For modifying DIIS directives, see section the DIIS directive.

For avoiding error exits in case of no convergence, see `CCSD(T)`

.

## Brueckner coupled-cluster calculations, BCCD

`BCCD`

,[`SAVEORB=`

*record*],[`PRINT`

],[`TYPE=`

,*type*]

In addition, the same options as for CCSD can be used.

`BCCD`

performs a Brueckner coupled-cluster calculation and computes Brueckner orbitals. With these orbitals, the amplitudes of the singles vanish at convergence. Using the `BCCD(T)`

command, the contributions of connected triples are also computed by perturbation theory. Normally, no further input is needed if the `BCCD`

card follows the corresponding HF-SCF. Otherwise, occupancies and orbitals can be specified as in the CI program. `BRUECKNER`

parameters can be modified using the `BRUECKNER`

directive.

The Brueckner orbitals and approximate density matrix can be saved on a `MOLPRO`

dump record using the `SAVEORB`

option. The orbitals are printed if the `PRINT`

option is given. `TYPE`

can be used to specify the type of the approximate density to be computed:

Compute and store density of reference determinant only (default). This corresponds to the BOX (Brueckner orbital expectation value) method of Chem. Phys. Lett.`TYPE=REF`

**315**, 248 (1999).Compute and store density with contribution of pair amplitudes (linear terms). Normally, this does not seem to lead to an improvement.`TYPE=TOT`

Compute and store both densities`TYPE=ALL`

Note: The expectation variables are stored in variables as usual. In the case that both densities are made, the variables contain two values, the first corresponding to `REF`

and the second to `TOT`

(e.g., DMZ(1) and DMZ(2)). If `TYPE=REF`

or `TYPE=TOT`

is give, only the corresponding values are stored.

For avoiding error exits in case of no convergence, see `CCSD(T)`

.

### The BRUECKNER directive

`BRUECKNER`

,*orbbrk,ibrstr,ibrueck,brsfak*;

This directive allows the modification of options for Brueckner calculations. Normally, none of the options has to be specified, and the `BCCD`

command can be used to perform a Brueckner CCD calculation.

if nonzero, the Brueckner orbitals are saved on this record.*orbbrk*:First iteration in which orbitals are modified (default=3).*ibrstr*:Iteration increment between orbital updates (default=1).*ibrueck*:Scaling factor for singles in orbital updates (default=1).*brsfak*:

## Singles-doubles configuration interaction, CISD

Performs closed-shell configuration interaction, CISD. The same results as with the CI program are obtained, but this code is somewhat faster. Normally, no further input is needed. For specifying DIIS directives, see section the DIIS directive

## Quasi-variational coupled cluster, QVCCD

Performs closed-shell quasi-variational coupled cluster, QVCCD Normally the effect of single excitations needs to be included through orbital optimisation, and this can be done either through the Brueckner condition (`BQVCCD`

), or through variational minimisation of the energy functional with respect to orbital rotations (`OQVCCD`

).

J. B. Robinson and P. J. Knowles, J. Chem. Phys. **136**, 054114 (2012)

The effects of triple excitations can be included using the standard perturbation theory, `BQVCCD(T)`

or `OQVCCD(T)`

. The renormalised, `OQVCCDR(T)`

, and asymmetric renormalised, `OQVCCDAR(T)`

, methods can also be used to include the effects of the triples, while avoiding the overcorrection of the standard triples correction due to small energy differences appearing in the denominator.

J. B. Robinson and P. J. Knowles, Phys. Chem. Chem. Phys. **14**, 6729-6732 (2012)

J. B. Robinson and P. J. Knowles, J. Chem. Phys. **138**, 074104 (2013)

J. A. Black and P. J. Knowles, Mol. Phys. **116**, 1421-1427 (2018)

## Distinguishable cluster, DCD

`DCSD`

, `BDCD`

, or `ODCD`

perform closed-shell distinguishable cluster calculations with singles, Brueckner orbitals, or optimized orbitals, respectively. All CCSD and BCCD options (the closed shell CCSD program) can also be used for distinguishable cluster calculations. The analytical gradients (analytical energy gradients), explicit correlation (explicitly correlated methods), open-shell versions (open-shell coupled cluster theories), and a linear-scaling implementation (PAO-based local correlation treatments) are all available for the DCSD method. All reported calculations using this feature should cite:

D. Kats and F. R. Manby, J. Chem. Phys. **139**, 021102 (2013)

D. Kats, J. Chem. Phys. **141**, 061101 (2014)

D. Kats, Mol. Phys. **116**, 1435 (2018) (for Spin-Component Scaled DC)

## The DIIS directive

[sec:ccdiis]

`DIIS`

,*itedis,incdis,maxdis,itydis*;

This directive allows to modify the DIIS parameters for CCSD, QCISD, or BCCD calculations.

First iteration in which DIIS extrapolation may be performed (default=2).*itedis*:Increment between DIIS iterations (default=1).*incdis*:Maximum number of expansion vectors to be used (default=6).*maxdis*:DIIS extrapolation type. itydis=1 (default): residual is minimized. itydis=2: $\Delta T$ is minimized.*itydis*:

In addition, there is a threshold *THRDIS* which may be modified with the `THRESH`

directive. DIIS extrapolation is only done if the variance is smaller than *THRDIS*.

## Saving the wavefunction

`SAVE`

,*record*;

record name for save of wavefunction. The wavefunction is saved in each iteration. This has the same effect as the*record*:`SAVE`

option.

## Starting wavefunction

`START`

,*record*;

record name from which the wavefunction is restored for a restart. This has the same effect as the*record*:`START`

option. If the option`MAXIT=1`

is set on the CCSD command line, no further iterations will be performed. (Note: in this case there is no check whether the previous wavefunction was converged.) It is then for example possible to compute the triples correction using the restarted CCSD wavefunction.

## Examples

### Single-reference correlation treatments for H2O

- examples/h2o_ccsd.inp
***,h2o test geometry={o;h1,o,r;h2,o,r,h1,theta} !Z-matrix geometry input basis=vtz !cc-pVTZ basis set r=1 ang !bond length theta=104 !bond angle hf !do scf calculation text,examples for single-reference correlation treatments ci !CISD using MRCI code cepa(1) !cepa-1 using MRCI code mp2 !Second-order Moeller-Plesset mp3 !Second and third-order MP mp4 !Second, third, and fourth-order MP4(SDTQ) mp4;notripl !MP4(SDQ) cisd !CISD using special closed-shell code ccsd(t) !coupled-cluster CCSD(T) qci(t) !quadratic configuration interaction QCISD(T) bccd(t) !Brueckner CCD(T) calculation ---

### Single-reference correlation treatments for N2F2

- examples/n2f2_ccsd.inp
***,N2F2 CIS GEOMETRY (C2h) rnn=1.223,ang !define N-N distance rnf=1.398,ang !define N-F distance alpha=114.5; !define FNN angle geometry={N1 N2,N1,rnn F1,N1,rnf,N2,alpha F2,N2,rnf,N1,alpha,F1,180} basis=vtz !cc-pVTZ basis set $method=[hf,cisd,ccsd(t),qcisd(t),bccd(t)] !all methods to use do i=1,#method !loop over requested methods $method(i) !perform calculation for given methods e(i)=energy !save energy in variable e enddo !end loop over methods table,method,e !print a table with results title,Results for n2f2, basis=$basis !title of table

This calculation produces the following table:

Results for n2f2, basis=VTZ METHOD E E-ESCF CISD -308.4634948 -0.78283137 BCCD(T) -308.6251173 -0.94445391 CCSD(T) -308.6257931 -0.94512967 QCISD(T) -308.6274755 -0.94681207

## Expectation values

`EXPEC`

[,[`METHOD`

=]*method*][,[`TYPE`

=]*type*][,*opname*]

One-electron properties in closed-shell QCISD, QCISD(T), CCD, CCSD, CCSD(T), and DCSD can be computed by using `EXPEC`

. (the `GEXPEC`

directive has no effect in this case). The following options can be specified

`METHOD`

:*method*=`EDERIV`

: Default case. Properties are obtained as analytical energy derivatives.

*method*=`XCCSD`

: expectation-value CCSD method (with the exponential ansatz used also for bra) – see section First- and second-order properties for CCSD from expectation-value CC theory (XCCSD).

*method*=`XCCSD(3)`

: A simplified version of XCCSD – see section First- and second-order properties for CCSD from expectation-value CC theory (XCCSD). It is also available for the local CCSD method.

`TYPE`

:*type*=`RELAX`

: Compute fully relaxed properties with taking into account the orbital relaxation (default). Nonrelaxed properties are also computed and printed.

*type*=`NORELAX`

: Compute properties without taking into account the orbital relaxation. The `TYPE`

option has only effect for the analytical energy derivatives.

Note that, in CCSD, the orbital relaxation effect on properties is less pronounced than, for example, in MP2. By default, only dipole moments are computed. The syntax for operators *opname* is explained in section One-electron operators and expectation values (GEXPEC). See also `DM`

in section saving the density matrix, and `NATORB`

in section natural orbitals.

## Saving the density matrix

`DM`

,*record.ifil*;

The first-order density matrix in AO basis is stored in record *record* on file *ifil*. For details, see section saving the density matrix. See also `NATORB`

in section natural orbitals.

## Natural orbitals

`NATORB`

,[`RECORD=`

]*record.ifil*,[`PRINT=`

*nprint*],[`CORE`

[=*natcor*]];

Calculate natural orbitals. The options are the same as for MP2. For details, see section natural orbitals.

## Dual basis set calculations

Dual basis set calculations are possible with the closed-shell MP2 and CCSD codes (conventional and local, also with density fitting where available). Normally this means that the Hartree-Fock calculation is done with a smaller basis set than the correlation calculation. In MOLPRO, two possibilities exist: the recommended one is to perform the HF and MP2 or CCSD calculations using the same basis set, and only omit higher angular momentum functions in the HF. This means that the resulting HF orbitals can be used directly in the correlation calculation. Alternatively, one can use entirely different basis sets in HF and the correlation calculation; in this case the orbitals in the correlation calculation are determined by a least square fit to the HF orbitals. This is less efficient (in particular in fully direct calculations) and somewhat less accurate. In any case, a new Fock matrix is computed in the MP2/CCSD program and block diagonalized in the occupied and virtual orbital subspaces. A perturbative singles correction is applied in the MP2 in order to reduce the HF basis set error.

Typically, the input is as follows:

basis=vtz(d/p) !triple zeta basis set without f on heavy atoms and without d on hydrogen hf !Hartree-Fock in the small basis basis=vtz !full cc-pVTZ basis set to be used in ccsd ccsd(t),dual=1 !ccsd calculation

The option `dual=1`

is required, otherwise the program will stop with an error message saying that the basis set of the reference orbitals is inconsistent. This is a precaution in order to avoid unexpected results.

Similarly, this works for other closed-shell single reference methods such as `MP2, QCISD, MP2-F12, CCSD-F12`

, and for the local variants `LMP2, LQCISD, LCCSD`

in either conventional or direct mode. Furthermore, dual basis set `DF-LMP2, DF-LCCSD, DF-LMP2-F12`

calculations are possible.