# Relativistic Quantum Calculations to Understand the Contribution of f-type
Atomic Orbitals and Chemical Bonding of Actinides with Organic Ligands

Andy D. Zapata-Escobar  [  Srimanta Pakhira  Department of Physics, Department
of Metallurgy Engineering and Materials Science (MEMS), Centre for Advanced
Electronics (CAE), Indian Institute of Technology Indore (IIT Indore), Simrol,
Khandwa Road, Indore, 453552, Madhya Pradesh (M.P.), India.  Joaquin Barroso-
Flores  Instituto de Química, Universidad Nacional Autónoma de México (UNAM),
Circuito Exterior SN, Ciudad Universitaria, Coyoacán CP 04510 Ciudad de
México, México  Gustavo A. Aucar  Physics Department, Natural and Exact
Science Faculty, Northeastern University of Argentina<EMAIL_ADDRESS>(G.A.A.) Jose L. Mendoza-Cortes  Department of Physics, Scientific Computing,
Material Sciences and Engineering, High-Performance Material Institute,
Condensed Matter - High Magnetic Field National Lab, Florida State University,
Tallahassee FL, 32310, USA<EMAIL_ADDRESS>(J.L.M-C.)

###### Abstract

The nuclear waste problem is one of the main interests of the rare earth and
actinide elements chemistry. Studies of Actinide-containing compounds are at
the frontier of the applications of current theoretical methods due to the
need to consider relativistic effects and approximations to the Dirac equation
in them. Here, we employ four-component relativistic quantum calculations and
scalar approximations to understand the contribution of f-type atomic orbitals
in the chemical bonding of actinides (Ac) to organic ligands. We studied the
relativistic quantum structure of an isostructural family made of Plutonium
(Pu), Americium (Am), Californium (Cf), and Berkelium (Bk) atoms with the
redox-active model ligand; DOPO (2,4,6,8-tetra-tert-
butyl-1-oxo-1H-phenoxazin-9-olate). Crystallographic structures were available
to validate our calculations for all mentioned elements except for Cf. In
short, state-of-the-art relativistic calculations were performed at different
levels of theory to investigate the relativistic effects and electron
correlations on geometrical structures and bonding energies of $Ac$-DOPO3
complexes ($Ac$=Pu, Am, Cf, Bk) : 1) the scalar relativistic zeroth order
regular approximation (ZORA) within the hybrid density functional theory (DFT)
and 2) the four-component Dirac equation with the Dirac-Hartree-Fock (4c-DHF)
and Lévy-Leblond (LL) Hamiltonians. We show that scalar DFT-ZORA could be used
as an efficient theoretical approximation to first approximate the geometry
and electronic properties of actinides which are difficult to synthesize or
characterize; but knowing that the higher levels of theory, like the 4c-DHF,
gives closer results to experiments than the scalar DFT-ZORA. We also
performed spin-free calculations of geometric parameters for the Americium and
Berkelium compounds. To the best of our knowledge, this is the first time that
such kind of large actinide compounds (the largest contains 67 atoms and 421
electrons) are studied with highly accurate four-component methods (all-
electron calculations with 6131 basis functions for the largest compound). We
show that relativistic effects play a key role in the contribution of f-type
atomic orbitals to the frontier orbitals of $Ac$-DOPO3 complexes. The analysis
of results obtained applying different theoretical schemes to calculate
bonding energies are also given.

###### keywords:

Magnetic shielding LRESC Relativistic effects …

unne]Institute for modeling and innovative technology, IMIT (CONICET-UNNE),
Corrientes, Argentina [imit]Institute for modeling and innovative technology,
IMIT (CONICET-UNNE), Corrientes, Argentina [imit]Institute for modeling and
innovative technology, IMIT (CONICET-UNNE), Corrientes, Argentina Department
of Chemical and Biomedical Engineering, FAMU-FSU Joint College of Engineering,
Tallahassee, FL 32310, USA. Current address: Department of Chemical
Engineering & Materials Science, Michigan State University, East Lansing,
Michigan 48824, USA.

##  1 Introduction

Our understanding of how chemistry and physical properties evolve across rows
and down groups in the modern periodic table has been at the forefront of
chemistry for the last 150 years. The latest additions the so-called actinide
elements, which are the heaviest elements, are at the frontier of our current
chemical and physical understanding.1, 2, 3, 4, 5 The radioactive actinide
elements residing underneath the main body of the periodic table conjure
images of mushroom clouds, glowing test tubes, and superheroes awakening.
While there is some truth to this view, as chemists working beyond the edge of
nuclear stability, we have peered into the abyss to find elements whose
realities are so complex and mystifying that we do not need to invent fanciful
tales to remain a captive audience. Our fundamental knowledge and
understanding of actinide bonding trends and coordination modes is
significantly limited compared to other elemental series in the periodic
table. Even though there are several and thorough studies involving actinide
elements,6-10 there are no studies that use state-of-the-art four-component
relativistic methods that permits to grasp on the contribution of f-type
atomic orbitals (AOs) in the electronic structures and chemical bonding of
actinides (Ac) to medium-size organic ligands. This is needed to understand
the changes in interactions and chemical bonds, which may impact speciation
and separation selectivity to different heavy elements.1, 2, 3, 4

Some chemical reactions are clearly explained by considering relativistic
effects on valence electrons of heavy-atoms, e.g. the electromotive force of
the cells of lead-batteries, where the relativistic effects on Pb atoms allows
for the battery to work.5 There are also few other examples which are
explained by relativistic effects: e.g., the yellow color of gold, the yellow
color of hexachloroplumbate(IV), the catalytic methane activation by Pt+, to
name a few.6 From these grounds, we have undertaken a research program whose
goals are oriented to learn more about the physics that underlies the
structure and electronic properties of actinides bonding, where one must
consider relativity. Here, we present our first results of this research
program, by studying actinide containing compounds having Pu, Am, Cf and Bk
and an organic ligand model. Accordingly, we show how the 5f electrons of
these actinides behaves differently when relativity is taken into account,
thus changing the frontier molecular orbitals of complexes. More specifically,
we estimate relativistic effects at different levels of theory: ZORA–HF,
ZORA–DFT, ZORA–DFT-D3, and LL and DHF. More specifically, a) the spin-orbit
(so) in addition to the scalar (sr) relativistic zeroth order regular
approximation in hybrid density functional theory (ZORA–DFT)7, 8 and without
DFT (ZORA-HF), and b) four-component (4c) Dirac equation with the Dirac-
Hartree-Fock (DHF), Lévy-Leblond (LL) and spin-free (SF) Hamiltonians.9 In
this manner, it was possible to learn about the influence of both,
relativistic and electron correlation effects on the geometry, the energy
involved in the bonding process and the change of 5f and 6d AOs in the
population of molecular orbitals. To the best of our knowledge, this is the
first time that such kind of large actinide compounds are studied with highly
accurate four-component methods.

The research for carbon-free energy has generated a tremendous interest in
alternative energy sources. Nuclear energy is a promising alternative because
it offers opportunities for low cost and high efficiency. However, it
generates a significant amount of radioactive waste, which we are not able to
process or reduce manly because of our lack of understanding of bonding in
heavy elements. The proper disposal of nuclear waste is still challenging and
one of the most difficult kinds of waste to handle because it is highly
hazardous and harmful for humans and the biosphere even for small
concentrations, and especially if they are radioactive. Advanced actinide
separation strategies reduce hazards and costs of managing radioactive waste.
Among them liquid fluoride thorium reactor and other molten salt reactor (MSR)
have been designed to use the thorium fuel cycle with a fluoride-based,
molten, liquid salt for fuel. These reactors are passively safe, produce less
radioactive waste, offer higher fuel burn up efficiency, and do not produce
industrial quantities of fissile 239Pu.1 Similarly, the development of
vitrification using borate chemistry was inspired, in part, from evidence that
a significant boron concentration exists in the brine solutions found at the
Waste Isolation Pilot Plant (WIPP) in Carlsbad, New Mexico.10, 11

Actinide elements are of considerable technological importance. In condensed
matter, as they may be high temperature superconductors or centers of lasing
activity.3, 12 Recently, various studies have been performed in lanthanides
and actinide elements to understand the f-orbitals expansion and contractions
when loosing or gaining electrons.3, 4, 13, 14 Therefore, the further
understanding of the chemical similarities and differences between the
trivalent lanthanides and the actinides analogues is of immense interest. One
indication of these advances is the rapid increase in f-block element
containing compounds for which high-resolution structures exist, while their
properties have been investigated both experimentally and theoretically during
the last decade.2, 3 This is promising because the f-electrons remain quite
localized in going from the atomic to the condensed state, thus, a lot of
knowledge gained from atoms might be transferable to the condensed state.3,
13, 14, 15 Therefore, understanding f-block elements and trends is essential,
yet challenging because of the intrinsic radioactivity of most of the actinide
elements.4, 16 However, this lack of knowledge poses challenges that prevents
advances in nuclear energy, radioactive waste management, separation, storage,
and next-generation f-electrons materials.16

There has been an increasing interest in fundamental research about actinide
elements since the 20th century. Understanding the chemical bonding, chemical
reactivity, interaction, and speciation of hazardous heavy actinide elements
in aqueous solutions are a first step toward modeling their behavior, both in
the laboratory and under natural conditions. Although the chemistry of the
actinides has been extensively explored for over a century, its theoretical
grounds and basic understandings still lags much compared to that of most
elements in the periodic table.17 The ability of actinides to form chemical
bonds and the role of 5f and 6d AOs in the contribution to molecular orbitals
(MO) are still subjects of investigation.18, 19 We are still developing our
understanding of the role of 5f- and 6d-shell electrons from actinides
elements when interacting with other atoms or molecules.

