Introduction

    The Thomas Relativistic Electronic Structure Calculation (TRESC) package is designed to compute the electronic structure of non-periodic polyatomic systems within the Born-Oppenheimer approximation. It supports Hartree-Fock and Kohn-Sham self-consistent field (HF/KS-SCF) calculations based on a two-component (2c) second-order Douglas-Kroll-Hess (DKH2) Hamiltonian for a given structure.
    TRESC exhibits strong performance across both main-group molecules and transition metal complexes, see /examples. The core computational engine (called tkernel) is implemented 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).
• All potential integrals are computed using the Obara-Saika scheme (2e-integrals can be switched to the Rys Quadrature scheme), and all 2e-integrals are handled using PRISM.
• 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. When the RMSDP is sufficiently small, the effect of density matrix will be considered in screening as: $$P_{\max}\sqrt{\left( \mu \nu |\mu \nu \right) \left(\sigma \lambda |\sigma \lambda \right)} \leq threshold$$. Similar algorithm has been employed in the CP2K program.
• 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), dynamic damping, virtual orbital level shifting can be used to facilitate HF/KS-SCF process.
• 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) combines the Lorentz vector feature of the spin 4-vector $$\left(0,\vec{s}\right)$$ with the Lorentz scalar feature of the magnitude of its spatial components ($$s=\hbar/2$$). The term 'Second Relativized' signifies that the magnitude of spin vector remains strictly independent of the chosen reference frame. To achieve this, we introduce a newly-defined frame transformation rule, which ensures that the observed $$\vec{s}/s$$ in any frame is identical to that derived under standard Lorentz transformations while its magnitude $$s$$ is invariantly preserved as $$\hbar /2$$.
    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.

observation of Lorentz transformation and newly-defined transformation

    Its mathematical form can be given directly as a nonlinear equation

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

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

    Although this newly proposed transformation is essentially kinematic, it dynamically alters the form of Thomas precession, as this precession is intrinsically linked to the properties of reference frame transformations.
    Following detailed derivation, the contribution of Thomas precession to the electronic energy in the low-velocity limit can be expressed as

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

    where

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

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

$$ \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 formula yields the modified electron spinor wavefunction through the DKH transformation. Furthermore, since the SRTP correction scales as $$\mathcal{O}(c^{-4})$$, one must first account for other terms up to this order, including radiation effects. Moreover, evaluating the lowest-order SRTP correction requires computing complex integrals such as $$\langle i|p_{x}^{3}V_{ij}p_y|j\rangle$$, while this yields only a minor impact on the final results, it significantly increases the computational cost of one-electron integrals.
    Currently, SRTP lacks empirical support; however, interested users can explore this effect by adding the pppVp keyword to the %Hamiltonian block when performing DKH2 calculations.

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 calculates the contributions of several spin-pure states to the converged spinor state, with all permissible spin-pure states 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. The deviation of $$\kappa$$ from 0 reflects the degree of TRS breaking within the system. For a scalar single-configuration wavefunction, the reference $$\kappa$$ value is $$\sqrt{N_{\alpha}-N_{\beta}}$$; deviation from this reference following a 2c-Hamiltonian HF/KS-SCF calculation indicates the extent of SOC-induced TRS breaking which is closely related to the zero-field splitting of the system.

Visualisation of complex spinor MOs

    Vis2c is a Python package designed to visualize scalar and spinor molecular orbitals, offering three visualization methods based on different projection spaces and grid partitioning strategies. Grid parameters of all these methods can be changed in $TRESC/gridsettings.ini.

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, users can execute the cub2c.py 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:

HOMO-3 of triplet [Ni-C2H4], <s^2>=0.67 , <s_z>=-0.41

    The visualization displays both the spin and orbital phases; the variance in the spin phase reflects the strength of the spin-orbit coupling (SOC) within the selected orbital. As demonstrated by the frontier spinor orbitals in the examples/, orbitals with more nodal surfaces near heavy atoms exhibit stronger SOC and are more susceptible to spin flipping.

    Furthermore, executing cub2c.py with the -slice flag generates cross-sectional slices of both the orbital amplitude and the spin phase. These visualizations reveal that the amplitude slices are typically solid, whereas the spin-phase slices appear hollow, indicating that spin perturbations predominantly occur at the orbital phase boundaries.

spin-phase slice of HOMO-3 of triplet [Ni-C2H4], a distinctive hollow, umbrella-shaped pattern

2. Untructured grid data: mog2c

    While structured grid data facilitates efficient post-processing and visualization, unstructured grids become necessary in the following three scenarios:

1. Visualizing a specific orbital localized within the inner shell of a particular atom in a large molecule;
2. Illustrating detailed phase changes induced by SOC in the vicinity of heavy nuclei;
3. Capturing the frame-dragging effect of a rapidly moving molecule;

    Vis2c supports visualization using unstructured data, specifically Becke's fuzzy grids. By utilizing the jobname.molden.d generated by TRESC, users can execute the mog2c.py to automatically generate the requisite Becke grid data and visualize the selected orbital.

