The non-adiabatic coupling vectors between the states can be calculated with the SSR method. By setting NonAdiabaticCoupling = Yes in REKS block, the user can easily calculate the non-adiabatic coupling vectors. When this option is turned on, SSR shows the gradients for all states and possible coupling vectors:

--------------------------------------------------
--------------------------------------------------
1 st state (SSR)
0.09734680    -0.00000000     0.00000000
-0.09734682    -0.00000000    -0.00000000
0.00639517     0.00417744    -0.00000000
0.00639517    -0.00417744     0.00000000
-0.00639516     0.00417743    -0.00000000
-0.00639516    -0.00417743     0.00000000

2 st state (SSR)
-0.12803706    -0.00000000    -0.00000000
0.12803704    -0.00000000     0.00000000
0.00639486     0.00417837    -0.00000000
0.00639486    -0.00417837     0.00000000
-0.00639485     0.00417836    -0.00000000
-0.00639485    -0.00417836     0.00000000
--------------------------------------------------
Coupling Information
--------------------------------------------------
between  1 and  2 states
--------------------------------------------------
0.11269193    -0.00000000     0.00000000
-0.11269193     0.00000000    -0.00000000
0.00000016    -0.00000046     0.00000000
0.00000016     0.00000046     0.00000000
-0.00000016    -0.00000046     0.00000000
-0.00000016     0.00000046     0.00000000

h vector - derivative coupling
-0.00031874    -0.00000000     0.00000000
-0.00031874     0.00000000    -0.00000000
0.00015937     0.00007242     0.00000000
0.00015937    -0.00007242     0.00000000
0.00015937    -0.00007242     0.00000000
0.00015937     0.00007242     0.00000000

G vector - GDV
0.11269193    -0.00000000     0.00000000
-0.11269193     0.00000000    -0.00000000
0.00000016    -0.00000046     0.00000000
0.00000016     0.00000046     0.00000000
-0.00000016    -0.00000046     0.00000000
-0.00000016     0.00000046     0.00000000

H vector - DCV - non-adiabatic coupling
-0.00152718    -0.00000000     0.00000000
-0.00152718     0.00000000    -0.00000000
0.00076359     0.00034700     0.00000000
0.00076359    -0.00034700     0.00000000
0.00076359    -0.00034700     0.00000000
0.00076359     0.00034700     0.00000000
--------------------------------------------------


The above gradients and coupling vectors are obtained with the planar structure of a ethylene molecule. The g vector is defined as gradient difference vector, thus it can be calculated from the difference of SA-REKS gradients. Similarly to this, G vector is calculated from the difference of SSR gradients. The h vector is defined as coupling gradient, so it can be simply calculated from the gradient of state-interaction terms. The H vector is defined as the derivative of the coupling vectors, thus its norm increases as the energy gap becomes smaller.

The g and h vectors can be regarded as the vectors defined through diabatic states, and the G and H vectors are defined through the adiabatic (SSR) states. In general, the non-adiabatic coupling vectors can be used for surface hopping molecular dynamics, while the g and h vectors can be used for minimum energy conical intersection (MECI) optimisation.

## Relaxed Density¶

[Input: recipes/reks/relaxed_density/]

The relaxed density is calculated for the TargetState when the user sets RelaxedDensity to Yes. The calculation of a relaxed density requires the information about gradient, thus it can be calculated when the input enables gradient calculation. When this option is turned on, the relaxed FONs are given in the bottom of standard output and the relaxed_charge.dat file is generated. It includes the total charge as well as the Mulliken charges of each atom for target state. This example shows relaxed charges for the ground state:

total charge:     -0.00000000 (e)

relaxed S0 atomic charge (e)
atom        charge
1      -0.18168616
2      -0.18168616
3       0.09084308
4       0.09084308
5       0.09084308
6       0.09084308


If one want to exploit SSR with a QM/MM approach, the RelaxedDensity option should be used to obtain the gradient of external point charges. Then, the gradient of external point charges for target state are given in detailed.out file. Therefore, this option can be used to run surface hopping dynamics with a QM/MM approach.

## Spin tuning constants¶

The DFTB/SSR method well describes equilibrium geometries and vertical excitation energies as compared with SSR/wPBEh results. (See the paper JCTC, 2019, 15, 3021-3032.) However, the behaviours at MECI points sometimes does not match those obtained with SSR/wPBEh. For example, the n/$$\pi^*$$ type MECI geometry of ethylene or methyliminium molecule cannot be located with DFTB/SSR, with an incorrect description of the relative stability of the PPS and OSS states being mostly responsible. Their relative stability depends on the splitting between the open-shell singlet microstates and the triplet microstates in the PPS and OSS energies.