Recently, Silver et al. experimentally characterized and synthesized the
berkelium coordination complex, Bk(III)tris(dipicolinate), and showed that it
had formed a chemically distinct Bk(III) borate material.20 A recent study has
synthesized nonaqueous isostructural family of f-element compounds (Ce, Nd,
Sm, Gd; Am, Bk, Cf) of the redox-active dioxophenoxazine ligand (DOPO =
2,4,6,8-tetra-tert-butyl-1-oxo-1H-phenoxazin-9-olate) for the first time.17
Their electronic structure was further studied by applying the zeroth order
regular approximation (ZORA) with hybrid DFT as well as CASSCF calculations to
predict and investigate the structure, oxidation states, and electronic
properties of both the lanthanides and actinides complexes in detail.17 More
recently, Albrecht-Schmitt et al. showed that the complexation of
berkelium(III) by carbonate results in spontaneous oxidation to berkelium(IV)
and that multiple species can be present in solution.21 They explored the
chemistry of Bk(IV) carbonate and carbonate-hydroxide complexes based on
theoretical comparisons with spectroscopic data previously reported for Bk(IV)
carbonate solutions.21 They performed molecular orbital calculations to try to
obtain a better understanding of the nature of their chemical bonding using
DFT at the level B3PW91. Brenner et al. synthesized and analyzed the
structure, and solid-state UV–vis–NIR spectroscopy of three new f-element
squarates, M2(C4O4)3(H2O)4, where M = Eu, Am, Cf and they showed that Cf has a
3+ ionic state with the nine–coordinate ionic radius of Cf3+ about 1.127 $\pm$
0.003 Å.4 Galley et al. synthesized the Cf(DOPO)2(pyridine)(NO3) complex and
they carried out a computational analysis of the electronic structure of the
complex to investigate the ligand-CfIII interactions.17 Goodwin et al.
experimentally synthesized Am3+ organometallic complex [Am(C5Me4H)3] and
applied DFT and ab initio wave function theory calculations to explore the
interaction and electronic structure of the complex.22

In the next sections we shall first explain the methodology which we have
implemented to calculate the geometry and bonding energies of the molecular
compounds. We have considered two different levels of theory: 1) ZORA23, 24,
25, 26, 27, 28 and 2) 4c-DHF and LL, with different schemes of basis sets.26,
27, 28 The main idea is to show how important relativistic effects are when
one is interested to reproduce those properties. Afterwards we shall focus on
one of the main aims of this work, namely, to learn about the dependence of
the participation of p-, d- and f- AOs of the actinides atoms on the frontier
electronic structure of the just mentioned family of compounds: $Ac$-DOPO
($Ac$=Pu, Am, Cf, Bk) with the above mentioned level of theory. We shall show
that relativistic effects together with electron correlation modify in a great
manner the contribution of those AOs, as well as the values of bonding
energies.

##  2 Methods and Computational details

All geometrical structures of the Pu-DOPO, Am-DOPO, Cf-DOPO, and Bk-DOPO
complexes reported in the present study were fully optimized with different
levels of theory: i) ZORA with hybrid DFT including dispersion corrections,
i.e., ZORA-B3LYP-D3 or DFT for short,18, 19, 20, 21, 22, 23 ii) ZORA-HF, which
means ZORA without electron correlation, calculated for only the Am-DOPO and
Bk-DOPO compounds, and iii) 4c DHF/LL Hamiltonians with basis sets that are
large enough for describing heavy atoms and those in the first coordination
sphere to the heavy atoms.29 We must mention here that we had performed 4c
calculations to obtain geometrical structures at DFT level of approach but in
this case we were unable to obtain converged results.

The first method essentially improves the usual non-relativistic hybrid DFT23,
24 methods and is useful for general chemistry applications. Long-range van
der Waals (vdW) dispersion correction were included.25 The Grimme’s–D3 terms
were used for all the light atoms, but due to availability, the dispersion
terms from the universal force field (UFF) were employed for actinides. The
Slater–Type Orbital (STO) basis sets were used instead of Gaussian-type
orbital (GTO) basis sets for the ZORA calculations. Valence triple-$\zeta$
with polarization functions (TZP) quality were used for Pu, Am, Cf, and Bk
atoms, and the valence double-$\zeta$ with polarization functions (DZP)
quality STO basis sets were used for the rest of the atoms (C, N, O, and H).30
In general, the STO basis sets give relatively consistent and rapidly
converging results.31 We have also made ZORA calculations without any DFT
functionals.

The second procedure is grounded on the 4c-DHF Hamiltonian, and then its Lévy-
Leblond model, to get appropriate non-relativistic (NR) results. In order to
compare results obtained with different relativistic methods, though not
including in any of them the spin-dependent contributions, we also performed
SF calculations to compare them with those of sr-ZORA. All studies were
performed by applying a local dense basis set scheme (LDBS),32 meaning basis
sets that are as large as possible for heavy atoms, and then, less complete
basis set for atoms that are not much involved in the main fragment of
interest (second coordination sphere and beyond). To accomplish highly
accurate four-component calculations we have designed an appropriate scheme of
basis sets. For atoms different to carbon and hydrogen the selected basis set
is dyall-cv3z,33 while carbon and hydrogen atoms are treated with 3-21G
Gaussian basis sets.34, 35 Geometries of the complexed systems and that of
each part (actinide atoms and DOPO molecules) were optimized using the DIRAC
code.36 Some LDBS schemes were used to describe more accurately the influence
of d-shell and f-shell AOs on the bonding of actinide atoms to the oxygen
atoms and nitrogen atoms. Only neutral systems (i.e., Am and Bk), were
considered at the DHF and LL levels of theory because it is not feasible to do
the geometry optimization of open-shell systems with analytical gradient in
the current implementation of the DIRAC code.

On the other hand, to calculate the bonding energies (BEs) at DHF and LL
levels, we first considered the closed–shell systems though ionized, meaning
Ac+3 and DOPO-1, and optimized their geometries. We did it by including at the
begining of the calculation a higher density of charges on the nitrogen atoms
(-1 each) and a lower charge density on the $Ac$ atoms (+3). Then, we
calculated the energy of the DOPO structure (which is an open–shell system) by
a single-point calculation. The energies of the actinides were also obtained
from a single–point calculation with the DIRAC code at the DHF and LL levels
of theory. On the other hand, the BEs of both actinides systems, Am–DOPO and
Bk–DOPO were calculated at the ZORA level of approach without any difficulty.
In all cases we have chosen to consider both, the close-shell configurations
and the open-shell configurations as shown in the supplementary information
(SI).

The whole set of model compounds considered in this work has been
experimentally synthesized and characterized, and are denoted as $Ac$–DOPO
(Figure 1). In all computations no constraints were imposed on the geometry.
The scalar relativistic ZORA-DFT calculations in addition to the so
relativistic ZORA-HF calculations were performed with the Amsterdam density
functional (ADF) code37, 38 while the 4c-DHF, 4c-LL and 4c-DHF-SF calculations
were performed with the Dirac code.36 The ADF-GUI38 program was used for the
visualization of the highest occupied molecular orbitals (HOMO) and lowest
unoccupied molecular orbitals (LUMO) from the ZORA-DFT calculations. The
ligand field diagram (LFD) analysis of the $Ac$-DOPO system ($Ac$ = Pu, Am, Cf
and Bk) for all actinides was performed using the ZORA approximation, while
DHF and LL levels of theory were used only for Am and Bk.

Figure 1: Optimized structures of $Ac$-DOPO compounds with a) $Ac$=Pu, b)
$Ac$=Am, c) $Ac$=Bk, and d) $Ac$=Cf.

##  3 Results and Discussion

We shall first show the influence of relativistic effects on the geometric
parameters of the actinide compounds, and then on BEs. Afterwards, we shall
show how important are the contributions of each of the valence-shell p-, d-
and f- AOs of actinides in the electronic framework of the frontier MOs, HOMO-
LUMO.

###  3.1 Geometrical structures

The optimized geometry of Pu–DOPO complex at the ZORA level is shown in Figure
1a. The root mean square deviation (RMSD) value of all the 1st coordination
bonds (Pu–N, Pu–O) w.r.t. to experiments is 0.104 Å, while the RMSD for the
2nd coordination sphere (C–C, C–N, and C–O) is 0.046 Å, as shown in Tables
S1-S3. In other words, the calculations at ZORA–DFT–D3 level of approach
estimate the experimental Pu–N and Pu–O bonds to 0.1 Åof error, which is the
order of magnitude for high quality x–ray structures.

The optimized geometry of the Am–DOPO complex is shown in Figure 1b. In this
case, we used four different levels of theory: ZORA–HF, ZORA–DFT, ZORA–DFT–D3,
4c–LL and 4c–DHF (see Table 1 and Table S4). The RMSD for the bond lengths of
Am–N and Am–O are 0.206, 0.024, 0.052, 0.023, 0.050 and 0.149, 0.120, 0.106,
0.103, 0.066 Å, respectively. This trend is in agreement with the increasing
accuracy of the level of theory, being 4c-DHF the highest and the smallest
values of RMSD is found for d(Am–O) though there are very small values of RMSD
for d(Am–N). In this last case only at the ZORA–HF level we obtain RMSD values
that are larger than 0.1 Å. Results of 4c-DHF-SF calculations gives the same
bond lengths as those of ZORA-DFT. Looking more closely at the bonds, the
calculations with ZORA–DFT–D3 are within a difference of 0.105 Åfor Am-O and
0.046 Åfor Am–N (Table S5). When they are calculated with the Lévy–Leblond
model, those differences are 0.102 Åand 0.029 Å, respectively. In the case of
DHF calculations such differences are within 0.066 Åand 0.054 Å, respectively.
This means that both, electron correlation and relativistic effects must be
included to best reproduce experimental measurements, though relativistic
effects showed by DHF calculations give the best reproduction of d(Am–N) and
d(Am–O). Remarkably, ZORA with dispersion corrections improves the RMSD only a
bit versus calculations without it.

Table 1: The equilibrium bond distances of Am-N and Am-O at different levels
of theory. All values are given in Å.

Bond length Difference w.r.t. experiments Bonds Expt.(a) ZORA-HF ZORA-DFT
ZORA-DFT-D3 4c-LL 4c-DHF ZORA-HF ZORA-DFT ZORA-DFT-D3 4c-LL 4c-DHF Am–N1 2.591
2.372 2.554 2.542 2.614 2.626 -0.219 -0.037 0.049 -0.023 -0.035 Am–N2 2.582
2.369 2.563 2.522 2.601 2.621 -0.213 -0.019 0.06 0.019 0.039 Am–N3 2.554 2.369
2.546 2.507 2.581 2.623 -0.185 -0.008 -0.047 0.027 0.069 RMSD 0.206 0.024
0.052 0.023 0.05 Am–O1 2.472 2.334 2.592 2.595 2.588 2.541 -0.138 0.12 0.123
0.116 0.069 Am–O2 2.465 2.321 2.614 2.573 2.571 2.536 -0.144 0.149 0.108 0.106
0.071 Am–O3 2.476 2.314 2.557 2.558 2.556 2.538 -0.162 0.081 0.082 0.08 0.062
Am–O4 2.462 2.333 2.579 2.579 2.58 2.545 -0.129 0.117 0.117 0.118 0.083 Am–O5
2.481 2.324 2.621 2.586 2.581 2.539 -0.157 0.14 0.105 0.1 0.058 Am–O6 2.493
2.332 2.592 2.586 2.583 2.544 -0.161 0.099 0.093 0.09 0.051 RMSD 0.149 0.12
0.106 0.103 0.066

