Efficient Relativistic Density-Matrix Renormalization Group

Subscriber access provided by UNIV OF DURHAM

Quantum Electronic Structure

An efficient relativistic density-matrix renormalization group implementation in a matrix-product formulation Stefano Battaglia, Sebastian Keller, and Stefan Knecht J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01065 • Publication Date (Web): 20 Mar 2018 Downloaded from http://pubs.acs.org on March 23, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

An efficient relativistic density-matrix renormalization group implementation in a matrix-product formulation Stefano Battaglia, Sebastian Keller, and Stefan Knecht∗ ETH Z¨ urich, Laboratorium f¨ ur Physikalische Chemie, Vladimir-Prelog-Weg 2, 8093 Z¨ urich, Switzerland E-mail: [email protected]

Abstract We present an implementation of the relativistic quantum-chemical density matrix renormalization group (DMRG) approach based on a matrix-product formalism. Our approach allows us to optimize matrix product state (MPS) wave functions including a variational description of scalar-relativistic effects and spin–orbit coupling from which we can calculate, for example, first-order electric and magnetic properties in a relativistic framework. While complementing our pilot implementation (S. Knecht et al., J. Chem. Phys., 140, 041101 (2014)) this work exploits all features provided by its underlying non-relativistic DMRG implementation based on an matrix product state and operator formalism. We illustrate the capabilities of our relativistic DMRG approach by studying the ground-state magnetization as well as current density of a paramagnetic f 9 dysprosium complex as a function of the active orbital space employed in the MPS wave function optimization.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

With the advent of the density-matrix renormalization group (DMRG) approach 1–4 in (nonrelativistic) quantum chemistry, 5–17 approximate solutions for a complete-active-space (CAS)type problem within chemical accuracy became computationally feasible on a routine basis albeit the fact that variational optimization merely involves a polynomial number of parameters. The combination with a self-consistent-field orbital optimization ansatz (DMRGSCF) 18–22 makes active orbital spaces accessible that surpass the CASSCF limit by about five to six times. Yet, the construction of a (chemically) meaningful active orbital space poses a fundamental challenge that calls for automation. 23–26 In contrast to early work on ab initio DMRG, 3 recent efforts on the quantum-chemical DMRG algorithm 4 focused on a variational optimization of a special class of ansatz states called matrix product states (MPSs) 27,28 in connection with a matrix-product operator (MPO) representation of the quantum-chemical Hamiltonian. 17,29–32 Combining the principles of quantum mechanics 33 and special relativity 34 into relativistic quantum mechanics, Dirac 35 put forward a (one-electron) theory which constitutes the central element of relativistic quantum chemistry. 36,37 In this framework, deviations in the description of the electron dynamics increase compared to (nonrelativistic) Schr¨odinger quantum mechanics as one considers electrons moving (in the vicinity of a heavy nucleus) at velocities close to the speed of light. Consequently, the differences between results of relativistic quantum calculations employing a finite and an infinite speed of light, respectively — which are for the latter case equivalent to results from conventional nonrelativistic calculations based on the Schr¨odinger equation — can serve as a definition for the chemical concept of “relativistic effects”. The latter are commonly split into kinematic relativistic effects (sometimes also referred to as scalar-relativistic effects) and magnetic effects both of which grow approximately quadratic with the atomic number Z. 38 A scalar-relativistic theory therefore considers only kinematic relativistic effects whereas a fully relativistic theory includes both kinematic 2


Page 2 of 55

Page 3 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


relativistic and magnetic effects. Kinematic relativistic effects reflect changes of the nonrelativistic kinetic energy operator while spin-orbit (SO) interaction, which is the dominant magnetic effect (for heavy elements), originates from a coupling of the electron spin to the induced magnetic field resulting from its orbital motion in the field created by the other charged particles, namely the nuclei and remaining electrons. 39 Encounters of relativistic effects are ubiquitous in the chemistry and physics of heavy elements compounds, see for example Refs. 36, 37 and 40 and references therein. Besides kinematic relativistic effects, SO interactions are also expected to strongly influence the chemical bonding and reactivity in heavy-element containing molecular systems. 41,42 For example, Ruud and co-workers 43 as well as Gaggioli et al. 44 recently highlighted distinguished cases of mercury- and goldcatalyzed reactions where SO coupling effects are driving the main reaction mechanisms by paving the way for catalytic pathways that would otherwise not even have been accessible. We ultimately aim in this work at the calculation of (static) magnetic properties for a prototypical molecular, open-shell lanthanide complex. The calculation of properties such as molecular g-factors and electron-nucleus hyperfine coupling, which are central parameters in electron paramagnetic resonance (EPR) spectroscopy, 45,46 either requires SO-coupled wave functions or assume a particular simple form in a relativistic framework. This includes also, for example, the calculation of magnetization and current densities. 37 We will pay in this work particular attention to the latter two properties since they play a particular role for the eligibility of open-shell lanthanide tags in spin-labeling of protein complexes. 47–51 As open-shell electronic structures are often governed by strong electron correlation effects, their computational description calls for a multiconfigurational ansatz 52,53 which usually set out from a spin-free (non-relativistic or scalar-relativistic) formulation. Here, electron correlation is commonly split into a static contribution and a dynamic contribution with the latter often being treated in a subsequent step by (internally-contracted) multi-reference perturbation or configuration interaction theories. 52,53 Their zeroth-order wave function is typically of CAS-type such as CASSCF or DMRG(-SCF) which are assumed to be able to

