# Energy and gradient calculations with REKS¶

## Single-state REKS¶

[Input: recipes/reks/single_state_reks/]

In single-state REKS, the ground state is calculated including multi-reference character. For example, a ethylene molecule with a torsional angle of 90 degrees has degenerate \({\pi}\) and \({\pi}^*\) orbitals, corresponding to the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) respectively. In this case a standard SCC-DFTB calculation can show oscillating phenomena when trying to obtain self-consistent charges and is missing the static electronic correlations required for a correct variational energy. Similarly, REKS can deal with bond-breaking, di-radicals and magnetic couplings which are difficult in standard calculations.

To deal with the static electronic correlation, REKS first determines the size of the space of active orbitals. Only two electrons in two frontier orbitals, (2,2), is supported in current version of DFTB+.

Note

From the active space, REKS automatically construct possible
electronic configurations, called microstate (See the paper JCTC, 2019, 15,
3021-3032.). A total
of 6 microstates are generated in a REKS(2,2) calculations, for example the
5^{th} and 6^{th} microstates express high-spin configurations in
the (2,2) space. To correctly describe the spin polarised configurations in
the space the `SpinConstants`

block is required in the `Hamiltonian`

. For
example, there are more up-electrons in the 5^{th} microstate than
down-electrons, thus atomic spin constants are needed to describe spin
contributions. Note that the user should not include the `SpinPolarisation`

block, only the `SpinConstants`

block is needed.

With this active space, single-state REKS can be carried out using the following input:

```
Hamiltonian = DFTB {
SCC = Yes
SpinConstants = {
ShellResolvedSpin = Yes
H = { -0.072 }
C = { -0.031 -0.025 -0.025 -0.023 }
}
...
}
REKS = SSR22 {
Energy = {
Functional = { "PPS" }
}
TargetState = 1
FonMaxIter = 100
Shift = 1.0
VerbosityLevel = 1
}
```

Note that the `REKS = SSR22`

block is a separate environment in the
calculation, so is not included in the `Hamiltonian`

block. In a REKS
calculation, the `Hamiltonian`

only controls detailed Hamiltonian settings
such as use of a range-separated hybrid functional or dispersion
corrections. Three singlet states can be generated from a (2,2) active
space. These are the perfectly spin-paired singlet state (PPS), an open-shell
singlet state (OSS) and the doubly excited singlet state (DES). Combining these
singlet states, we can find the many-particle singlet state which minimises the
energy functional. In the single-state REKS case we treat only one singlet state
of PPS symmetry, thus this is the only possibility that can be selected in the
`Functional`

block of the `Energy`

in this case.

The energy for the PPS state is expressed as

where \(C_L^{PPS}\) are the weighting factors of each microstate in the resulting PPS state and \(E_L\) are their energies. Here the weighting factors are calculated using the fractional occupation numbers (FONs) of frontier orbitals, thus the energy of the resulting PPS state is minimised with respect to both the Kohn-Sham orbitals and the FONs in single-state REKS.

[Input: recipes/reks/single_state_reks/dftb_in.hsd-0deg]

With this inputs, the standard output shows information about the energy of the PPS state and the FONs contributing to it. The minimised geometry for ethylene is a planar structure and the resulting energy and FONs are given in the standard output:

```
--------------------------------------------------
Final REKS(2,2) energy: -4.89357215
State Energy FON(1) FON(2) Spin
PPS -4.89357215 1.999990 0.000010 0.00
--------------------------------------------------
```

The energy of PPS state is -4.8936 Hartree and the FONs of frontier orbitals are ~2.0 and ~0.0, respectively. The orbital character for HOMO is \(\pi\) and the character for LUMO is \(\pi^*\), so the two electrons are occupied in the HOMO.

[Input: recipes/reks/single_state_reks/dftb_in.hsd-90deg]

When the torsional angle begins to rotate to 90 degree, the molecular orbitals change from \(\pi\) and \(\pi^*\) orbitals into two localised single \(\pi\) orbitals. This means the frontier orbitals are now almost degenerate, and this is well described with REKS(2,2) calculations. The resulting energy and FONs are:

```
--------------------------------------------------
Final REKS(2,2) energy: -4.77774313
State Energy FON(1) FON(2) Spin
PPS -4.77774313 1.000004 0.999996 0.00
--------------------------------------------------
```

Now, the energy of PPS state becomes -4.7777 Hartree and the FONs become ~1.0
and ~1.0, respectively. In addition, the user can check the energy contributions
to the PPS state from *detailed.out* file. The energy from the spin contribution
is ~0.0 Hartree for the planar structure, while the energy becomes -0.0188
Hartree for the 90 degree rotated structure. Overall, we can conclude that the
planar structure almost consists of only the first microstate, while the rotated
structure consists of several microstates.

## SA-REKS¶

[Input: recipes/reks/sa_reks/]

Single-state REKS can treat only the ground state, thus the vertical excitation energy cannot be calculated with this method. From the restricted open-shell Kohn-Sham scheme, we can construct the OSS state, which is expressed as