(a) Taken from Ref. 3.

All this analysis suggest that dispersion corrections might not be necessary
to consider for this type of compounds, as they were not included in the 4c-LL
and 4c-DHF calculations either. Considering that electron correlation effects
can be obtained as the difference between ZORA–DFT and ZORA–HF calculations,
we found that these effects increase both distances: in average 0.166 Åfor
d(Am–N) and 0.267 Åfor d(Am–O). On the other hand, considering that
relativistic effects can be obtained as the difference between DHF and LL
calculations we found that they increase the values of d(Am–N) (0.024 Åin
average) but diminish that of d(Am–O) (0.037 Åin average) though in a much
smaller amount as compared with the influence of those distances by the
electron correlation effects.

The geometry of the Bk–DOPO complex was calculated with the three levels of
theory mentioned above. In Figure 1c we show its equilibrium geometry. The
bond distances Bk–N and Bk–O at equilibrium and for all methods used here are
collected in Table 2. The RMSD is 0.082, 0.114 and 0.063 for ZORA–DF–D3, 4c-LL
and 4c-DHF, respectively. The distances at 4c-DHF level of theory gives again
the smallest value of RMSD, which one might expect since it is the most
accurate. Similarly to what happens with $Ac$=Am, the bond lengths with 4c-DHF
have the same consistent sign, though calculations at ZORA and 4c-LL levels
underestimate the Bk–N bonds. Once both effects, electron correlation and
relativistic effects are introduced, the behavior of those bond distances
follows a similar trend as observed for Am. In this case, electron correlation
effects are quite similar to those of relativistic effects for d(Bk–N).

Table 2: The equilibrium bond distances Bk–N and Bk–O in the Bk–DOPO complex.
All values are given in Å.

Bond length Difference w.r.t. experiments Bonds Expt.(a) ZORA-HF ZORA-DFT-D3
LL DHF ZORA-HF ZORA-DFT-D3 LL DHF Bk–N1 2.521 2.338 2.494 2.547 2.585 0.183
-0.027 0.026 0.064 Bk–N2 2.571 2.406 2.441 2.432 2.584 0.165 -0.13 -0.139
-0.013 Bk–N3 2.544 2.354 2.492 2.486 2.587 0.178 -0.04 -0.046 0.055 RMSD 0.175
0.08 0.086 0.049 Bk–O1 2.459 2.346 2.59 2.579 2.478 0.114 0.13 0.119 0.018
Bk–O2 2.436 2.32 2.495 2.592 2.509 0.126 0.049 0.146 0.063 Bk–O3 2.454 2.336
2.505 2.569 2.492 0.117 0.052 0.116 0.039 Bk–O4 2.468 2.323 2.583 2.569 2.529
0.145 0.115 0.101 0.061 Bk–O5 2.447 2.329 2.532 2.555 2.534 0.118 0.085 0.108
0.087 Bk–O6 2.452 2.313 2.528 2.626 2.549 0.139 0.076 0.174 0.097 RMSD 0.127
0.09 0.13 0.066

(a) Taken from Ref. 3.

The last system studied was Cf–DOPO for which the oxidation state is +3. Its
equilibrium structure is shown in Figure 1d. The bond distances at
equilibrium, Cf–N and Cf–O, are collected in Table S6. There are no
experimental nor previous theoretical results for comparison of Cf complexes.
Our calculations at the ZORA–DFT–D3 level of theory shows that the average
Cf–N bond distance is slightly higher than the average Pu–N bond distance by
0.037 Å, and similarly the average Cf–O bond distance is slightly higher than
the average Pu–O bond distance by 0.083 Å. The average of bond distances Cf–N
and Cf–O are reported in Table S7, where, as happens with Am, the introduction
of dispersion correction to the ZORA calculations changes the bond length
slightly.

###  3.2 Bonding energies

For these calculations we shall consider the Am and Bk containing systems as
model compounds. The bonding energies (BE) are calculated, at DHF and LL
levels of approach, as: BE = Ac[DOPO]3 \- (Ac(III) + 3DOPO-1). Instead, at
ZORA level the BE is calculated as the difference of energies among the whole
molecule and its fragments. The main results are given in Table 3 (see also
Tables 10-17 of SI).

We have performed calculations at four-component levels in such a way that we
are able to analize the influence of relativistic effects (DHF - LL), SO
effects at DHF level of approach (DHF - DHF-SF) and electron correlation
(4c-B3LYP - DHF). We have also performed calculations at ZORA level of theory
and there we can analyze the influence of SO effects and the electron
correlation. Given that we include calculations that consider the actinide
atom with different multiplicities we are able to learn about what kind of
calculations give results that are closer to the four-component ones.

The behavior of relativistic effects is opposite for Am–(DOPO)3 and Bk–(DOPO)3
complexes. For the first actinide relativistic effects on BE are positive
(0.139 au) of which SO effects are the most important (0.105 au). On the
contrary, for berkelium relativistic effects on BE are of opposite sign
(-0.319 au) and in this case the SO effect is smaller (-0.083 au) when it is
compared with its equivalent in Am.

On the other hand, at sr-ZORA level of approach we found that BEs are much
dependent on the multiplicity of the actinide atom. One should compare results
of DHF-SF calculations with that of sr-ZORA/HF. For both actinides the more
approximated results of ZORA calculations are found for the quintuplet
multiplicity. Besides, given that the so-ZORA calculations are performed with
the closed-shell system its results must be comparable with that of DHF.

When electron correlation is considered at each of both levels of theory we
found different behaviors. At ZORA level the electron correlation for Am is
equal to: -2.128 au + 1.335 au = -0.793 au. For Bk it is -2362 + 1.850 =
-0.712. In the case of 4c calculations, for Am the electron correlation is:
-2.097 au + 1.475 au = -0.622 au which is close to that of ZORA. But in the
case of Bk the electron correlation is: -2.121 au + 1.784 au = -0.337 au which
is almost halve of that of the ZORA calculation.

Furthermore one may assume that when electron correlation effects are included
at 4c-B3LYP level of theory their results must be close to those of so-
ZORA/B3LYP calculations. In the case of Am-DOPO compound this is what we
found: -2.097 au vs -2.128 au (-5505.67 kJ/mol vs -5587.06 kJ/mol) but for Bk-
DOPO compound a little large difference appears: -2.121 au vs -2.362 au
(-5568.69 kJ/mol vs -6201.43 kJ/mol). We shall be aware that both
methodologies, ZORA/B3LYP and 4c/B3LYP have different grounds. The first one
use frozen core description of the electronic behavior but the other is an all
electron methodology.

We could shed some more light on the reason why those differences arises by
observing the energy-spectra of the frontier orbitals shown in Figure 2) for
the Am–(DOPO)3 complex. We show the energy of MOs that belongs to HOMO-LUMO
structures of sr-ZORA/B3LYP, 4c/B3LYP and 4c-DHF. There is a clear difference
among the energy-gap between HOMO-LUMO orbitals and also the influence of
f-type AOs in them.

Another finding is that by including electron correlation both molecular
systems become more stable.

Our previous analysis shows that higher levels of theory, such as DHF, might
be needed to obtain a better estimation not only of the geometry but also a
more accurate estimation of reaction energies. Additionally, the basis set
superposition error is less than 6% for the BE of both, Am–(DOPO)3 and
Bk–(DOPO)3 complexes.

Table 3: Bonding energies (BE) for Am and Bk, calculated at DHF, DHF-spin-free
(DHF-SF), Lévy-Leblond (LL) and sr (so)-ZORA levels of approach, expressed in
a.u. (between parentheses in kJ/mol).

Am(III).3DOPO-1 Am(III) DOPO-1 BEa $\Delta$Eb DHF -32694.990 -30488.849
-734.889 -1.475 (-3872.61) -0.0536 (140.81) DHF-SF -32619.097 -30412.630
-734.929 -1.680 (-4410.80) LL -29897.818 -27692.640 -734.521 -1.614 (-4237.56)
-0.0878 (-230.45) 4c-B3LYP -32722.070 -30503.042 -738.977 -2.097 (-5505.67)
so-ZORA/HF -32.023 -0.056 -10.211 -1.335 (-3505.04) sr-ZORA/HFc -31.835 1.229
-10.211 -2.432 (-6385.22) sr-ZORA/HFd -31.474 0.569 -10.211 -1.411 (-3704.58)
sr-ZORA/B3LYPc -20.814 1.457 -6.609 -2.445 (-6418.56) sr-ZORA/B3LYPd -20.756
1.296 -6.609 -2.226 (-5844.36) sr-ZORA/B3LYPe -20.593 1.096 -6.609 -1.863
(-4890.08) so-ZORA/B3LYP -20.997 1.001 -6.609 -2.128 (-5587.06)
Bk(III).3DOPO-1 Bk(III) DOPO-1 BEa $\Delta$Eb DHF -34397.417 -32190.967
-734.889 -1.784 (-4683.89) -0.060 (158.38) DHF-SF -34308.076 -32101.587
-734.929 -1.701 (-4465.98) LL -31301.618 -29096.589 -734.521 -1.465 (-3846.36)
-0.105 (-274.53) 4c-B3LYP -34425.524 -32206.471 -738.977 -2.121 (-5568.69) so-
ZORA/HF -32.172 0.310 -10.211 -1.850 (-4857.20) sr-ZORA/HFc -31.916 1.377
-10.211 -2.661 (-6986.46) sr-ZORA/HFd -31.707 0.780 -10.211 -1.855 (-4870.73)
sr-ZORA/B3LYPc -20.846 1.562 -6.609 -2.581 (-6776.42) sr-ZORA/B3LYPd -20.804
1.390 -6.609 -2.367 (-6215.83) sr-ZORA/B3LYPe -20.624 1.176 -6.609 -1.974
(-5182.32) so-ZORA/B3LYP -21.004 1.185 -6.609 -2.362 (-6201.43)

a At DHF and LL levels of approach: BE = Ac[DOPO]3 \- (Ac(III) + 3DOPO-1); at
ZORA level, BE is calculated as the difference of energies among the whole
molecule and its fragments.
b $\Delta E$ = E(DHF with optimized ZORA geometry) - E(DHF with optimized DHF
geometry) or E(LL with ZORA geometry) - E(LL with LL optimized geometry),
respectively.
c Singlet multiplicity.
d Quintuplet multiplicity.
e Septuplet multiplicity.