3



adequately grasp static electron correlation effects. In this context, additional challenges may be encountered in particular for heavy-element complexes which arise from the large number of (unpaired) valence electrons to be correlated, i.e., the (n − 2)f (n − 1)d ns np manifold of f -elements, as well as the occurrence of near-degeneracies of electronic states. To arrive at SO-coupled eigenstates in such a spin-free setting, requires then to invoke a so-called two-step SO procedure (see for example Refs. 54–58 for recent developments based on CAS- and DMRG-like formulations). Notably, this procedure relies on an additivity of electron correlation and SO effects or a weak polarization of orbitals due to SO interaction, or both. Moreover, it entails the need to take into account a priori a sufficiently large number of spin-free electronic states for the evaluation of the spin-orbit Hamiltonian matrix elements in order to ensure convergence of the SO-coupled eigenstates. Hence, the predictive potential of such two-step approaches is inevitably limited for a finite set of spin-free states. By contrast, a genuine relativistic DMRG model, initially proposed for a “traditional” (i.e. non-MPO based formulation) DMRG model in Ref. 59 and proposed in this work within an efficient MPS/MPO framework of the DMRG algorithm, addresses simultaneously all of the above issues by providing access to large active orbital spaces combined with a variational description of relativistic effects in the orbital basis. The paper is organized as follows: We first briefly review the concept of representing a wave function and an operator in a matrix-product ansatz, respectively, including a short summary of the variational optimization of an MPS wave function in such a model. In Section 3, a relativistic Hamiltonian framework suitable for a variational description of SO coupling is introduced and a detailed account of its implementation in an MPS/MPO setting of the DMRG ansatz is given. In addition, we illustrate the calculation of first-order molecular properties within our relativistic DMRG model. Subsequently, after providing computational details in Section 4.1, we assess in Section 4.2 the numerical performance of our present implementation with respect to the calculation of absolute energies for the sample molecule thallium hydride. In Section 4.3, we demonstrate the applicability of our

4


Page 4 of 55

relativistic MPS/MPO-based DMRG model for the calculation of magnetic properties in a real-space approach at the example of a large Dy(III)-containing molecular complex. Conclusions and an outlook concerning future work based on the presented relativistic DMRG model are summarized in Section 5.

2

Matrix product states and matrix product operators

We briefly introduce in this section the concepts of expressing a quantum state as an MPS and a (Hermitian) operator as an MPO in a non- or scalar-relativistic framework which will constitute the basis for our genuine relativistic DMRG model — including a variational account of SO coupling — discussed in Section 3. Our notation closely follows the presentations of Refs. 30 and 60. A comprehensive review of the DMRG approach in an MPS/MPO formulation can be found in Ref. 4.

2.1

MPS “left-to-right sweep”

left subsystem

active sites

right subsystem

“right-to-left sweep” = orbital (“site”)

Figure 1: Alignment of orbitals (“sites”) of a given active orbital space on a fictitious (as indicated by the connecting · · · ) one-dimensional lattice. The latter is partitioned into a left and right subsystem as well as active sites on which DMRG many-electron basis states are constructed and successively optimized in a variational approach. While traversing (“sweeping”) through the lattice from left to right, the left subsystem is gradually enlarged by the left active site whereas the leftmost site on the right subsystem becomes active. Reaching the end of the lattice initiates a right-to-left sweep and the procedure is reversed.

Consider an arbitrary state |Ψi in a Hilbert space spanned by (N electrons distributed 5


Page 6 of 55

in) L spatial orbitals that are assumed to be ordered in a favorable fashion 61 on a fictitious one-dimensional lattice. In a traditional (CI-like) ansatz, |Ψi is commonly expressed as a linear superposition of occupation number vectors (ONVs) |ki with the expansion coefficients ck1 ...kL |Ψi =

X

ck |ki =

k

X

ck1 ...kL |k1 . . . kL i .

(1)

k1 ,...,kL

