Introduction

    Thomas Relativistic Electronic Structure Calculation (TRESC) is developed to calculate the electronic structure of non-periodic polyatomic systems under the Born-Oppenheimer approximation, it supports Hartree-Fock/Kohn-Sham self-consistent field (HF/KS-SCF) calculation based on 2-component DKH2 Hamiltonian (2c-Hamiltonian) of a given molecule.<br>     The main part of the code is written in Fortran 2008 free format.

Algorithms

• HF/KS-SCF calculation is based on spherical-harmonic fragment-contracted Gaussian type orbital.
• Initial guess load from MOLDEN format file (generated by other programs, same or different basis set) or .ao2mo binary file (generated by TRESC HF/KS-SCF earlier, same basis set).
• Symmetric orthogonalisation of overlap matrix by default, canonical orthogonalisation will be used if basis linear dependence reaches the threshold.
• Relativistic spinor integrals based on optimal-parametrized 2nd order Douglas- Kroll-Hess(DKH2) transformation proposed by Hess et al., including DKH2 1e-integrals and DKH2 2e-integrals.
• DKH2 1e-integrals are calculated to the order $$c^{-4}$$, which contains scalar realativistic terms and one-electron spin-orbit coupling terms; Gaussian finite nucleus model is considered, affecting the integral value of $$\left<V\right>,\left<pVp\right>, \left<pppVp\right>$$.
• DKH2 2e-integrals are calculated to the order $$c^{-2}$$, which contains realativistic Coulomb, Exchange terms and so-called spin-same-orbit coupling terms; all four-center integrals utilize permutation symmetry and Cauchy-Schwarz screening technique.
• Construct Fock matrix via direct way, which is time consuming but less demanding on memory and disk r&w.
• Grid integration are based on the Becke's fuzzy partitioning, the exchange-correlation energy and the partial derivative terms of the exchange-correlation potential are obtained by external library Libxc.
• Support for density functionals: (hybrid) LDAs, (hybrid) GGAs and (hybrid) meta-GGAs excluding long-range corrected functionals and non-local correlation, the results of non-relativistic calculation differ negligibly from the Psi4 program.
• 2c-Hamiltonian causes mixing between alpha and beta orbitals, sometimes we want this mixing as little as possible (maintaining the initial spin state as much as possible); in this case, should try cspin=f or cspin=d (the latter is for nearly degenerate frontier orbitals), but these methods may cause variational instability, so be sure to check the convergence of the last few iterations.
• DIIS (Pulay mixing) can be used to accelerate HF/KS-SCF, dynamic damping can be used to enhance convergence.
• Basic linear algebra is computed using LAPACK subroutines.
• 1e and 2e Fock matices construction and grid-based integration support OpenMP parallel computation and all parallel zone are thread safe.
• Outprint $$\left< s^2 \right> \left( L\ddot{o}wdin \right)$$, energy and orbital components.
• Support dispersion correction via DFT-D4 program (stand-alone) developed by Grimme's group.

For mathematical and algorithmic details, see docs/Mathematics_and_Algorithms.pdf.

Characteristic

A special Hamiltonian: SRTP

    Second Relativized Thomas Precession (SRTP) is to conbine the Lorentz vector feature of the spin 4-vector $$\left(0,\vec{s}\right)$$ and the Lorentz scalar feature of the magnitude of its spatial components ($$s=\hbar/2$$). 'Second Relativized' means the magnitude of spin vector is independent of the reference frame choice.<br>     To accomplish it, we start with a newly-defined reference frame transformation rule, which makes the observed $$\vec{s}/s$$ from any frame identical with the observed $$\vec{s}/s$$ from corresponding frame under the Lorentz transformation rule, but magnitude $$s$$ always $$\hbar /2$$.<br>     Assuming that frame O' is moving along the x-axis in frame O, the Lorentz transformation and the newly-defined transformation lead to different observation.<br>

<p align="center"> <img src="docs/Observation_of_2_transformations.png" alt="docs/Observation_of_2_transformations.png" width="500"> <br> <em>observation of Lorentz transformation and newly-defined transformation</em> </p>

    Its mathematical form can be given directly as a nonlinear equation<br>

$$ \left( \begin{array}{c} 0\ s_1\prime\ s_2\prime\ s_3\prime\ \end{array} \right) =\left( \begin{matrix} 1& 0& 0& 0\ -\gamma \beta _1\zeta& \left[ 1+\frac{\left( \gamma -1 \right) \beta _{1}^{2}} {\beta ^2} \right] \zeta& \frac{\left( \gamma -1 \right) \beta _1\beta _2} {\beta ^2}\zeta& \frac{\left( \gamma -1 \right) \beta _1\beta _3}{\beta ^2} \zeta\ -\gamma \beta _2\zeta& \frac{\left( \gamma -1 \right) \beta _1\beta _2} {\beta ^2}\zeta& \left[ 1+\frac{\left( \gamma -1 \right) \beta _{2}^{2}} {\beta ^2} \right] \zeta& \frac{\left( \gamma -1 \right) \beta _2\beta _3} {\beta ^2}\zeta\ -\gamma \beta _3\zeta& \frac{\left( \gamma -1 \right) \beta _1\beta _3} {\beta ^2}\zeta& \frac{\left( \gamma -1 \right) \beta _2\beta _3}{\beta ^2} \zeta& \left[ 1+\frac{\left( \gamma -1 \right) \beta _{3}^{2}}{\beta ^2} \right] \zeta\ \end{matrix} \right) \left( \begin{array}{c} 0\ s_1\ s_2\ s_3\ \end{array} \right) $$

    where $$s_i$$ represent spin components and $$\beta _i$$ represent velocity components, $$\gamma$$ represent Lorentz factor and<br>