where \(C_L^{OSS}\) is weighting factors of each microstate making up the
OSS state. By minimising the energy for the ensemble of PPS and OSS states, we
can calculate the vertical excitation energy between them with state-average
REKS (SA-REKS). Again we calculate the energy of the PPS and OSS states of the
two geometries of an ethylene molecule with SA-REKS. The `REKS`

block has now
an additional `Gradient`

block to calculate the gradient and optimise for the
target state:

```
REKS = SSR22 {
Energy = {
Functional = { "PPS" "OSS" }
}
TargetState = 1
FonMaxIter = 100
Shift = 1.0
Gradient = ConjugateGradient {
CGmaxIter = 100
Tolerance = 1.0E-8
Preconditioner = Yes
SaveMemory = Yes
}
VerbosityLevel = 1
}
```

The user also now has to include the OSS state in addition to the PPS in the
`Functional`

block so that the energy of both are now calculated. The
resulting energy and additional information is given by:

```
--------------------------------------------------
Final SA-REKS(2,2) energy: -4.78921495
State Energy FON(1) FON(2) Spin
PPS -4.89357215 1.999990 0.000010 0.00
OSS -4.68485776 1.000000 1.000000 0.00
Trip -4.73085776 1.000000 1.000000 1.00
--------------------------------------------------
Lagrangian Wab: 0.00000000 0.00000000
--------------------------------------------------
SSR: 2SI-2SA-REKS(2,2) Hamiltonian matrix
PPS OSS
PPS -4.89357215 0.00000000
OSS 0.00000000 -4.68485776
--------------------------------------------------
unrelaxed SA-REKS FONs for S0: 1.999990 0.000010
```

The final SA-REKS(2,2) energy is from the ensemble of the PPS and OSS states,
the energy of which is being minimised. For the planar structure of an ethylene
molecule, the energies of two states are -4.8936 and -4.6849 Hartree,
respectively. The FONs for the PPS state are ~2.0 and ~0.0, while those for OSS
state are ~1.0 and ~1.0. In addition, the energy of the triplet configuration
which corresponds to contributions from the 5^{th} and 6^{th}
microstates is now given in the standard output. Note that this energy is *not*
the energy of the triplet state. The user can check for successful convergence
by comparing the two Lagrangian *Wab* values, these will become almost the same
if the energy is well converged.

After the energy calculation is finished, the gradient for the target state
(`TargetState = 1`

, which is the PPS state in this example) is calculated and the
final gradient appears at the bottom of the standard output. The keywords in the
`Gradient`

block affect coupled-perturbed REKS (CP-REKS) equations which are
used to calculate the gradient of target state. Here we choose a conjugate
gradient solver and the keywords `Preconditioner`

and `SaveMemory`

are used
to accelerate the computational speed of CP-REKS. These two keywords may be
switched off depending on the user’s system.

Similar to the case above, the distorted structure can be calculated using SA-REKS and the results are given in the following:

```
--------------------------------------------------
Final SA-REKS(2,2) energy: -4.75910919
State Energy FON(1) FON(2) Spin
PPS -4.77753765 1.000000 1.000000 0.00
OSS -4.74068073 1.000000 1.000000 0.00
Trip -4.77753837 1.000000 1.000000 1.00
--------------------------------------------------
Lagrangian Wab: -0.00004046 0.00004230
--------------------------------------------------
SSR: 2SI-2SA-REKS(2,2) Hamiltonian matrix
PPS OSS
PPS -4.77753765 -0.00000000
OSS -0.00000000 -4.74068073
--------------------------------------------------
unrelaxed SSR FONs for S0: 1.000000 1.000000
```

Now the energy of OSS state is -4.7407 Hartree. The energy of the triplet
configuration is similar to the energy of PPS state. Since REKS can consider the
static electronic correlations, it can produce the correct shape for the
potential energy curve with respect to the torsional angle of C=C bond. If you
want to calculate the energy of the DES state, then `IncludeAllStates = Yes`

keyword in the `Energy`

block will produce the energy of DES state as well as
its FONs.

## SI-SA-REKS¶

[Input: recipes/reks/si_sa_reks/]

State-interaction SA-REKS (SI-SA-REKS, briefly SSR) energies are obtained by solving a 2 \(\times\) 2 secular equation with the possible couplings between the SA-REKS states.

By considering the state-interaction terms, the SSR states become more reliable
when the excited states are included. The SSR states can be calculated with the
`StateInteractions = Yes`

in `Energy`

block. For the planar structure, the
resulting energies are given by:

```
----------------------------------------------------------------
SSR: 2SI-2SA-REKS(2,2) states
E_n C_{PPS} C_{OSS}
SSR state 1 -4.89357215 -1.000000 0.000000
SSR state 2 -4.68485776 -0.000000 -1.000000
----------------------------------------------------------------
```

In this case, the ground state consists of the PPS state, while the lowest excited state is of OSS type. As the coupling term increases, the difference between the SA-REKS and SSR energies becomes larger.