In principle, DFTB/SSR employs spin-polarised DFTB formalism, in which the spin-polarisation contribution is obtained from the second-order expansion of the magnetisation density with respect to zero magnetisation electronic structure. At the n/$$\pi^*$$ type MECI geometries, both frontier orbitals are located on the same atom. In such a case, the second-order expansion of magnetisation may not be suitable for the triplet microstates, as the spin density becomes too large. As a simple solution, the stability between the PPS and OSS states can be adjusted by scaling the atomic spin constants. For most molecules the FONs for the PPS state become $$n_a$$ ~ 2.0 and $$n_b$$ ~ 0.0, hence the energy of the PPS state is determined by the 1st microstate alone and it is only the energy of the OSS state that depends on the atomic spin constants.

If the user runs the test calculation included with the main DFTB+ repository in the test/prog/dftb+/reks/PSB3_2SSR_rangesep_tuning directory, the following results are given in the standard output:

----------------------------------------------------------------
SSR: 2SI-2SA-REKS(2,2) states
E_n       C_{PPS}    C_{OSS}
SSR state  1  -16.39950035  -0.955182  -0.296020
SSR state  2  -16.38979921   0.296020  -0.955182
----------------------------------------------------------------

H vector - DCV - non-adiabatic coupling
-0.11443527     0.16066530     0.18737790
0.15477051    -0.03677914     0.02047645
-0.81366262     0.56406749    -1.41489688
0.69877935    -2.68846832     1.70292281
-0.34893713     0.81739627     0.50771020
-0.08309257     0.13356401     0.02704612
0.99884931     0.52013559    -1.04847579
-0.46122482     0.81158082    -0.36230963
-0.03763611     0.00213152     0.02451782
0.65298995     0.42710364    -0.20222417
-0.52172598    -0.75236687     1.00033386
-0.24407958    -0.11018546    -0.50337767
0.04971872    -0.03074284     0.01547016
0.06968625     0.18189801     0.04542884


It shows the energies at MECI point of a PSB3 molecule, thus the non-adiabatic coupling vectors show large elements near to the centre C=C bond. The atomic spin constants can be modified by using SpinTuning keyword in REKS block as follows:

Reks = SSR22 {
Energy = {
Functional = { "PPS" "OSS" }
StateInteractions = Yes
}
TargetState = 2
FonMaxIter = 30
shift = 0.3
SpinTuning = { 3.2 3.2 3.2 }
TransitionDipole = Yes
CGmaxIter = 100
Tolerance = 1.0E-8
Preconditioner = Yes
SaveMemory = Yes
}
RelaxedDensity = Yes
VerbosityLevel = 1
}


## Microstate calculation¶

[Input: recipes/reks/microstate/]

Obviously, the SSR method treats only singlet states like PPS or OSS. If one want to compare the energy of singlet and triplet states, SSR provides the energy of a triplet configuration as an alternative which corresponds to the 5th or 6th configuration in the (2,2) active space. Thus, the user can easily compare the energy of singlet and triplet microstates.:

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


In this example, for a distorted structure of ethylene, the energy of triplet microstate is close to that of the PPS state, since the frontier orbitals are the localised $$\pi$$ orbitals in this system. If one want to know the gradient of the triplet microstate as well as its relaxed density, these can be obtained by using TargetMicrostate keyword in REKS block. If the value for this keyword is to 5, then the properties will be calculated according to the index of the microstate. In a (2,2) active space, the 5th microstate indicates a triplet configuration, thus the output shows the quantities for this microstate. The following results are obtained from the distorted structure of ethylene molecule:

--------------------------------------------------
--------------------------------------------------
5 microstate
-0.00926010    -0.00000000     0.00000000
0.00926011     0.00000000    -0.00000000
0.00096561     0.00066166     0.00000008
0.00096561    -0.00066166    -0.00000008
-0.00096562    -0.00000009     0.00066176
-0.00096562     0.00000009    -0.00066176
--------------------------------------------------


The gradient is now calculated for the 5th microstate. In addition, the energy of spin contribution in the detailed.out file is -0.023, Hartree which corresponds to the spin constant $$W_{pp}$$ for a carbon atom. In this case the frontier orbitals consisted of only p orbitals of a carbon atom, thus the energy of the spin contribution mostly consists of interactions between these orbitals. With this option, one can run the molecular dynamics simulation for the triplet microstate.

This example shows an input file for calculation of a triplet microstate:

Reks = SSR22 {
Energy = {
Functional = { "PPS" "OSS" }
}
TargetState = 1
TargetMicrostate = 5
FonMaxIter = 100
shift = 20.0

Note that TargetMicrostate keyword can be used only with the SA-REKS input settings discussed above.