$$ \zeta =\left( 1+\gamma ^2\left( \vec{\beta}\cdot \frac{\vec{s}}{s} \right) ^2 \right) ^{-1/2} $$

    This newly-defined transformation is kinematic, but it will change the form of Thomas precession dynamically since Thomas precession is related to the intrinsic property of reference frame transformation.<br>     After some derivation, the contribution of the Thomas precession to electron energy at low speed can be represented as<br>

$$ H_{\mathrm{SRTP}}=\frac{1}{2}\vec{s}_{\gamma}\cdot \left( \dot{\vec{\beta}}\times \vec{\beta} \right) $$

    where<br>

$$ s_{\gamma ,i}=\frac{1}{\sqrt{1-\beta _i^{2}}}s_i $$

    Then quantization and use Pauli vector rule to modify Dirac matrices as<br>

$$ \alpha _i=\left( \frac{\left( 1-\beta _{j}^{2} \right)\left( 1-\beta _{k}^{2} \right)}{1-\beta _{i}^{2}} \right) ^{\frac{1}{4}}\left( \begin{matrix} & \sigma _i\ \sigma _i& \ \end{matrix} \right) $$

    This formular leads to the modified electron spinor wavefunction through DKH transformation.<br>     In addition, SRTP effect is of order $$c^{-4}$$, one have to consider other terms of order $$\geqslant c^{-4}$$ before it, including radiation effect. Moreover, the lowest order of SRTP still requires the computation of integrals like $$\langle i|p_{x}^{3}V_{ij}p_y|j\rangle$$, it has a small effect on results but will significantly increases the one-electron integral cost.<br>     SRTP is currently has no evidence support, if you're interested, try keyword pppVp in %Hamiltonian when performing DKH2 calculation.<br>

Characterization of spinor states

Decomposition to spin-pure states

    A given spinor state $$|\psi \rangle$$ can always be decomposed into series of spin-pure states $$|S,M\rangle$$ by rotation group integration(RGI):

$$ \langle \varPsi |S,M\rangle =\frac{2S+1}{8\pi ^2} \int{\mathrm{d}\varOmega}D_{MM}^{S*}\left( \varOmega \right) \langle \varPsi | R\left( \varOmega \right) |\varPsi \rangle $$

    where $$D_{MM}^{S}$$ is the Wigner D-matrix and $$R$$ is the SU(2) rotation matrix. TRESC will calculates the contribution of several spin-pure states to the converged spinor state, any possible spin-pure state is predicted by the Wigner-Eckart theorem.

Measure of Time-Reversal Symmetry(TRS) breaking

    TRESC calculates the $$\kappa$$ parameter by:

$$ \begin{aligned} \kappa&=\left| MM^*+I_N \right|\ M_{ij}&=\langle \phi _i|-i\sigma_y|\phi _j\rangle \end{aligned} $$

    where $$N$$ is the number of electrons, $$|\phi_i\rangle$$ denotes the i-th occupied orbital.<br>     The deviation of $$\kappa$$ from 0 reflects the degree of TRS breaking of the system. For scalar single-configuration wavefunction, the reference $$\kappa$$ is $$\sqrt{N_{\alpha}-N_{\beta}}$$; deviation from this value after 2c-Hamiltonian HF/KS-SCF calculation reflects the degree of TRS breaking caused by SOC, which is closely related to the zero-field splitting of the system.

Visualisation of complex spinor MOs

    Vis2C is a Python package for visualising spinor (and scalar) molecular orbitals. It supports three visualization methods, which differ in terms of projection space and grid partitioning. All grid parameters of these methods can be set in $TRESC/gridsettings.ini.<br>

1. Structured grid data: cub1c/cub2c

    When TRESC finishes its HF/KS-SCF calculation, canonical orbitals will be dumped to jobname.molden.d. With it, one can use cub2c function to generate two GAUSSIAN cube format files (contain grid data of real and imaginary part of alpha and beta components of the selected orbital) and then visualise the selected orbital based on the grid data automatically. The visualisation are as follows:<br>

<p align="center"> <img src="docs/Ni-C2H4-3-AlphaHOMO-1.png" alt="docs/Ni-C2H4-3-AlphaHOMO-1.png" width="500"> <br> <em>Ni-C2H4-3-AlphaHOMO-1</em> </p>

    The visualisation shows spin phase and orbital phase, variance in spin phase reflects the strength of spin-orbit coupling (SOC) of the selected orbital. Some frontier spinor orbitals in examples/ shows the fact that orbital with more nodal surfaces in regions near heavy atoms exhibit stronger SOC and are mor