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:

• NOCHECK Ignore convergence checks.
• DIRECT Do calculation integral direct.
• NOSING Do not include singly external configurations.
• MAXIT=value Maximum number of iterations.
• SHIFTS=value Denominator shift for update of singles.
• SHIFTP=value Denominator shift for update of doubles.
• THRDEN=value Convergence threshold for the energy.
• THRVAR=value Convergence threshold for CC amplitudes. This applies to the square sum of the changes of the amplitudes.
• SAVE=record Save CCSD wavefunction to this record
• START=record Restart CCSD wavefunction from a previously written record
• HO=operator Specifies the zeroth-order Hamiltonian in closed-shell CCSD calculations with KS orbitals. 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.
• KSFOCK Forces the use of the KS Fock operator in the subsequent calculation. This is required for double-hybrid and RPA calculations.

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.

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). If KS or other non-HF reference states (like the Brueckner determinant) are employed in the calculation of the triples (T) additional singles terms are computed, see J. D. Watts, J. Gauss, R. Bartlett J. Chem. Phys. 98, 8718 (1993). This can be switched on/off with the option TRIP_FIA=1/0. If TRIP_FIA=1 the second order singles energy is computed in the triples program and if greater in magnitude than THR_E2SING the additional terms will be computed. THR_E2SING=1d-6 is the default setting for this threshold, which can also be altered via giving this keyword as an option to the CCSD program.

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).

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:

• TYPE=REF Compute and store density of reference determinant only (default). This corresponds to the BOX (Brueckner orbital expectation value) method of Chem. Phys. Lett. 315, 248 (1999).
• TYPE=TOT Compute and store density with contribution of pair amplitudes (linear terms). Normally, this does not seem to lead to an improvement.
• TYPE=ALL Compute and store both densities

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).

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.

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

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

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)

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)

[sec:ccdiis]

DIIS,itedis,incdis,maxdis,itydis;

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

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

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.

SAVE,record;

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

START,record;

• record: record name from which the wavefunction is restored for a restart. This has the same effect as the 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/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
---

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

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.

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.

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 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.