###  3.3 Influence of p-, d- and f- AOs on the HOMO-LUMO electronic structure

After getting the optimized structures of the Pu–(DOPO)3, Am–(DOPO)3,
Cf–(DOPO)3 and Bk–(DOPO)3 complexes, we have further carried out an
investigation of the contribution of the valence AOs of the heavy elements
(Pu, Am, Cf, and Bk) to the frontier MOs by analyzingthe pattern of the
HOMO–LUMO electronic structure and the LFD diagram. We got reliable results
for both, Am and Bk compounds at 4c level of approach (Table S8 and Table S9).
The other compounds are not closed–shell and so they cannot be studied with
the current version of the DIRAC code. On the other hand, the Lévy-Leblond
calculation of the structure of frontier orbitals of Bk–(DOPO)3 is not as
reliable for Am–(DOPO)3. In this last case the HOMO - LUMO orbitals of
Bk–(DOPO)3 have some inconsistencies in the LUMO structure, likely arising
from linear dependence.

Figure 2: Ligand Field Diagram (LFD) of the Am–DOPO system calculated at the
ZORA, DHF and LL levels of theory. Green lines represents the OM with f-type
primitives of the Am with significant Mulliken population.

To further understand the bonding structure of the electronic framework in the
Am(DOPO)3 complex, LFD computations were first performed at the ZORA level
(Figure 2). This calculation shows that the 5f–AOs of Am and DOPO ligands are
taking part to create MOs with bonding character. The HOMO and LUMO orbitals
in which those AOs are more involved are HOMO - 9, HOMO - 11, LUMO + 3 and
LUMO + 6 (Figures S6 and S7). The study of 5f– and 5p–AOs of Am contributions
to frontier MOs was also performed at both, the DHF and LL levels of theory to
corroborate the previous findings at the ZORA level. There are significant
contributions of Am 5f– and 6p–AOs to HOMO - 6 (fxyy = 0.271, fxxy = 0.271,
fxxx = 0.091) and LUMO + 3 (py = 0.633, px = 0.241, pz = 0.043) molecular
orbitals (Table S9). Thus, we found a similar scheme to that of ZORA–DFT–D3
though with little differences (see Figure 2). It is important to point out
that those relativistic effects are largely involved. At the LL level of
theory, the contribution of 5f–AO of Am to HOMO–6 is lost, and there is no
significant contribution of them to any HOMO - n (n = 1 - 12). On the other
hand, the pattern of the 5f–AOs contributions to LUMO + n orbitals (n = 9) are
just slightly modified, e.g., LUMO + 3 (py = 0.587, px = 0.234, pz = 0.043).
Thus, without the relativistic approximations (the LL reference method), there
are no contributions to the HOMO - n orbitals, which are present when ZORA or
DHF are used. This is a remarkable finding because it implies that
relativistic effects play a central role in the involvement of f–orbitals in
the frontier orbitals. Intuitively, this is an expected result due to the
well–known influence of direct/indirect relativistic effects on the size of
AOs, which means the contraction/expansion of the p–AOs/f–AOs. This is now
numerically observed on these actinide compounds. Furthermore, the pattern of
energies of the set of HOMO - n (n= 1 - 8) orbitals at both, four–component
and NR level of theory are quite similar.

A similar interaction symmetry can be seen in LUMO + 4 and LUMO + 5.

Figure 3: Ligand Field Diagram (LFD) of the Bk-DOPO system calculated at the
ZORA, DHF and LL levels of theory. The lines in green represents the OM with
f-type primitives of the Bk with significant Mulliken magnitudes.

The analysis of the contributions from Bk 5f–AOs to the frontier MOs was also
performed at the LL and DHF level of theory for comparison with ZORA
calculations. There are significant contributions of Bk 5f– and 6p–AOs to HOMO
- 8 (fxyy = 0.171, fxxy = 0.114, fxxx = 0.063) and LUMO + 3 (py = 0.578, px =
0.268, fxyy = 0.054) molecular orbitals at the DHF level (Table S10). The
non–relativistic calculations were performed with LL Hamiltonian, but in this
case some inconsistencies appears because the LUMO and LUMO + 1 orbital
energies are negative, and few linear dependencies arises in the wave function
calculation. In summary, as was also observed for the complexed Am compound
the symmetry, magnitude and types of Bk–AOs are affected by relativistic
effects, which in turn will affect how we interpret bonding and reaction on
these compounds. This calls for caution when interpreting standard DFT
calculations without scalar or 4–component relativistic effects for bonding or
reactions.

First, we focus on the Pu AOs involved in the bonding. The LFD of the
Pu–(DOPO)3 complex shows that the electronic 5f–orbitals of Plutonium are
interacting with the 2p–orbitals of both, the N and O atoms (Figure S3). To
provide further insight about the interaction between Pu and N or O atoms, we
have computed the HOMO - LUMO and HOMO - n and LUMO + n MOs (Figure S4 and
S5). These diagrams shows that the Pu AOs interact with the DOPO ligand AOs to
create molecular orbitals with bonding character. To give an example, three of
the highest bonding molecular orbitals, HOMO, HOMO - 6 and HOMO - 9 show that
5f–orbital electrons of Plutonium interact with the 2p–AOs of the N, O and
even C atoms. More specifically in the HOMO, the 5f–AO is delocalized with a
2p AO from N atom in the organic ligand (Figure S4). Similarly, we can observe
in LUMO, LUMO + 1, LUMO + 2 and LUMO + 6 that the Pu 5f–AOs are delocalized
with the 2p–AO from N and O atoms, but not C, in the ligand (Figure S5). This
is evidence that at the ZORA level of theory the 5f–AOs of the Pu plays an
important role in creating a bond with the organic ligand, i.e., 5f-AOs are
crucial to describe the bonding of the Pu(DOPO)3 complex.

The bonding interactions between Cf and DOPO (Figure S8) was studied by
analyzing the LFD of Cf–(DOPO)3. The LFD shows that near the frontier
orbitals, the 5f–AOs of Cf interact with the N, O and even with the C atoms in
the HOMO orbitals. However only C and N participate in the LUMO near the
frontier. The HOMO and LUMO orbitals were calculated at the ZORA level of
theory (Figures S9-10). In HOMO–7, HOMO–8 and HOMO–10, the overlap between
5f–AOs and AOs from the ligand can be observed, which have some bonding
character. A similar study was performed on Bk containing compounds. The LFD
shows that 5f–AOs of Bk are interacting with the carbon atom in the DOPO
ligand of the complex (Figure 3 and S11), which is also observed qualitatively
in the HOMO - LUMO representations at the ZORA level (Figures S12-13). The
overlap of the 5f–AOs can be seen in HOMO - 7, HOMO - 8 and HOMO - 9.

##  4 Conclusions

We computationally investigated the molecular structure at equilibrium of the
f-block element (Pu, Am, Bk, and Cf) compounds of the redox-active DOPO ligand
by using i) four–component with Dirac–Hartree–Fock (DHF), Lévy-Leblond (LL)
and spin-free (SF) Hamiltonians, together with ii) the scalar relativistic
(sr) and spin-orbit relativistic (so) zeroth–order regular approximation
(ZORA) with hybrid density functional theory (DFT).

We found that the RMSD values for bond distances are the smallest when the DHF
method is applied, followed by results of calculations with ZORA and LL levels
of theory. Taking both $Ac$–DOPO ($Ac$ = Am, Bk) compounds as model compounds,
we found that relativistic effects shorten the bond distances on average by
0.04Å for Am–O and 0.07Å for Bk–O. Furthermore, all six Am–O bond distances
and the three Am–N bond distances become closer to each other when
four–component relativistic effects are included, and they are also closer to
experimental distances as compared with sr-ZORA results. Total relativistic
effects shortened the $Ac$–O bond distances but expand the $Ac$–N bond
distances. When considering the scalar relativistic effects, we found that
they shorten $Ac$–N bonds lengths but the spin–dependent relativistic effects
expand them. The opposite behavior is found for the $Ac$–O bond lengths.
Furthermore, electron correlation stretches bond distances of $Ac$–N and
$Ac$–O though the effect on stretching the distance d($Ac$–O) is larger than
that effect on d($Ac$–N).

On the other hand, we found that the BE is more negative when it is calculated
at 4c-DHF level of theory than for so-ZORA/HF level of theory for Am-DOPO3
complex (-3872.61 kJ/mol vs -3505.04 kJ/mol, respectively). For the Bk-DOPO3
complex such a tendence is opposite (-4683.89 kJ/mol vs -4857.20 kJ/mol,
respectively). In both cases the inclusion of electron correlation at any
level of theory makes those systems more stables. This illustrates that
working with more accurate levels of theory, such as the one that consider
four-component Hamiltonians and all-electron calculations, might be needed to
obtain a better estimation of the geometrical structure and even a more
accurate estimation of reaction energies for $Ac$-DOPO3 complexes.

We have also analized the contributions of 5f– and 6p–AOs of central actinides
($Ac$ = Am, Bk) to frontier MOs at both levels of theory, ZORA and DHF. We
found that for $Ac$ = Am, there are significant contributions of Am 5f– and
6p–AOs to the HOMO - 6 (fxyy = 0.271, fxxy = 0.271, fxxx = 0.091) and LUMO + 3
(py = 0.633, px = 0.241, pz = 0.043). However, with the non-relativistic LL
Hamiltonian the Am 5f–AO contribution to HOMO - 6 is lost and also any other
HOMO - n (n = 1-12). This means that, without including relativistic effects
at ZORA or DHF level, no contributions of the Am 5f–AOs to the HOMO - n
orbitals are found. This shows that relativistic effects play a key role in
the population of f–orbitals in the frontier orbitals. A similar behavior
occurs for the Bk–containing compound, where also the symmetry and the
magnitude of Bk–AOs are affected by relativistic effects. In this case the
HOMO - 8 (fxyy = 0.171, fxxy = 0.114, fxxx = 0.063) and LUMO + 3 (py = 0.578,
px = 0.268, pz = 0.054) have Bk-AO contributions.

We can state now that calculations at both, ZORA and DHF (with a given local
dense basis set) levels of approach give good estimation of geometries, which
can be useful when the experimental synthesis is challenging or not cost-
effective, and also that the reaction of these complexes can only be
adequately reproduced if they are described within a relativistic framework,
given that the change depends on the frontier MOs. This is a call for caution
for using standard DFT calculations when the complex studied involves actinide
elements. One must include either, scalar or four–component relativistic
effects on the calculations of bondings and reactions. In addition, we want to
emphasize that one must include an adequate local dense basis set scheme for
making that results of calculations at DHF level be reliable.