Note that, in a non-relativistic formulation, each local space is of dimension four corresponding to the possible local basis states

kl = {|↑↓i , |↑i , |↓i , |0i} ,

(2)

for the l-th spatial orbital. By contrast, in an MPS representation of |Ψi, the CI coefficients l {ck1 ...kL } are written as a product of ml−1 × ml -dimensional matrices M kl = {Makl−1 al }, i.e.,

four matrices corresponding to the local space dimension,

|Ψi =

X

X

k1 L M1a Mak12a2 · · · MakL−1 1 |k1 . . . kL i = 1

k1 ,...,kL a1 ,...,aL−1

X

M k1 M k2 · · · M kL |ki .

(3)

k

Given the {ck1 ...kL } being scalar coefficients, requires that the final contraction of the matrices M kl must yield a scalar number. This further implies that the first and the last matrices in Eq. (3) are in practice 1 × m1 -dimensional row and mL−1 × 1-dimensional column vectors, respectively. The introduction of some maximum dimension m for the matrices M kl , with m commonly referred to as number of renormalized block states, 2 constitutes the fundamental idea that facilitates to approximate a full CI-type wave function (cf. Eq. (1)) to chemical accuracy while merely requiring to optimize a polynomial number of parameters in an MPS wave function ansatz. Moreover, the partitioning of the CI coefficients in products of the M kl matrices given in Eq. (3) is not unique, i.e., |Ψi may be expressed (through successive applications of a

6


Page 7 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


singular value decomposition) in different ways such as in a left-canonical form 4

|Ψi =

X

Ak1 Ak2 · · · AkL |ki ,

(4)

Akl † Akl = I ,

(5)

k

with X kl

where I is the unit matrix and the Akl matrices are referred to as left-normalized MPS tensors. Similarly, the MPS can be written in right-canonical form 4

|Ψi =

X

B k1 B k2 · · · B kL |ki ,

(6)

k

where the B kl matrices are right-normalized, satisfying the condition X

B kl B kl † = I .

(7)

σk

The non-uniqueness, which allows us to express an MPS in different forms originates from the existence of a gauge degree of freedom. 4 Considering for a given MPS two adjacent sets of matrices M kl and M kl+1 with matching bond dimensions m, it can be shown 4 that the MPS is invariant with respect to a right (left) multiplication with matrix X of dimension m × m, M kl → M kl X ∀kl ,

M kl+1 → X −1 M kl+1 ∀kl ,

(8)

provided that X is invertible and non-singular. The gauge freedom may therefore be further exploited to express the MPS of Eq. (3) in mixed-canonical form at sites (orbitals) {l, l + 1} 4

|Ψi =

X

X

k

a1 ,··· ,aL−1

kL kl+2 l kl+1 Ak1a1 1 · · · Akal−1 Makl−1 ,al+1 Bal+1 al+2 · · · BaL−1 1 |ki , l−2 al−1

7


(9)


Page 8 of 55

k

k

l+2 kl−1 = {Bal+1 where the MPS tensors Akl−1 = {Aal−1 al+2 } are left- and rightl−2 al−1 } and B

kk

l l+1 normalized, respectively, as shown in Eqs. (5) and (7) and Mal−1 ,al+1 is a two-site MPS

tensor l kl+1 Makl−1 ,al+1 =

X

kl+1 l Makl−1 ,al Mal ,al+1 .

(10)

al

The mixed-canonical form will be a central element of the two-site DMRG optimization algorithm outlined in Section 2.3.

2.2

MPO

Exploiting the matrix-product formulation introduced in the previous section for operators, c can be expressed in MPO form as an N -electron operator W

c = W

X

X

X

k k0

k k0

k k0

L L 0 0 W1b11 1 Wb12b22 · · · WbL−1 1 |k1 . . . kL i hk1 . . . kL |

0 b ,...,b k1 ,...,kL k10 ,...,kL 1 L−1

≡

X

wkk0 |ki hk0 | ,

(11)

kk0

with the incoming and outgoing physical states kl and kl0 and the virtual indices bl−1 and bl . Rearranging the summations and changing the order of contraction in Eq. (11) leads to 30,31 X

c= W

1 W1b · · · Wbll−1 bl · · · WbLL−1 1 , 1

(12)

b1 ,...,bL−1

where the entries of the

n

Wbll−1 bl

o

matrices comprise only local creation (annihilation) oper-

ators acting on the l-th orbital. For example, consider the operator a ˜†↑l which can be written as a linear combination of the local basis states a ˜†↑l = |↑↓i h↓| + |↑i h0| .

(13)

Hence, Eq. (13) implies that the matrix representation of a ˜†↑l corresponds to a (4 × 4)dimensional matrix with two non-zero entries equal to one. Likewise, similar considerations 8


