Coupled cluster calculation

The ground state

To perform a coupled cluster calculation, you must specify a reference wavefunction along with the specific coupled cluster model in the method section. For instance, to do a standard CCSD calculation, you add the following method section:

method
  hf
  ccsd
end method

Furthermore, you must specify which type of calculation you wish to run. The options here are ground state, excited state, mean value, and response. Suppose we wish to simply find the ground state energy. Then we write:

do
  ground state
end do

This is all the information you have to specify to run with default settings. A minimal working example is the following:

system
  name: methanol
end system

method
  hf
  ccsd
end method

do
  ground state
end do

geometry
basis: cc-pVDZ
H        0.8475370000      0.0347360000      1.0345270000
C        0.3504820000      0.0067450000      0.0608350000
H        0.6350750000      0.8927640000     -0.5200650000
H        0.6629360000     -0.8933800000     -0.4828330000
O       -1.0108350000     -0.0082340000      0.3643870000
H       -1.4851950000     -0.0326310000     -0.4568520000
end geometry

Save this as methanol.inp and type the following command in your terminal (here shown for four threads)

eT_launch --omp 4

You should now have a file called methanol.out. If the calculation exited succesfully (look for eT terminated successfully! at the bottom of the file), the output file should show iteration information for three solvers: one to decompose the electron repulsion integrals, one to converge the Hartree-Fock equations, and one to converge the CCSD ground state equation. You should see something like this in the coupled cluster solver:

 ...
 8          -115.416194297528     0.1673E-04     0.3289E-05
 9          -115.416194207952     0.4979E-05     0.8958E-07
10          -115.416194185137     0.1589E-05     0.2281E-07
11          -115.416194095902     0.4640E-06     0.8924E-07
---------------------------------------------------------------
Convergence criterion met in 11 iterations!

When running calculations in eT, each specific equation has a corresponding solver section. In these sections, you can specify which algorithm to use, how tight the thresholds should be, and so on. For example, if you wish to use the newton-raphson solver instead of the default diis solver for the above calculation, you can add a solver cc gs section (solver CC ground state) to your methanol.inp:

solver cc gs
  algorithm: newton-raphson
end solver cc gs

When running this calculation, you may want to specify a different output-file to not overwrite the existing file. This is done by giving the - -output-file flag (or in short: -of) and a filename to the launch script

eT_launch --omp 4 --output-file methanol_newton_raphson.out

If successful, you’ll see something like this in the convergence information for the coupled cluster ground state:

   Micro-iter.  Residual norm
   --------------------------
     1           0.6875E-06
     2           0.2180E-06
     3           0.9140E-07
     4           0.2428E-07
     5           0.4667E-08
   --------------------------

Macro-iter.    Energy (a.u.)        |omega|       Delta E (a.u.)
----------------------------------------------------------------
    4         -115.416194137363     0.6568E-08     0.2980E-07
----------------------------------------------------------------
Convergence criterion met in 4 iterations!

By default, eT provides two main output files, one .out and one .timing.out. The latter contains detailed information regarding the time spent computing the most expensive terms. For example, in the methanol_newton_raphson.timing.out, you should find multiple instances of the time to construct the ground state equation (called “omega”):

Construct ccsd omega
wall time (sec):                 0.04
cpu time (sec):                  0.11
cpu/wall ratio:                  2.92

Note

The Newton-Raphson algorithm is typically more expensive than the standard algorithm. Micro- and macro-iterations are equally expensive. Below we switch back to the default algorithm.

Excitation energies

Suppose we want two excitation energies in addition to the ground state. Then we need to change the type of calculation. Enter the methanol.inp file and edit your do section to ask for excited states:

do
  excited state
end do

To perform an excited state calculation, you also need to specify the number of states. This information must be specified in the solver section for the excited state calculation. Add solver cc es section to your input file, and specify two singlet states:

solver cc es
  singlet states: 2
end solver cc es

Now run the calculation:

eT_launch --omp 4 -of methanol_2_es.out

If successful, you should find a summary of the excitation energies towards the end of the output file methanol_2_es.out, given in Hartree and eV:

   - Electronic excitation energies:

                                 Excitation energy
                     ------------------------------------------
 State                (Hartree)             (eV)