As a continuation of our line of research and based on these results, we
started the study of some response properties like the NMR spectroscopic
parameters, which shall proportion new insights about the physics that is
behind these Actinide–containing systems.

##  5 Supporting Information

The equilibrium geometries with the x, y, z coordinators using the scalar
(ZORA) and the 4–component DHF and LL Hamiltonian are given in the Supporting
Information. G.A.A. acknowledge support from the Argentinian National Founding
for Science and Technique, FONCYT (Grant PICT2016-2936). S.P. thanks the
Science and Engineering Research Board, Department of Science and Technology
(SERB-DST), Govt. of India for providing the Ramanujan Faculty Fellowship
under the scheme no. SB/S2/RJN-067/2017 and the Early Career Research Award
(ECRA) under the project number ECR/2018/000255. The computing for this
project was performed partly on the High-Performance Computer cluster at the
Research Computing Center at Florida State University (FSU). The large portion
of this work was supported in part through computational resources and
services provided by the Institute for Cyber–Enabled Research at Michigan
State University (MSU).

##  6 Conflicts of interest

There are no conflicts to declare.

## References

  * Y. Ohtsuka 2019 Y. Ohtsuka, Y. T. Y. I. M. H. K. H. S. H., M. Aoyama _Sci. Rep._ 2019, _9_ , 1–8
  * Greenwood and Earnshaw 1997 Greenwood, N.; Earnshaw, A. _Chemistry of the Elements 2nd Edition_ ; Butterworth-Heinemann, 1997
  * Galley et al. 2019 Galley, S. S.; Pattenaude, S. A.; Gaggioli, C. A.; Qiao, Y.; Sperling, J. M.; Zeller, M.; Pakhira, S.; Mendoza-Cortes, J. L.; Schelter, E. J.; Albrecht-Schmitt, T. E., et al. Synthesis and characterization of tris-chelate complexes for understanding f-orbital bonding in later actinides. _Journal of the American Chemical Society_ 2019, _141_ , 2356–2366
  * Brenner et al. 2020 Brenner, N.; Sperling, J. M.; Poe, T. N.; Celis-Barros, C.; Brittain, K.; Villa, E. M.; Albrecht-Schmitt, T. E.; Polinski, M. J. Trivalent f-element squarates, squarate-oxalates, and cationic materials, and the determination of the nine-coordinate ionic radius of Cf (III). _Inorganic Chemistry_ 2020, _59_ , 9384–9395
  * Leszczyk et al. 2019 Leszczyk, A.; Tecmer, P.; Boguslawski, K. _Transition Metals in Coordination Environments_ ; Springer, 2019; pp 121–160
  * Pyykkö 2012 Pyykkö, P. Relativistic effects in chemistry: more common than you thought. _Annual review of physical chemistry_ 2012, _63_ , 45–64
  * Chang et al. 1986 Chang, C.; Pelissier, M.; Durand, P. Regular two-component Pauli-like effective Hamiltonians in Dirac theory. _Physica Scripta_ 1986, _34_ , 394
  * van Lenthe et al. 1994 van Lenthe, E.; Baerends, E.-J.; Snijders, J. G. Relativistic total energy using regular approximations. _The Journal of chemical physics_ 1994, _101_ , 9783–9792
  * Saue 2011 Saue, T. Relativistic Hamiltonians for chemistry: A primer. _ChemPhysChem_ 2011, _12_ , 3077–3094
  * A. Goel 2019 A. Goel, R. P. A. A. K., J. S. McCloy _J. Non–Cryst. Sol.: X_ 2019, _4_ , 100033
  * Vienna 2010 Vienna, J. D. Nuclear waste vitrification in the United States: recent developments and future options. _International Journal of Applied Glass Science_ 2010, _1_ , 309–321
  * Fisk et al. 1988 Fisk, Z.; Hess, D.; Pethick, C.; Pines, D.; Smith, J.; Thompson, J.; Willis, J. Heavy-electron metals: New highly correlated states of matter. _Science_ 1988, _239_ , 33–42
  * van Schilfgaarde et al. 1999 van Schilfgaarde, M.; Abrikosov, I.; Johansson, B. Origin of the Invar effect in iron–nickel alloys. _Nature_ 1999, _400_ , 46–49
  * Boring and Smith 2000 Boring, A. M.; Smith, J. L. Plutonium condensed-matter Physics. _Los Alamos Science_ 2000, _26_ , 90
  * Cao et al. 2021 Cao, C.; Vernon, R. E.; Schwarz, W.; Li, J. Understanding periodic and non-periodic chemistry in periodic tables. _Frontiers in Chemistry_ 2021, _8_ , 813
  * Cooper 2000 Cooper, N. G. Challenges in plutonium science. _Los Alamos Sci_ 2000, _26_ , 1–493
  * Galley et al. 2020 Galley, S. S.; Gaggioli, C. A.; Zeller, M.; Celis-Barros, C.; Albrecht-Schmitt, T. E.; Gagliardi, L.; Bart, S. C. Evidence of Alpha Radiolysis in the Formation of a Californium Nitrate Complex. _Chemistry–A European Journal_ 2020, _26_ , 8885–8888
  * Tobisu et al. 2012 Tobisu, M.; Kinuta, H.; Kita, Y.; Remond, E.; Chatani, N. Rhodium (I)-catalyzed borylation of nitriles through the cleavage of carbon–cyano bonds. _Journal of the American Chemical Society_ 2012, _134_ , 115–118
  * Kundu et al. 2009 Kundu, S.; Choliy, Y.; Zhuo, G.; Ahuja, R.; Emge, T. J.; Warmuth, R.; Brookhart, M.; Krogh-Jespersen, K.; Goldman, A. S. Rational design and synthesis of highly active pincer-iridium catalysts for alkane dehydrogenation. _Organometallics_ 2009, _28_ , 5432–5444
  * Silver et al. 2016 Silver, M. A.; Cary, S. K.; Johnson, J. A.; Baumbach, R. E.; Arico, A. A.; Luckey, M.; Urban, M.; Wang, J. C.; Polinski, M. J.; Chemey, A., et al. Characterization of berkelium (III) dipicolinate and borate compounds in solution and the solid state. _Science_ 2016, _353_ , aaf3762
  * Albrecht-Schmitt et al. 2020 Albrecht-Schmitt, T. E.; Hobart, D. E.; Páez-Hernández, D.; Celis-Barros, C. Theoretical examination of covalency in berkelium (IV) carbonate complexes. _International Journal of Quantum Chemistry_ 2020, _120_ , e26254
  * Goodwin et al. 2019 Goodwin, C. A.; Su, J.; Albrecht-Schmitt, T. E.; Blake, A. V.; Batista, E. R.; Daly, S. R.; Dehnen, S.; Evans, W. J.; Gaunt, A. J.; Kozimor, S. A., et al. Back Cover:[Am (C5Me4H) 3]: An Organometallic Americium Complex (Angew. Chem. Int. Ed. 34/2019). _Angewandte Chemie International Edition_ 2019, _58_ , 11924–11924
  * Becke 1993 Becke, A. D. Becke’s three parameter hybrid method using the LYP correlation functional. _J. Chem. Phys_ 1993, _98_ , 5648–5652
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Physical review B_ 1988, _37_ , 785
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _Journal of computational chemistry_ 2006, _27_ , 1787–1799
  * Lenthe et al. 1993 Lenthe, E. v.; Baerends, E.-J.; Snijders, J. G. Relativistic regular two-component Hamiltonians. _The Journal of chemical physics_ 1993, _99_ , 4597–4610
  * Van Lenthe et al. 1996 Van Lenthe, E.; Van Leeuwen, R.; Baerends, E.; Snijders, J. Relativistic regular two-component Hamiltonians. _International Journal of Quantum Chemistry_ 1996, _57_ , 281–293
  * Van Lenthe et al. 1996 Van Lenthe, E. v.; Snijders, J.; Baerends, E. The zero-order regular approximation for relativistic effects: The effect of spin–orbit coupling in closed shell molecules. _The Journal of chemical physics_ 1996, _105_ , 6505–6516
  * Cremer et al. 2014 Cremer, D.; Zou, W.; Filatov, M. Dirac-exact relativistic methods: the normalized elimination of the small component method. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2014, _4_ , 436–467
  * Van Lenthe and Baerends 2003 Van Lenthe, E.; Baerends, E. J. Optimized Slater-type basis sets for the elements 1–118. _Journal of computational chemistry_ 2003, _24_ , 1142–1156
  * Guell et al. 2008 Guell, M.; Luis, J. M.; Sola, M.; Swart, M. Importance of the basis set for the spin-state energetics of iron complexes. _The Journal of Physical Chemistry A_ 2008, _112_ , 6384–6391
  * Dyall 2007 Dyall, K. G. Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the actinides Ac–Lr. _Theoretical Chemistry Accounts_ 2007, _117_ , 491–500
  * Dyall 2012 Dyall, K. G. Core correlating basis functions for elements 31–118. _Theoretical Chemistry Accounts_ 2012, _131_ , 1–11
  * Binkley et al. 1980 Binkley, J. S.; Pople, J. A.; Hehre, W. J. Self-consistent molecular orbital methods. 21. Small split-valence basis sets for first-row elements. _Journal of the American Chemical Society_ 1980, _102_ , 939–947
  * Gordon et al. 1982 Gordon, M. S.; Binkley, J. S.; Pople, J. A.; Pietro, W. J.; Hehre, W. J. Self-consistent molecular-orbital methods. 22. Small split-valence basis sets for second-row elements. _Journal of the American Chemical Society_ 1982, _104_ , 2797–2803
  * 36 Gomes, A.; Saue, T.; Visscher, L.; Jensen, H.; Bast, R.; Aucar, I.; Bakken, V.; Dyall, K.; Dubillard, S.; Ekström, U., et al. DIRAC19. _Zenodo_
  * Te Velde et al. 2001 Te Velde, G. t.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. _Journal of Computational Chemistry_ 2001, _22_ , 931–967
  * te Velde et al. 2001 te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. _J. Comput. Chem._ 2001, _22_ , 931–967


# Deep Learning Based Model Identification System Exploits the Modular
Structure of a Bio-Inspired Posture Control Model for Humans and Humanoids