Page 9 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hold for the remaining local operators while the matrix representations of the corresponding local, relativistic operators of our relativistic DMRG model will be discussed in Section 3.2. For completeness, we conclude this section by considering the action of the local Hamilˆ given in MPO form on an MPS state in mixed-canonical form (cf. Eq. (9)) at sites tonian H {l, l + 1} which can be written as 4

ˆ |Ψi = H

X

al−1 ,a0l−1

X

Lbl−1

0 kl kl+1 ,k0 kl+1

Wbl−1 ,bl+1l

al+1 ,a0l+1

Rbl+1

k0 ,k0

l l+1 Mal−1 ,al+1 |al−1 iA |kl i|kl+1 i|al+1 iB , (14)

bl−1 ,bl a0l−1 ,kl0 ,a0l

with the left and right basis states given by

|al−1 iA =

X

X

Ak1a1 1 · · · Akal−1 |k1 , · · · , kl−1 i , l−2 al−1

(15)

kL l+2 Bakl+1 al+2 · · · BaL−1 1 |kl+1 , · · · , kL i .

(16)

k1 ,··· ,kl−1 a1 ,··· ,al−1

and X

X

|al+1 iB =

kl+1 ,··· ,kL al+1 ,··· ,aL

Here, we introduced the left and right boundaries as 4,30 

 al−1 ,a0l−1

Lbl−1

=

X

k1 ,k10 1,b1

X

1∗ Ak1,a W 1

 {ai ,bi ,a0i ;il+1}

0

l+1

0 kl+1 ,kl+1

 ×

0

kl+2 ,kl+2 kl+2 l+2 ∗ ··· Bakl+1 ,al+2 Wbl+1 ,bl+2 Ba0 ,a0 l+2

 X

k ,k0

k0

L L L L∗ BakL−1 ,a1 WbL−1 ,b1 Ba0

0 L−1 ,a1

0 kL ,kL

9

 ,


(18)


Page 10 of 55

as well as the two-site MPO tensor 4,30 0 kl kl+1 k0 kl+1

Wbl−1 ,bl+1l

=

X

k k0

kl+1 k0

l l l+1 Wbl−1 . ,bl Wbl ,bl+1

(19)

bl

2.3

Variational optimization of an MPS

ˆ expressed in MPO form, the In a variational optimization of |Ψi with the Hamiltonian H D E ˆ objective is to minimize the energy expectation value Ψ|H|Ψ with respect to the entries of the MPS tensors under the constraint that the wave function is normalized, i.e., hΨ|Ψi = 1. To this end, a Lagrangian multiplier λ is introduced to ensure normalization. The optimization of the MPS then corresponds to find the extremum of the Lagrangian 62

ˆ |Ψi − λ (hΨ| Ψi − 1) . L = hΨ| H

(20)

The optimization parameters are the matrices M kl of the MPS. To solve this highly nonlinear optimization problem, the same idea as in the original (non-MPO) DMRG algorithm l of one matrix is adopted: iterating through the lattice and optimizing the entries Makl−1,a l

(single-site variant) or of the two-site MPS tensor

kl ,kl+1 Mal−1 ,al+1

(two-site variant) at a time,

while keeping all the others fixed. This approach not only greatly simplifies the optimization problem to one of lower complexity but also ensures to find a better approximation to the true ground state in a variational sense. In the following, we consider the two-site variant (corresponding to treating simultaneously two sites as “active” as indicated by a red color code in Figure 1). Taking the derivative k ,k

l l+1 with respect to the complex conjugate of Mal−1 ,al+1

∂ kl ,kl+1 ∗ ∂Mal−1 ,al+1