3. Replace phase with momentum: pro1c/pro2c

    Momentum space is the dual of real space, where spatial phase information is transformed into momentum amplitude. Consequently, pro2c.py avoids complex phase plotting and instead utilizes intuitive isosurface representations to convey all orbital information. The visualizations are as follows:

HOMO-3 of triplet [Ni-C2H4], d-orbital in momentum space also has a quatrefoil shape.

    The momentum-space projection of a bound-state orbital is inherently centrally symmetric, exhibiting symmetry about the origin (indicated by the central black dot). Furthermore, mapping to momentum space emphasizes the particle-like nature of electrons; for instance, the conventional spatial phase-matching principle for bonding orbitals is effectively replaced by a momentum-matching principle.

Options

    TRESC allows to adjust computation by providing keywords in each module, now list the keywords currently supported.

%ATOMS (modeling)

• basis: choice of basis set
• geom: XYZ geometry file in working directory
• charge: charge of molecule
• spin: spin multiplicity of molecule

%HAMILTONIAN (electronic integrations)

• pVp1e: one-electron pVp integrals (DKH2 spinor)
• pVp2e: two-electron pVp integrals (DKH2 spinor)
• pVp: equal to pVp1e + pVp2e
• pppVp: one-electron pppVp integrals (SRTP effect)
• cutS: threshold of eigenvalue of overlapping matrix for orthogonalisation
• threads: number of threads in multi-threaded calculation

%SCF (SCF settings)

• guess: initial guess of wavefunction
• schwarz: threshold of Cauchy-Schwarz screening
• dmschwarz: RMSDP threshold for considering density matrix in Cauchy-Schwarz screening
• maxiter: maximum number of SCF iterations
• convertol: SCF convergence tolerance
• damp: SCF damping coefficient
• diisdamp: DIIS damping coefficient
• nodiis: number of initial iterations without DIIS
• subsp: size of subspace of DIIS
• cutdamp: cut damp when threshold is reached
• cutdiis: cut DIIS when threshold is reached
• prtlev: minimum AO coefficient output when SCF done
• cspin: constrained spin multiplicity calculation
• molden: dump MOs to MOLDEN files when SCF done
• emd4: DFT-D4 dispersion correction
• lshift: virtual orbital energy level shift (Eh), better be in range [0.1,1]

%FUNCTIONAL (exchange-correlation functionals)

Identity of functionals can be found at Libxc
All functionals are available except range-separated and non-local correlation

• xid: identity of exchange functional
• cid: identity of correlation functional
• xcid: identity of exchange-correlation functional

Build

If AVX2 / AVX512 instruction set is supported, please modify the compilation options and compiler directives (e.g. align_size) manually.
If DFT-D4 dispersion correction is involved, make sure you have install DFT-D4.

Linux (Debian as example)

Deploy build tools by root or sudo user:

sudo apt install ninja-build
sudo apt install cmake

Install Intel oneAPI HPC Toolkit, and append the following to ~/.bashrc:

export ONEAPI_ROOT="/path/to/oneapi"
export PATH="$ONEAPI_ROOT/compiler/latest/bin:$PATH"

Download the stable release of Libxc and build it (in oneAPI environment) by:

cmake -S . \
-B build \
-G Ninja \
-DCMAKE_BUILD_TYPE=Release \
-DENABLE_FORTRAN=ON \
-DBUILD_TESTING=OFF \
-DBUILD_FPIC=ON \
-DCMAKE_Fortran_COMPILER="${ONEAPI_ROOT}/compiler/latest/bin/ifx" \
-DCMAKE_C_COMPILER="${ONEAPI_ROOT}/compiler/latest/bin/icx" \
-DCMAKE_BUILD_WITH_INSTALL_RPATH=ON \
-DCMAKE_C_FLAGS="${CMAKE_C_FLAGS} -O3 -xCORE-AVX2"

then

cd build && ninja

append to ~/.bashrc after a successful build:

export LIBXC_ROOT="/path/to/libxc"

Change directory to TRESC root and build it by:

chmod +x release.sh && ./release.sh

You need to create an environment with

conda create --name vis2c python=3.10 numpy scipy matplotlib \
qcelemental mayavi traits traitsui pyqt

grant all .py in /vis2c executable permissions and replace the shebangs with the path to the Python interpreter of env vis2c, it will enable you to use scripts in any scenario. Append to ~/.bashrc when everything is done:

export TRESC="/path/to/TRESC"
export PATH="$TRESC/build:$PATH"
export PATH="$TRESC/vis2c:$PATH"
alias TRESC='tshell.sh'
alias tshell='tshell.sh'

Upcoming

• high-performance vectorized code for 2e-integrals
• perturbation calculation;