Vittorio Lippi1
1Neurological University Clinic, University of Freiburg, Freiburg im Breisgau,
Germany
<EMAIL_ADDRESS>
https://orcid.org/0000-0001-5520-8974

###### Abstract

This work presents a system identification procedure based on Convolutional
Neural Networks (CNN) for human posture control using the DEC (Disturbance
Estimation and Compensation) parametric model. The modular structure of the
proposed control model inspired the design of a modular identification
procedure, in the sense that the same neural network is used to identify the
parameters of the modules controlling different degrees of freedom. In this
way the presented examples of body sway induced by external stimuli provide
several training samples at once.

##  1 INTRODUCTION

The application of convolutional neural networks (CNN) in human movement
analysis has produced promising results in recent experiments [Batchuluun et
al., 2018, Karatzoglou et al., 2018, Abdu-Aguye and Gomaa., 2019a, Abdu-Aguye
and Gomaa., 2019b], and in general deep learning is starting to be applied to
system identification [De la Rosa and Yu, 2016, Andersson et al., 2019,
Miriyala and Mitra, 2020, De la Rosa et al., 2015, Ljung et al., 2020]. This
work aims to apply CNN to the identification of human posture control models,
that are used in several studies for the analysis of human data and the
control of humanoids, e.g. [van der Kooij et al., 2007, van der Kooij et al.,
2005, van Asseldonk et al., 2006, Goodworth and Peterka, 2018, Mergner, 2010,
Engelhart et al., 2014, Pasma et al., 2014, Jeka et al., 2010, Boonstra et
al., 2014]. Most of the studies on posture control use linear models like the
independent channel (IC) model [Peterka, 2002], and make the assumption of
linear and time-invariant behavior [Engelhart et al., 2016]. In this work a
nonlinear model will be considered. This makes the use of deep learning more
interesting, as the identification of nonlinear systems is more complex and in
general it is performed using computationally expensive iterative algorithms
as in [Assländer et al., 2015]. In a previous work, CNNs proved to be suitable
to identify the parameters of a human posture control model [Lippi. et al.,
2020]. Specifically, the CNN presented in [Lippi. et al., 2020] was applied to
a single inverted pendulum (SIP) model, but it proved to work also in
identifying the parameters of a double inverted pendulum system. In this work
a triple inverted pendulum (TIP) model of posture control will be used to show
how the identification process can be extended to an arbitrary number of
degrees of freedom. The model used here, the Disturbance Estimation and
Compensation (DEC) has a modular structure [Lippi and Mergner, 2017], i.e.
each degree of freedom is controlled by a module and all the modules share the
same structure. Exploiting this characteristic, one CNN is used to identify
the parameters of the three modules controlling the three degrees of freedom.
In this way each simulation provides three training samples at once.

##  2 METHODS

###  2.1 Posture Control Scenario

The scenario considered here models a human (or humanoid) balancing on a
tilting support surface as a TIP. The three degrees of freedom considered in
the sagittal plane are the ankles, the knees and the hips (Fig. 1A). The
support surface tilt $\alpha_{FS}$ represents the input of the system and it
is the same for all the simulations. The profile of the tilt of the support
surface is the pseudo-random ternary sequence, PRTS, shown in Fig. 1B. Such
stimulus is used in human experiments because it is not predictable for the
subject [Peterka, 2002]. It is composed by a sequence of velocity steps
suitable to excite the dynamics of the system over several frequencies. The
output of the system is the sway of the body segments: shank, thigh (leg) and
trunk addressed as $\alpha_{SS}$, $\alpha_{LS}$, and $\alpha_{TS}$
respectively.

###  2.2 Human and Humanoid Posture Control: The DEC Model

Figure 1: Balance scenario and posture control system. (A) The triple inverted
pendulum model used to simulate human posture control together with the
control modules. The $\alpha_{BS}$ represent the sway of the body COM around
the ankle joint. (B) The Pseudo-Random Ternary Signal, PRTS, the time profile
for support surface tilt. (C) Schema of a control module. The disturbances are
compensated feeding them in the form of an angle equivalent as input to the
servo controller (PID).

The DEC is a model of the human postural control mechanisms [Mergner, 2010,
Lippi and Mergner, 2017]. It has a modular structure that can be applied to
multiple DoF robots [Lippi and Mergner, 2017, Zebenay et al., 2015, Ott et
al., 2016, Lippi, 2018, Hettich et al., 2013, Hettich et al., 2015]. Here it
is implemented specifically on the TIP model, using three control modules. A
block diagram of the DEC control is shown in Fig.1. A complete description of
the model is beyond the scope of this paper and can be found in [Lippi and
Mergner, 2017]. In the following a general description of the model is
provided with particular emphasis on the components that have an impact on the
definition of the machine learning problem. A control module based on the DEC
concept is implemented as servo loop, here implemented as a PID controller
(the neural controller in Fig. 1C). The controlled variable can consist in the
COM sway of all the segments above the joint controlled by the module with
respect to the gravitational vertical passing through the controlled joint
(e.g. $\alpha_{BS}$ in Fig 1 for the ankle joint). The sensory channels shown
in Fig.1 (C) as Vestibular, Proprioceptive, and Force are used for disturbance
estimation. The disturbance estimates are fed into the servo so that the joint
torque on-line compensates for the disturbances. The lumped delay in Fig. 1
(C) represents all the delay effects that are distributed in humans, and
humanoids [Antritter et al., 2014, Hettich et al., 2014]. The model used in
this work considers gravity and support surface tilt as disturbances. The
gravity torque to be compensated is assumed to be proportional to the sway of
the COM of groups of segments above the controlled joints. Ideally the gain
required to compensate gravity is $mgh$, where $m$ and $h$ are the total mass
and the COM height of the controlled segments respectively, and $g$ is the
gravity acceleration. The integral component of the PID controller is not
applied to the compensation as in [Ott et al., 2016]. The support surface tilt
estimator includes a non-linearity, introduced to represent the behavior
observed in humans [Mergner et al., 2009, Mergner et al., 2003], defined as:

$\alpha_{FS}=\int_{0}^{t}f_{\theta}\left(\frac{d}{dt}\alpha_{BS}^{vest}+\frac{d}{dt}\alpha_{BF}^{prop}\right)$
(1)

where $\alpha_{BF}^{prop}$ is the ankle joint angle signal from
proprioception. $BF$ stands for Body-to-Foot. The function $f_{\theta}$ is a
dead-band threshold defined as