---------------------------------------------------------------
   1                  0.288688760560        7.855621303236
   2                  0.357786659905        9.735870915469
---------------------------------------------------------------
eV/Hartree (CODATA 2014): 27.21138602

As for the ground state, you can specify different thresholds. Say you only want to converge the excited states to \(10^{-4}\) in the residual norms. Then you can add two keywords to your methanol.inp file:

solver cc es
  singlet states:     2
  energy threshold:   1.0d-4
  residual threshold: 1.0d-4
end solver cc es

Rerunning the calculation, we again obtain excitation energies, but this time less tightly converged:

- Electronic excitation energies:

                                 Excitation energy
                     ------------------------------------------
 State                (Hartree)             (eV)
---------------------------------------------------------------
   1                  0.288689280999        7.855635465111
   2                  0.357786036365        9.735853948081
---------------------------------------------------------------
eV/Hartree (CODATA 2014): 27.21138602

It is also straight-forward to run core excitations instead of UV excitations. To do so, you can ask for the core-valence separation (CVS) approximation. Suppose we wish to look at excitations out of the orbital with the lowest energy (orbital number 1, i.e., the 1s orbital on the oxygen in methanol). Then we specify 1 using the core excitation keyword:

solver cc es
  singlet states:     2
  energy threshold:   1.0d-4
  residual threshold: 1.0d-4
  core excitation:    {1}
end solver cc es

Save and run:

eT_launch --omp 4 -of methanol_2_core_es.out

You should then see the following summary in your output file:

- Electronic excitation energies:

                                  Excitation energy
                      ------------------------------------------
  State                (Hartree)             (eV)
 ---------------------------------------------------------------
    1                 19.801727040165      538.832438352602
    2                 19.909586934852      541.767455583007
 ---------------------------------------------------------------
 eV/Hartree (CODATA 2014): 27.21138602

Oscillator strengths

If we want ground-to-excited-state intensities, we need to compute the transition moments. To do so, we first need to change the type of calculation. Enter methanol.inp and change the do section to the following:

do
  response
end do

In addition, we should specify that we want to use equation of motion (EOM), the dipole operator, and request transition moments (not polarizabilities). This is specified in the cc response section:

cc response
  eom
  dipole length
  transition moments
end cc response

Now run the calculation again:

eT_launch --omp 4 -of methanol_core_es_transitions.out

In the summary, you should now see the state energies along with transition moments (left and right), their components, and the oscillator strengths. For state 1, you should get something like the following:

State 1:
----------
Calculation type:                             EOM
Excitation energy [E_h]:          19.801727040165
Excitation energy [eV]:          538.832438352604
Hartree-to-eV (CODATA 2014):          27.21138602

              Transition moments [a.u.]         Transition strength [a.u.]
--------------------------------------------------------------------------
Comp. q     < k |q| 0 >       < 0 |q| k >        < 0 |q| k > < k |q| 0 >
--------------------------------------------------------------------------
X           0.0074932967      0.0160811966            0.0001205012
Y           0.0006087560      0.0012881757            0.0000007842
Z           0.0237667088      0.0501232307            0.0011912642
--------------------------------------------------------------------------
Oscillator strength [a.u.]:      0.017327165800

Polarizabilities

To get polarizabilities in addition to the transition moments, we can further specify this in the cc response section. Let us evaluate the \(xx\) and \(xz\) components of the polarizability at the frequencies \(0.1\) and \(0.2\). Then we can change the section to:

cc response
  eom
  dipole length
  transition moments
  polarizabilities: {11, 13}
  frequencies: {0.1d0, 0.2d0}
end cc response

Here, 11 and 13 refers to the \(xx\) and \(xz\) components.

At the very end of the output, you should see the evaluated polarizabilities:

6) Calculation of the EOM polarizabilities
   << mu_x, mu_x >>(0.10E+00):    -17.622255018069
   << mu_z, mu_x >>(0.10E+00):     -9.099903505128
   << mu_x, mu_x >>(0.20E+00):    -19.265671981087
   << mu_z, mu_x >>(0.20E+00):     -9.921342684364