ˆ |Ψi − λ hΨ| Ψi) = 0 , (hΨ| H

10


(21)

Page 11 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


yields 4 X

X

al−1 ,a0l−1

Lbl−1

0 kl kl+1 ,k0 kl+1

Wbl−1 ,bl+1l

a0

l+1 Rbl+1

,al+1

0 k0 kl+1 0 l−1 ,al+1

Ma0l

0 a0l−1 a0l ,bl−1 bl+1 kl0 kl+1

l kl+1 = λMakl−1 ,al+1 ,

(22)

where we assumed that the left and right boundaries (cf. Eqs. (17) and (18)) were calculated from left- and right-normalized MPS tensors, respectively. The resulting Eq. (22) can be recast in a matrix eigenvalue equation 4,27

Hv − λv = 0 ,

(23)

with the local Hamiltonian matrix H at sites {l, l + 1} after reshaping given by

0 H(kl,l+1 al−1 al+1 ),(kl,l+1 a0l−1 a0l+1 ) =

X

al−1 ,a0l−1

Lbl−1

0 kl kl+1 ,k0 kl+1

Wbl−1 ,bl+1l

a0

l+1 Rbl+1

,al+1

,

(24)

bl−1 ,bl+1

and the vector v collecting l kl+1 vkl kl+1 al−1 al+1 = Makl−1 ,al+1 .

(25)

Since we are often interested in only a few of the lowest eigenvalues λ, Eq. (23) is best solved by an iterative eigensolver such as the Jacobi-Davidson procedure. For example, having obtained the lowest eigenvalue λ0 and the corresponding eigenvector vk0l kl+1 al−1 al+1 , the kk

l l+1 latter can be reshaped back to Mal−1 ,al+1 which is then subject to a left- or right-normalization

k0

into Akall−1 al or Ball+1 al+1 by a singular value decomposition (discarding the 3m smallest singular values) in order to maintain the desired normalization structure and the dimensionality of the MPS tensors. Given the optimized MPS tensors for sites l and l + 1, the complete algorithm now sweeps sequentially forth and back from left-to-right and right-to-left through the lattice (cf. Figure 1) consisting of the L spatial orbitals ordered in a suitable form. At each step l the MPS tensors {Makl−1 al } are optimized until convergence is reached.

11



3 3.1

Page 12 of 55

Formulation of a relativistic DMRG model Hamiltonian framework

In accordance with the pilot relativistic DMRG implementation presented in Ref. 59, our present formulation of a relativistic DMRG model in a matrix-product ansatz is based on the electronic, time-independent, external magnetic-field-free Dirac-Coulomb(-Breit) Hamiltonian (note that any two-component Hamiltonian is equally well applicable) which encompasses all important relativistic “corrections” 36,37,40 relevant for chemistry. Invoking the no-pair approximation, this operator can be written in second quantization as, 36,63

ˆ = H

X pq

ˆ+ hpq X pq

1 + + ˆ + hp¯q X ˆ hp¯q X + p¯q p¯ q 2

1 X (pq|rs) xˆ++ pq|rs) xˆ++ q |rs) xˆ++ pq,rs + (¯ p¯q,rs + (p¯ p¯ q ,rs 2 pqrs 1X + (¯ pq|r¯ s) xˆ++ p¯q,r¯ s 4 pqrs 1 X ++ (¯ pq|¯ rs) xˆ++ + (p¯ q |r¯ s ) x ˆ + p¯q,¯ rs p¯ q ,r¯ s 8 pqrs 1 X −− −− − (pαq|rαs) xˆ−− + (¯ p αq|rαs) x ˆ + (pα¯ q |rαs) x ˆ pq,rs p¯q,rs p¯ q ,rs 2 pqrs 1X (¯ pαq|rα¯ s) xˆ−− − p¯q,r¯ s 4 pqrs 1 X − (¯ pαq|¯ rαs) xˆ−− q |rα¯ s) xˆ−− p¯q,¯ rs + (pα¯ p¯ q ,r¯ s , 8 pqrs

+

(26)

where the last three rows on the right-hand side of Eq. (26) comprise the two-electron terms originating from the Breit (Gaunt) interaction and the summation indices p, q, r, s strictly refer to positive-energy spinors. Moreover, following earlier works, 59,63–66 we assume in Eq. (26) a Kramers pair basis {ϕl , ϕ¯l } built from pairs of fermion spinor functions {ϕl } and {ϕ¯l }, respectively. This implies that the summations in Eq. (26) run over Kramers pairs and not single spinors. According to Kramers theorem, which holds in the absence

12


Page 13 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of an external magnetic field, states of a single fermion function {ϕl } are (at least) doubly degenerate where each state is related to its energy-degenerate partner by time reversal, 36,37 i.e., ˆ l = ϕ¯l , Kϕ

(27)

ˆ is the time-reversal operator. In addition, we employ in Eq. (26) Kramers singlewhere K s 1 ,s2 ˆ pq (with s1 , s2 = ±), 67 respectively. (with s = ±) and double replacement operators xˆspq,rs X

They can be expressed in terms of creation and annihilation operators as 36,63,67,68 ˆ s = a† aq + sa†q¯ap¯ X pq p

ˆ s = a†p¯aq − sa†q¯ap X p¯q

ˆ s = a† aq¯ − sa† ap¯ , X p¯ q p q

(28)

and † † † † 1 ,s2 2 ,s1 ˆ s1 X ˆ s2 xˆspq,rs =X ˆsrs,pq , pq rs − δrq ap as − s1 δrp¯aq¯as − s2 δs¯q ap ar¯ − s1 s2 δp¯s¯aq¯ar¯ = x

(29)

where the sign indices s and s1 , s2 in Eqs. (28) and (29) indicate the symmetry of the operators under time reversal and Hermitian conjugation, respectively. The remaining double replacement operators appearing in Eq. (26) follow from Eq. (29)] through the application of time-reversal symmetry to the creation (annihilation) operators. 36 An efficient account of double-group symmetry is a crucial aspect of our relativistic DMRG formulation. In contrast to a spin-free approach, where spin and spatial degrees of freedom can be considered independently, SO coupling introduced through the relativistic molecular Hamiltonian (cf. Eq.(26)), entails a coupling of both degrees of freedom. Hence, rather than having to deal with simple spatial point group symmetry we have to consider the corresponding double groups. In contrast to the single groups, these groups comprise extra irreps originating from an introduction of a rotation 2π about an arbitrary axis. 36,37,69 These additional irreps are commonly referred to as fermion irreps and are spanned by spinors whereas the “regular” ones are called boson irreps spanned, for example, by spinor products and operators. Following the symmetry handling of our pilot relativistic DMRG implementation, 59 we 13



also take advantage of a quaternion symmetry scheme introduced by Saue and Jensen 70 in the DIRAC program package 71 to which our DMRG software QCMaquis is interfaced. In addition to the Abelian subgroups of the binary double D∗2h , we introduced in this work the double groups C∗16h and C∗32v for linear molecules with and without inversion symmetry as closed, finite Abelian approximations to the full groups D∗∞h and C∗∞v , respectively, along the lines presented in Ref. 72 and exploited in Ref. 73 for a relativistic linear-response coupled cluster approach. Exploiting a quaternion symmetry scheme allows us to resort to either real (for the double groups D∗2h , D∗2 and C∗2v ), complex (C∗2h , C∗2 and C∗s ) or quaternion algebra (C∗i and C∗1 ) in the construction of a matrix representation of the Hamiltonian operator (cf. Eq.(26)) in a finite Kramers-paired spinor basis and the subsequent (self-consistent) solution of the corresponding eigenvalue problem, 36,70 respectively. Moreover, working in a Kramers-paired spinor basis, entails that, by symmetry, all matrix elements of a timereversal symmetric one-electron operator with an odd number of barred spinor indices are zero. 70 Furthermore, for the two-electron integrals, we exploit the so-called (NZ,3) matrix representation 64,74 where the number of nonzero rows for a given binary double group corresponds to the rank NZ of the matrix. For example, as shown in Eq. (26), a two-electron integral (pq|rs) may generally comprise a mixed number of barred and unbarred spinors. Following the (NZ,3) classification discussed in detail in Ref. 64, for real- and complex-valued double groups, labeled as NZ=1 and NZ=2, respectively, only integrals with an even number (nbarred = 0, 2, 4) of barred spinors are non-vanishing and are either real-valued (NZ=1) or complex-valued (NZ=2). By contrast, for quaternion-valued double groups with NZ=4, all integrals are in general non-zero and complex-valued. In summary, for a given double group we exploit a quaternion algebra scheme solely in the DIRAC program package to reduce the number of non-zero one- and two-electron (complex-valued) integrals in the Hamiltonian which in turn serve as input coefficients for the MPO construction of the subsequent DMRG step that is implemented in the QCMaquis software package. Consequently, the latter step does not involve a formulation that makes use of a quaternion algebra scheme.

14


Page 14 of 55

Page 15 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2

Relativistic DMRG within an MPS/MPO framework

In the previous section 3.1, we briefly introduced a Hamiltonian framework suitable for a relativistic DMRG model. In the following we will elaborate on details of a relativistic DMRG model which exploits an MPS wave function and MPO Hamiltonian representation. The basic concepts are identical to those discussed in Section 2 for a non- or scalar-relativistic framework. Hence, we will in particular focus on differences that originate from the underlying use of a relativistic Hamiltonian model. We first provide a brief introduction to the implementation of (double group) symmetry in our MPS/MPO framework. Subsequently, details concerning the actual implementation of a relativistic DMRG model in our QCMaquis software packages are presented. 3.2.1

Symmetry

As indicated in Section 3, exploiting symmetries in a (relativistic) DMRG algorithm is a key element to enhance both accuracy and computational efficiency (see for example also Ref. 31 for the spin-symmetry adaptation of the MPS/MPO representation of wave function and operators in our QCMaquis software package). At present, we only consider the application of (binary) Abelian double groups for our relativistic DMRG model including in addition closed, finite Abelian double group approximations for the linear (non-Abelian) double groups C∗∞v and D∗∞h . For a brief introduction on the use of double-group symmetry in quantum chemistry we refer the interested reader to the textbook by Dyall and Fægri. 36 Concerning the scope of applications for our relativistic DMRG model towards the study of large transition-metal and f -element complexes that often exhibit either low- or even no point-group symmetry, such a restriction will, however, be less critical with respect to the computational efficiency of our resulting implementation. To illustrate the account of symmetry in our relativistic DMRG model, we consider an

15



Page 16 of 55

ˆ that commutes with the Hamiltonian of Eq. (26), i.e., operator Q h

i ˆ Q ˆ =0. H,

(30)

Hence, the eigenfunctions |Ψi of the Hamiltonian can always be chosen as eigenfunctions of ˆ with Q ˆ |Ψi = q |Ψi . Q

(31)

Hence, it is a sine qua non in the DMRG optimization algorithm that in each step any local operation transforms the basis according to the group it belongs to such that any element of the basis remains an element of the group after the transformation. To further ˆ 1 and Q ˆ 2 which both commute with H, ˆ illustrate this concept, we consider two operators Q e.g., h i ˆ ˆ H, Q1 = 0, i h ˆ Q ˆ 2 = 0, H,

(32) (33)

with the corresponding eigenvalues (in the following denoted as quantum numbers) given by

ˆ 1 |Ψi = q |Ψi , Q 1

(34)

ˆ 2 |Ψi = q |Ψi . Q 2

(35)

This allows us to label the eigenstate |Ψi according to the quantum numbers given by Eqs. (34) and (35), respectively,

|Ψi → :j ,

(36)

where j is an integer counting the possible number of realizations for this state with quantum numbers q1 and q2 during the DMRG optimization procedure. Assume, for example, in a 16


Page 17 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DMRG sweep starting with the left subsystem composed by only the first site and with the symmetries of the total system defined by the particle number N and a given double group ˆ1 = N ˆ and Q ˆ 2 = G, ˆ where G ˆ is any allowed operation of the double symmetry. Hence, Q group. The quantum numbers are therefore relabeled as q1 → N and q2 → g. At each site, the local space has dimension two, with one state corresponding to an empty spinor (labeled :1) and the other state corresponding to an occupied spinor (labeled :1), where gl refers to the irreducible representation of the spinor ϕl on site l. Note that the lattice is entered from the left with the vacuum state |a0 i = :1. After optimization of the first two sites, the first site is merged into the left subsystem block consisting of the vacuum state such that the new basis states are given by

{|a1 i} = |a0 i ⊗ {|k1 i}.

(37)

The two states defined on the first site are both eigenstates of the symmetry operators and can be labeled accordingly as

|k1 = 0i = :1 ,

(38)

|k1 = 1i = :1 .

(39)

Since the state |a0 i is an eigenstate of the symmetry operators, the tensor product will also be an eigenstate. This condition can always be enforced by the following equalities

ˆ |a0 i ∗ N ˆ |k1 i = N ˆ |a1 i , N

(40)

ˆ |a0 i ∗ G ˆ |k1 i = G ˆ |a1 i , G

(41)

where ∗ denotes the group operation. Including the second site into the left subsystem results

17



Page 18 of 55

again in a basis of eigenstates of the symmetry operations,

ˆ |a1 i ∗ N ˆ |k2 i = N ˆ |a2 i , N

(42)

ˆ |a1 i ∗ G ˆ |k2 i = G ˆ |a2 i , G

(43)

where, by construction, the states {|a1 i} are eigenstates of the symmetry operators. The new states of the enlarged left block comprising the first two sites are therefore eigenstates ˆ and G. ˆ This procedure can be iterated until the end of the lattice is reached which will of N lead to a basis set of many-particle states that satisfy the global symmetry constraints. The above considerations hold for any type of parametrization of the quantum states, hence, also for an MPS wave function representation. To illustrate the latter, we consider a lattice of L spinors and impose as symmetry constraints a total number of particles of N = 4 and the totally symmetric irreducible representation for the target quantum state, i.e., g=0. Following the procedure outlined above for all sites 4 < l ≤ L all states arising from the tensor product of the left subsystem with the local l-th site basis that would lead at the final site L to N > 4 will be automatically discarded because they do not satisfy our imposed symmetry constraints for the target state. The same restrictions hold for all intermediate states that satisfy N = 4 but, combined in all possible ways with the corresponding local states up to site L, would yield g 6= 0 because the resulting quantum state would no longer transform as the totally symmetric irreducible representation. Hence, exploiting symmetries leads to a considerable reduction of the number of possible states that need to be taken into account which in turn lowers the associated numerical effort of an MPS wave function optimization.

3.2.2

Implementation aspects

In this section we discuss selected important aspects of the actual implementation of our relativistic Hamiltonian model within the existing (nonrelativistic) framework of QCMaquis,

18


Page 19 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


while further details on the building blocks of QCMaquis can be found in Refs. 30,31,60,75. Figure 2 illustrates the main classes that constitute the essential building blocks for an MPS wave function optimization based on a given Hamiltonian model within the framework of QCMaquis. dmrg_measure

dmrg_sim

MPS Sim MPO

Lattice

Model

Figure 2: UML class diagram outlining the logic of the QCMaquis program. The central objects in QCMaquis are the MPS tensors M kl , introduced in Section 2.1, which are represented by the class MPS in Figure 2. By taking into account the symmetry considerations outlined in the previous section, we can associate each MPS tensor a-index 60 l of Makl−1 al (cf. Eq. (3)) with a quantum number (or charge)

ql =< Nl , gl > ,

(44)

l which motivates an extended notation for each MPS tensor Mqkl−1 al−1 ;ql al where the a-indices

are supplemented with the associated quantum numbers ql−1 and ql , respectively. Because of the symmetry constraint ql ∈ ql−1 ⊗ kl ,

19


(45)


the MPS tensor partitions into symmetry blocks as indicated by the extended index ql−1 al−1 ; ql al . Hence, an MPS tensor can uniquely be characterized by the quantum numbers ql−1 and ql and kl which are also referred to as left, right, and physical index, respectively. 60 The block structure of the MPS tensor consists of one block matrix for each of the local basis states of kl which will hold in the case of Abelian double-group symmetry as is generally assumed in the present work. For an efficient storage, the physical index is fused in QCMaquis either to the left or right index resulting in a rank-two tensor, i.e., a matrix, as depicted in Figure 3.

Mkla

q a q

sha1_base64="jZG8y0x+T479RCu/s4LQfMy08Tw=">AAAI3nicjVVbb9s2FFbTJbOzW9M97oVYYqABJE3SGsdG4aHYXvawAdmwrAVsz6AoyiZMXUbSTV1Bz3sa9rj9tf2ZYYeUnUm23JVAYvLcvnO+c0SGOWdSed7fD44evnd88n6ne/rBhx99/Mmjs8c/y2wlCL0lGc/EyxBLyllKbxVTnL7MBcVJyOmLcPmN1r94RYVkWfqTWud0muB5ymJGsALR7Oz4z+4kpHOWFoot3+SMqJWg5VgSzOnIn552u5NI4Ds0jhnnU/TEs71LRJggnMLBDS6fbU3GALNGasHIsrIDJXKcatu/RGkW0TEOs1cU1JLNL1FxsZzxi9JEqGux0k6+C0jFzWItIVmOZI4JLQ+huYHOy6C5fb018QSbLxSoRYX2axNto9VogKWdih+16O1QTg3LqYNxGmss+KmwCu74ZR2vMtBwju8+NXjfgagG10PaBuWYCRrtUe8EthO8E/lO4Po2gAyvNnkGkGd1fqtH4HoND3NuLWCrrpWKJlgqtNPRkPLsDu07BtviYXaZPMi17w5MzSYlTZo+1Ju311QDeM+yv8WqGvsfWA+ZAIeofmemd4j+X553aG5juTaWdZYdYLaiWPNd7pO879dWeEtOdZIbHFdta53l3b5WMdra2kNYiOxO3mMbWNv5qvqWvrT1v/sOX2v6rttnzrTT3yA4unMsnQMG2qzerofnDkz4gSFwouhrZS7MAlpeFt//UgCl5ayYxFmm0kxRyd7QAlfFbmp+hjDQDjNWltv6d4to1tBSQuOe2ZRg2lKroTuhadS4gmePzj3XD/yr/hDBAA76wfUVbJ4OPX/YR3BhmXVubdbN7Ozoi0mUkVVCU0U4lnLse7maFlgoBmNcnk5WksJVs8RzOoZtihMqp4VhpEQ9kEQozgT8pQoZad2jwImU6yQEywSrhdzVaWGrTtdUdhuifB7nPFOybJUqHO5lq+LBtGBpvlI0JVWy8YojlSH9pqEIPmKi+Bo2mAgG9SKywAITBS/fqQml0+AsFFisizyTTL9+QL0dM2XD80IAUJtItea0iGG4w+x1iUZorCPjdM6pnbCUJasE3bFILUaBe0USG22FC6o7OvJBpucMEegChTmzkZ6XUcihlGkdQ65yKvZBDjrbEZYLGjVimC/LBLifx7pae8x2baowxhRGbDtH6PDmNnCHrv9DcP78682sdazPrM+tJ5ZvXVvPrW+tG+vWIsf/nFycOCduB3d+6/ze+aMyPXqw8fnUaqzOX/8Cnb6wNg==

l-1 l-1 l l

Efficient Relativistic Density-Matrix Renormalization Group

Recommend Documents