$f_{\theta}(\alpha)=\left\\{\begin{array}[]{llc}\alpha+\theta&if&\alpha<-\theta\\\
0&if&-\theta<\alpha<\theta\\\ \alpha-\theta&if&\alpha>\theta\\\
\end{array}\right.$ (2)

The estimated $\alpha_{FS}$ is then up-channeled through the control modules
and used to control body position. In this work the threshold is set as
$\theta=0.0003rad/s$. The disturbance compensation and the total torque
commanded by the servo controller for the ankle joint is:

$\tau(s)=-e^{-s\Delta}\left(K_{p}+sK_{d}+K_{i}\frac{1}{s}\right)\left(\epsilon\right)+(K_{p}+sK_{d})T_{g}$
(3)

where $K_{p}$ and $K_{d}$ is the derivative coefficient for the PD controller,
and $\Delta$ is the lump delay. Notice that the derivative component is also
acting on gravity $T_{g}$, representing a sort of anticipation of the
disturbance. There is also a passive torque acting on the joints represented
as a spring-damper model:

$\tau_{passive}=-(K_{p}^{pass}+sK_{d}^{pass})\left(\alpha_{BF}^{prop}\right)$
(4)

The coefficients $K_{p}^{pass}$ and $K_{d}^{pass}$ are kept fixed for all the
simulations. The standard parameters are shown in Table 1. In general the
gains are proportional to the $mgh$ for the groups of segments above the
controlled joint, this will be used in the next section to normalize the
parameters. The anthropometric of the TIP model is taken from previous works
[Lippi et al., 2019b, Lippi and Mergner, 2020].

###  2.3 The Training Set

#### Control parameters and target set.

The training set has been generated simulating the posture control scenario
with different parameters, these parameters represent the target for the
neural network and the sway of the segments represents the input. In order to
exploit the modular structure of the DEC control the CNN has been designed to
identify the parameters of a specific module. In this way each simulation
provides three training samples associated with the three control modules. As
explained in the previous section, the gains of the controllers are
proportional to the $mgh$ of the segments above the controlled joint, and the
lumped delay is larger for the ankle module and smaller for the hip module. In
order to use such parameters for the training of the same CNN their value is
expressed as deviation from the default value and normalized dividing it by
the default value, e.g. $\tilde{K_{p}}=(K_{p}^{sample}-K_{p})/K_{p}$ in Table
1. Besides the PID gains the target vector includes the lumped delay and a
variable $C$ that is set to $-1$ when the controlled variable is the joint
angle and to $+1$ when the controlled variable is the COM. This leads to a
sample with the following form:

$\tilde{T}=\left[\tilde{K_{p}}\ \ \tilde{K_{i}}\ \ \tilde{K_{d}}\ \
\tilde{\Delta}\ \ C\right]$ (5)

where the $\tilde{}$ indicates that the values are normalized.

The training samples are generated with random parameters from normal
distributions $X\sim\mathcal{N}(0,\,0.5)$ that is summed to the normalized
parameters. The variable $C$ is sampled with equal probability ($0.5$) between
the two cases. In order to avoid negative values the absolute value of the
obtained parameters is used ”‘warping”’ the normal distribution on positive
values. A set of parameters is used as a sample only if the behavior it
produces is stable: simulations with body sway larger than $6^{\circ}$ are not
considered realistic balancing scenarios and are discarded. Overall the
obtained data consist of $20730$ samples, $14000$ used as training set, $3000$
as validation set and $3730$ as test set. The three sets are normalized
subtracting the average divided by the variance of the training set. The
resulting distribution of the normalized $\tilde{K_{p}}$ is shown in Fig. 2
together with the distribution of maximum body sway amplitudes.

Figure 2: Output of the performed simulations. Above the distribution of the normalized $\tilde{K_{p}}$ parameter produced during the stable simulations. Below the distribution of the body sway amplitudes. The last column of the histogram on the right includes the simulations that were aborted because the body sway reached the threshold of $6^{\circ}$. Parameter | Symbol | Ankle | Knee | Hip | unit
---|---|---|---|---|---
Active proportional gain | $K_{p}$ | 465.98 | 245.25 | 73.575 | $\frac{N\cdot m}{rad}$
Active derivative gain | $K_{d}$ | 116.49 | 18.394 | 18.394 | $\frac{N\cdot m\cdot s}{rad}$
Active integrator gain | $K_{i}$ | 11.649 | 6.1312 | 1.8394 | $\frac{N\cdot m}{rad\cdot s}$
Passive stiffness | $K_{p_{pass}}$ | 232.5000 | 61.2500 | 36.5000 | $\frac{N\cdot m}{rad}$
Passive damping | $K_{d_{pass}}$ | 145.000 | 11.2500 | 11.2500 | $\frac{N\cdot m\cdot s}{rad}$
Foot rotation velocity threshold | $\theta_{vfs}$ | 0.03 |  |  | $rad/s$
Lumped delay | $\Delta$ | 0.10 | 0.07 | 0.1210 | $s$
Table 1: Default parameters

#### The CNN input.

The Input is a representation of the system output, i.e the sway of the
segments. The simulation was performed with a fixed integration step of $2ms$
and sampled at $50\ Hz$, producing $6051$ samples with a resolution of
$20\;ms$. In order to adapt the signal to the CNN the input was transformed
into an image. The data considered were the sway with respect to the vertical
of the two segments under and above the controlled joint, e.g. for the ankle
module they are $\alpha_{FS}$ and $\alpha_{LS}$. The Spectrogram of the two
signals is computed with short-time Fourier transform over windows of 250
samples overlapping for 135 computing 250 frequency points. This results in a
$51\times 250$ matrix of complex values. The first 51 columns of the matrix
are here considered, corresponding to a bandwidth of approximately $10\ Hz$.
Due to the low-pass dynamics of the body sway, the power spectrum of the
signal at higher frequency is very poor and hence not giving important
informations (i.e. almost black images independently of the parameters) . The
resulting $51\times 51$ matrix is used to define the input image as follows:
the module and the phase of the matrices associated to the two segments sway
are computed. The first channel of the image is the module of the matrix
describing the sway of the segment above, the second is the matrix describing
the segment below and the third channel is the difference between the phases.
The process is summarized in Fig 3. Also the input images are normalized
subtracting the average and dividing by the variance of the training set
through element-wise operation (i.e. pixel by pixel, channel by channel).

Figure 3: Design of the input features for the CNN. In A the spectrogram of a
signal is given, the dashed line shows where it is cut. In B it is shown how
the spectrograms associated with the sway of two body segments are combined to
form a three channel $51\times 51$ pixel image. The image C shows the
normalization of the samples using the mean and the variance of the training
set, computed element-wise (pixel by pixel).

###  2.4 CNN Architecture

The neural network architecture is schematized in Fig 4. The network is
implemented with Matlab™ Deep Learning Toolbox™. The network has been trained
using stochastic gradient descent with momentum as policy. The training was
set to a limit of 200 epochs. The loss function is the Mean Squared Error MSE
as expected with a regression task, although the target vector includes a
categorical feature, the controlled variable, that implies a classification
problem. The categorical variable produced by the CNN is interpreted
considering the sign (positive = COM sway, negative = joint angle). The
performance in regression and classification is discussed separately in the
Results section.

Figure 4: Neural Network architecture. The group of layers addressed with x 5
is repeated 5 times.

##  3 Results

###  3.1 Training and Test

The results of the training set are reported in Table 2. Notice that there is
a total RMSE on the training set, including the error on the categorical
variable, and a fit RMSE computed on the continuous variables. The average
absolute error and absolute variance is shown in Fig. 5. The error is larger
on $K_{i}$. The identification error in this context can be defined as the
norm of the difference between body sway obtained with the target parameters
($\bm{\alpha}_{BS}$) and the one associated with the identified parameters
($\bm{\tilde{\alpha}}_{BS}$) divided by the number of samples $N=6051$,

$E_{id}=\frac{\left\|\bm{\alpha}_{BS}-\bm{\tilde{\alpha}}_{BS}\right\|}{N}$
(6)

, where the bold text represents the fact that $\mathbf{\alpha}_{BS}$ is a
vector of samples. The MSE on the single sample is computed considering the
target error on the three control modules used in the simulation. The
prediction error $\epsilon_{p}$ on the 15 normalized parameters is computer
for each sample, leading to an MSE of

$MSE=\sqrt{\epsilon_{p}^{T}\epsilon_{p}}/15$ (7)

. The identification error plotted versus the MSE is shown in Fig. 6. On
average the identification error on the training set is $0.0024^{\circ}$

Figure 5: Average and variance of the error on the parameters.

Dataset total RMSE fit RMSE Accuracy Training 1.3664 1.3296 85.79% Validation
1.5186 1.4810 83.90% Test 1.4846 1.4486 84.72%

Table 2: Training Results Figure 6: Identification Error for a set of 60
training samples. The identification error is plotted versus the MSE error on
the parameters (target output). The arrow addresses a point that is out of the
scale of the axes (i.e MSE = 0.1, error=0.07). The correlation between the two
measures of errors is $-0.1265$, which means in general larger MSE on the
parameters are not associated to the larger identification errors. This may
suggest a certain degree of redundancy in the model, allowing the system to
reproduce a similar output with different sets of parameters.

###  3.2 Identification of Human Posture Control Parameters

The neural network is used to identify the control parameters from a human
experiment performed with the same scenario simulated in the training set. The
human experiment had the user stand with straight legs and hence the
simulation was applied just to the ankle and the hip module. The identified
parameters are reported in table 3 and the simulated sways of the segments are
shown in Fig. 7, the result shows a good similarity between the simulation and
the original data. For both the control modules, i.e. hip and ankle, the
network proposed the sway in space of the segments as controlled variable, in
according to previous studies on support surface tilt [Hettich et al., 2014].

Module $K_{p}$ $K_{i}$ $K_{d}$ $\Delta$ Var Ankle 421.8574 74.6664 12.0254
0.0685 Sway Angle Hip 74.5148 8.2559 1.8854 0.0219 Sway Angle

Table 3: Parameters identified for an example from a human experiment. Figure
7: Human data compared with simulated responses.

###  3.3 Identification with a Monolithic model

In order to evaluate the impact of the modular approach the experiment is
repeated with a network that considers the three DOF together. The target
vector incorporates the parameters for the 3 control modules, i.e. 15 values;
the input is encoded as a 5 channel picture with with spectrograms as
described in §2.3. Support surface tilt $\alpha_{FS}$ is not taken in account
as it is the same for all the samples. The channels represent the modulus of
the spectrogram for $\alpha_{LS}$, $\alpha_{MS}$, and $\alpha_{TS}$, and the
phase difference between $\alpha_{LS}$ and $\alpha_{MS}$ as well as between
$\alpha_{MS}$ and $\alpha_{TS}$. The dataset, rearranged in this way accounted
for $2264$ samples, i.e. one third of the original set. The dataset was split
into a training set and a validation set of $1000$ samples and a test set of
$264$ samples.The network had the same structure of the one used in the
modular case, but the input and the output layers were modified according to
the dimensionality of the new dataset. Such network has the disadvantage of
having less training samples and more parameters, but potentially the
advantage of integrating more information about the global structure of the
system (e.g. distinguishing explicitly between different DOFs). The results
are shown in Table 4.

Dataset total RMSE fit RMSE Accuracy Training 4.0731 3.2723 60.20% Validation
4.1597 3.3297 53.50% Test 4.0531 3.2586 51.51%

Table 4: Training Results with a monolithic model

The performance is overall worse than in the modular case, especially for the
classification accuracy, suggesting that a modular approach is advantageous
with this problem.

##  4 CONCLUSIONS AND FUTURE WORK

In this work a method for posture control parameter identification based on
CNN is presented. The system provides an efficient way to fit a model of the
non-linear bio-inspired control system DEC on experimental data. This
represents an advantage with respect to previous solutions relying on
iterative methods. the training set is produced with parameters from normal
distributions (although only the parameters producing a stable simulation are
selected). Future work will focus on the distribution of human data. An a
posteriori test can be performed comparing the distribution of the parameters
identified on the test set with the expected distribution on real data. This
can help the process of choosing between different possible network
hyperparameters sets as shown in [Sforza et al., 2011, Sforza and Lippi,
2013].

Some parameters were better identified than others as shown in Fig. 5. This
may be due to the kind of experimental set up or due to the choice of the
input features. For example the effect of $K_{i}$ is mainly visible just in
the low pass components of the spectrogram. Furthermore, simulations and robot
experiments are able to reproduce human behavior in the considered set up
(PRTS support surface tilt) without the integral component of the neural
controller [Mergner, 2010, Hettich et al., 2015, Hettich et al., 2014]. The
integral component of the controller is more important in analyzing transient
behavior, e.g. reaching a desired position, where it can be used to guarantee
zero tracking error [Ott et al., 2016]. A tracking task may be better to
identify $K_{i}$ more precisely.

The proposed CNN for posture control modeling can find application in bio-
inspired humanoid control, e.g. [Choi and Kim, 2007, Abedi and Shoushtari,
2012, Zebenay et al., 2015, Mergner and Lippi, 2018]. The CNN can also be
helpful in setting up wearable robots using the control parameters identified
on the user [Chugo et al., 2019, Mergner and Lippi, 2019]. The parameters can
be a tool to benchmark humanoids and wearable devices [Torricelli et al.,
2020], in particular, in the framework of the COMTEST project [Lippi et al.,
2019a, Lippi et al., 2020] that aims to make a posture control testbed
available for the humanoid robotics community, and to define performance
metrics.

## ACKNOWLEDGEMENTS

This work is supported by the project COMTEST, a sub-project of EUROBENCH
(European Robotic Framework for Bipedal Locomotion Benchmarking,
www.eurobench2020.eu) funded by H2020 Topic ICT 27-2017 under grant agreement
number 779963.

## REFERENCES

  * Abdu-Aguye and Gomaa., 2019a Abdu-Aguye, M. G. and Gomaa., W. (2019a).  On the feasibility of on-body roaming models in human activity recognition.  In Proceedings of the 16th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO,, pages 507–516. INSTICC, SciTePress.
  * Abdu-Aguye and Gomaa., 2019b Abdu-Aguye, M. G. and Gomaa., W. (2019b).  Versatl: Versatile transfer learning for imu-based activity recognition using convolutional neural networks.  In Proceedings of the 16th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO,, pages 507–516. INSTICC, SciTePress.
  * Abedi and Shoushtari, 2012 Abedi, P. and Shoushtari, A. L. (2012).  Modelling and simulation of human-like movements for humanoid robots.  In Proceedings of the 9th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO,, pages 342–346. INSTICC, SciTePress.
  * Andersson et al., 2019 Andersson, C., Ribeiro, A. H., Tiels, K., Wahlström, N., and Schön, T. B. (2019).  Deep convolutional networks in system identification.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 3670–3676. IEEE.
  * Antritter et al., 2014 Antritter, F., Scholz, F., Hettich, G., and Mergner, T. (2014).  Stability analysis of human stance control from the system theoretic point of view.  In Control Conference (ECC), 2014 European, pages 1849–1855. IEEE.
  * Assländer et al., 2015 Assländer, L., Hettich, G., and Mergner, T. (2015).  Visual contribution to human standing balance during support surface tilts.  Human movement science, 41:147–164.
  * Batchuluun et al., 2018 Batchuluun, G., Naqvi, R. A., Kim, W., and Park, K. R. (2018).  Body-movement-based human identification using convolutional neural network.  Expert Systems with Applications, 101:56–77.
  * Boonstra et al., 2014 Boonstra, T. A., van Vugt, J. P., van der Kooij, H., and Bloem, B. R. (2014).  Balance asymmetry in parkinson’s disease and its contribution to freezing of gait.  PLoS One, 9(7):e102493.
  * Choi and Kim, 2007 Choi, Y. and Kim, D. (2007).  On the balancing control of humanoid robot.  In Proceedings of the Fourth International Conference on Informatics in Control, Automation and Robotics - Volume 2: ICINCO,, pages 248–252. INSTICC, SciTePress.
  * Chugo et al., 2019 Chugo, D., Koyama, M., Yokota, M., Kawazoe, S., Muramatsu, S., Yokota, S., Hashimoto, H., Katayama, T., Mizuta, Y., and Koujina, A. (2019).  Sitting assistance that considers user posture tolerance.  In Proceedings of the 16th International Conference on Informatics in Control, Automation and Robotics - Volume 2: ICINCO,, pages 489–496. INSTICC, SciTePress.
  * De la Rosa and Yu, 2016 De la Rosa, E. and Yu, W. (2016).  Randomized algorithms for nonlinear system identification with deep learning modification.  Information Sciences, 364:197–212.
  * De la Rosa et al., 2015 De la Rosa, E., Yu, W., and Li, X. (2015).  Nonlinear system identification using deep learning and randomized algorithms.  In 2015 IEEE International Conference on Information and Automation, pages 274–279. IEEE.
  * Engelhart et al., 2016 Engelhart, D., Boonstra, T. A., Aarts, R. G., Schouten, A. C., and van der Kooij, H. (2016).  Comparison of closed-loop system identification techniques to quantify multi-joint human balance control.  Annual reviews in control, 41:58–70.
  * Engelhart et al., 2014 Engelhart, D., Pasma, J. H., Schouten, A. C., Meskers, C. G., Maier, A. B., Mergner, T., and van der Kooij, H. (2014).  Impaired standing balance in elderly: a new engineering method helps to unravel causes and effects.  Journal of the American Medical Directors Association, 15(3):227–e1.
  * Goodworth and Peterka, 2018 Goodworth, A. D. and Peterka, R. J. (2018).  Identifying mechanisms of stance control: a single stimulus multiple output model-fit approach.  Journal of neuroscience methods, 296:44–56.
  * Hettich et al., 2014 Hettich, G., Assländer, L., Gollhofer, A., and Mergner, T. (2014).  Human hip—ankle coordination emerging from multisensory feedback control.  Human Movement Science, 37:123–146.
  * Hettich et al., 2013 Hettich, G., Lippi, V., and Mergner, T. (2013).  Human-like sensor fusion mechanisms in a postural control robot.  In Londral, A. E., Encarnacao, P., and Pons, J. L., editors, Proceedings of the International Congress on Neurotechnology, Electronics and Informatics. Vilamoura, Portugal, pages 152–160.
  * Hettich et al., 2015 Hettich, G., Lippi, V., and Mergner, T. (2015).  Human-like sensor fusion implemented in the posture control of a bipedal robot.  In Neurotechnology, Electronics, and Informatics, pages 29–45. Springer.
  * Jeka et al., 2010 Jeka, J. J., Allison, L. K., and Kiemel, T. (2010).  The dynamics of visual reweighting in healthy and fall-prone older adults.  Journal of motor behavior, 42(4):197–208.
  * Karatzoglou et al., 2018 Karatzoglou, A., Schnell, N., and Beigl, M. (2018).  A convolutional neural network approach for modeling semantic trajectories and predicting future locations.  In International Conference on Artificial Neural Networks, pages 61–72. Springer.
  * Lippi, 2018 Lippi, V. (2018).  Prediction in the context of a human-inspired posture control model.  Robotics and Autonomous Systems.
  * Lippi and Mergner, 2017 Lippi, V. and Mergner, T. (2017).  Human-derived disturbance estimation and compensation (dec) method lends itself to a modular sensorimotor control in a humanoid robot.  Frontiers in neurorobotics, 11:49.
  * Lippi and Mergner, 2020 Lippi, V. and Mergner, T. (2020).  A challenge: Support of standing balance in assistive robotic devices.  Applied Sciences, 10(15):5240.
  * Lippi. et al., 2020 Lippi., V., Mergner., T., and Maurer., C. (2020).  Deep learning for posture control nonlinear model system and noise identification.  In Proceedings of the 17th International Conference on Informatics in Control, Automation and Robotics - Volume 1: ICINCO,, pages 607–614. INSTICC, SciTePress.
  * Lippi et al., 2020 Lippi, V., Mergner, T., Maurer, C., and Seel, T. (2020).  Performance indicators of humanoid posture control and balance inspired by human experiments.  In The International Symposium on Wearable Robotics (WeRob2020) and WearRAcon Europe.
  * Lippi et al., 2019a Lippi, V., Mergner, T., Seel, T., and Maurer, C. (2019a).  COMTEST project: A complete modular test stand for human and humanoid posture control and balance.  In 2019 IEEE-RAS 19th International Conference on Humanoid Robots (Humanoids) Toronto, Canada. October 15-17.
  * Lippi et al., 2019b Lippi, V., Molinari, F., and Seel, T. (2019b).  Distributed bio-inspired humanoid posture control.  In 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 5360–5365. IEEE.
  * Ljung et al., 2020 Ljung, L., Andersson, C., Tiels, K., and Schön, T. B. (2020).  Deep learning and system identification.  In 21st IFAC World Congress, page 8.
  * Mergner, 2010 Mergner, T. (2010).  A neurological view on reactive human stance control.  Annual Reviews in Control, 34(2):77–198.
  * Mergner and Lippi, 2018 Mergner, T. and Lippi, V. (2018).  Posture Control–Human-Inspired Approaches for Humanoid Robot Benchmarking: Conceptualizing Tests, Protocols and Analyses.  Frontiers in Neurorobotics, 12:21.
  * Mergner and Lippi, 2019 Mergner, T. and Lippi, V. (2019).  Integrating posture control in assistive robotic devices to support standing balance.  In Carrozza M., Micera S., P. J., editor, Wearable Robotics: Challenges and Trends. WeRob 2018., volume 22 of Biosystems & Biorobotics. Springer, Cham.
  * Mergner et al., 2003 Mergner, T., Maurer, C., and Peterka, R. J. (2003).  A multisensory posture control model of human upright stance.  Progress in Brain Research, 142:189–201.
  * Mergner et al., 2009 Mergner, T., Schweigart, G., and Fennell, L. (2009).  Vestibular humanoid postural control.  Journal of Physiology - Paris, 103:178–194.
  * Miriyala and Mitra, 2020 Miriyala, S. S. and Mitra, K. (2020).  Deep learning based system identification of industrial integrated grinding circuits.  Powder Technology, 360:921–936.
  * Ott et al., 2016 Ott, C., Henze, B., Hettich, G., Seyde, T. N., Roa, M. A., Lippi, V., and Mergner, T. (2016).  Good posture, good balance: Comparison of bioinspired and model-based approaches for posture control of humanoid robotsg.  IEEE Robotics & Automation Magazine, 23(1):22–33.
  * Pasma et al., 2014 Pasma, J., Engelhart, D., Schouten, A., Van der Kooij, H., Maier, A., and Meskers, C. (2014).  Impaired standing balance: the clinical need for closing the loop.  Neuroscience, 267:157–165.
  * Peterka, 2002 Peterka, R. (2002).  Sensorimotor integration in human postural control.  Journal of neurophysiology, 88(3):1097–1118.
  * Sforza and Lippi, 2013 Sforza, F. and Lippi, V. (2013).  Support vector machine classification on a biased training set: Multi-jet background rejection at hadron colliders.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 722:11–19.
  * Sforza et al., 2011 Sforza, F., Lippi, V., and Chiarelli, G. (2011).  Rejection of multi-jet background in pp→ e$\nu$\+ jj channel through a svm classifier.  In Journal of Physics: Conference Series, volume 331, pages 32–45. IOP Publishing.
  * Torricelli et al., 2020 Torricelli, D., Mizanoor, R. S., Lippi, V., Weckx, M., Mathijssen, G., Vanderborght, B., Mergner, T., Lefeber, D., and Pons, J. L. (2020).  Benchmarking human likeness of bipedal robot locomotion: State of the art and future trends.  In Metrics of Sensory Motor Coordination and Integration in Robots and Animals, pages 147–166. Springer.
  * van Asseldonk et al., 2006 van Asseldonk, E. H., Buurke, J. H., Bloem, B. R., Renzenbrink, G. J., Nene, A. V., van der Helm, F. C., and van der Kooij, H. (2006).  Disentangling the contribution of the paretic and non-paretic ankle to balance control in stroke patients.  Experimental neurology, 201(2):441–451.
  * van der Kooij et al., 2005 van der Kooij, H., van Asseldonk, E., and van der Helm, F. C. (2005).  Comparison of different methods to identify and quantify balance control.  Journal of neuroscience methods, 145(1-2):175–203.
  * van der Kooij et al., 2007 van der Kooij, H., van Asseldonk, E. H. F., Geelen, J., van Vugt, J. P. P., and Bloem, B. R. (2007).  Detecting asymmetries in balance control with system identification: first experimental results from parkinson patients.  Journal of Neural Transmission, 114(10):1333.
  * Zebenay et al., 2015 Zebenay, M., Lippi, V., and Mergener, T. (2015).  Human-like humanoid robot posture control.  In 2015 12th International Conference on Informatics in Control, Automation and Robotics (ICINCO), volume 2, pages 304–309. INSTICC, SciTePress.