# HAIR: Head-mounted AR Intention Recognition

David Puljiz  Intelligent Process Automation and Robotics Lab (IPR), Institute
for Anthropomatics and Robotics, Karlsruhe Institute of
TechnologyKarlsruheGermany<EMAIL_ADDRESS>,  Bowen Zhou  Intelligent
Process Automation and Robotics Lab (IPR), Institute for Anthropomatics and
Robotics, Karlsruhe Institute of TechnologyKarlsruheGermany ,  Ke Ma
Intelligent Process Automation and Robotics Lab (IPR), Institute for
Anthropomatics and Robotics, Karlsruhe Institute of TechnologyKarlsruheGermany
and  Björn Hein  Intelligent Process Automation and Robotics Lab (IPR),
Institute for Anthropomatics and Robotics, Karlsruhe Institute of
TechnologyKarlsruheGermany

(2021)

###### Abstract.

Human teams exhibit both implicit and explicit intention sharing. To further
development of human-robot collaboration, intention recognition is crucial on
both sides. Present approaches rely on a vast sensor suite on and around the
robot to achieve intention recognition. This relegates intuitive human-robot
collaboration purely to such bulky systems, which are inadequate for large-
scale, real-world scenarios due to their complexity and cost. In this paper we
propose an intention recognition system that is based purely on a portable
head-mounted display. In addition robot intention visualisation is also
supported. We present experiments to show the quality of our human goal
estimation component and some basic interactions with an industrial robot.
HAIR should raise the quality of interaction between robots and humans,
instead of such interactions raising the hair on the necks of the human
coworkers.

Human Intention Estimation, Augmented Reality, Human-robot Collaboration, Head
Mounted Displays

††copyright: acmcopyright††journalyear: 2021††conference: Virtual, Augmented
and Mixed reality human-robot interaction Workshop; Virtual, Augmented and
Mixed reality human-robot interaction Workshop††booktitle: ; Figure 1. The
view from the HoloLens on the left, showing the user looking at the goal
”Sphere” and approaching it with the hand. On the right the current and
previous probabilities. The user rotated from the first goal, the sphere,
through the other two goals and back towards the sphere. To note is how the
output proclaims the user irrational when they are facing away from all the
defined goals

##  1\. Introduction

(a)

(b)

Figure 2. YOLOv4 classification of HoloLens camera data. Bounding box
classification approaches together with the known HoloLens egomotion and depth
data can be used to define the interaction objects/goals as well as to
outsource part of the environmental sensing from the robot to the human.

Communicating intentions between members of a team is paramount for successful
cooperation and task completion. Previous work in the field of Augmented
reality (AR) human-robot interaction (HRI) focused on either improving robot
programming (Quintero et al., 2018) or visualising robot motions (Walker et
al., 2018). Although quite important for collaboration, such systems still
lack the estimation of the human intention from the robot’s side. Several such
systems have been proposed, such as (Bascetta et al., 2011) where the human is
tracked and their goal estimated inside a robot cell. Such systems however
require a big overhead in complexity and cost of the robot cells. With the
advent of Head-mounted Displays (HMDs) the possibility arises of a fully
portable, completely worn system possessing both the robot intention
visualisation and human intention estimation.

A similar system based on a HMD and intended for human-robot collaborative
task planning was presented by Chakraborti et al. (Chakraborti et al., 2018).
In the presented system, however, the human coworker had to explicitly select
and reserve objects it wished to interact with, slowing down task execution
and increasing physical and mental workloads on the human worker.

Here we propose a system that implicitly evaluates the intentions of the
human, thus minimising the increase in workload. The proposed system is based
on the Microsoft’s HoloLens HMD and is aimed at a collaborative scenario
between a single human and an industrial manipulator. The system is robot
agnostic and completely portable requiring a very short set-up at the
beginning of interaction. This guarantees that a single human worker can
interface with multiple robots one after another, without the need for
specialised robot cells or sensors around any of those robots.

The system takes as input the pose of the HMD in the world coordinate system,
the position of the hand joints on the world coordinate system as well as a
set of possible spatial goals, which can be added and removed during the
interaction itself. The output is a set of probabilities of the goal the human
wants to approach as well as the action they wish to perform.

This paper will present our current work and tests aimed mostly at having a
robust goal estimation. To the best of the authors knowledge such an intention
estimation algorithm using a completely worn system has not yet been
developed.

##  2\. Methodology

###  2.1. Referencing

First and foremost a common coordinate system must be establish between the
HMD and the robotic manipulator. Referencing can be done in a variety of ways,
perhaps the most popular is the use of QR codes or other preset visual markers
(Krupke et al., 2018). Although these offer continuous instead of one-shot
referencing, as well as very good precision, they require a setup step that we
would like to avoid. Manually selecting the robot base such as presented in
(Quintero et al., 2018) is more flexible yet also more imprecise. We have
proposed several referencing methods in (Puljiz et al., 2019), with the semi-
automatic one, consisting of a rough user guess followed by a refinement step,
offering the best balance between accuracy and computational time. The
refinement step consists of filtering a point cloud captured with the HMD and
using a registration algorithm to fit the model of the manipulator into the
filtered point cloud, using the user guess as the start point of the
registration algorithm. The user guess prevents the common problem of
registration algorithms being stuck in the local minima, and we found that
even a basic ICP algorithm performs a good job of refining the guess of the
user. Another approach is a fully automatic one without a user guess. Such a
referencing algorithm, similarly to our automatic method proposed in (Puljiz
et al., 2019), was proposed by Ostantin et al. (Ostanin et al., 2020). It
clusters the point cloud captured by the HMD using the DBSCAN clustering
algorithm and then performs model matching between the clusters and a model of
the robot.

###  2.2. Defining Spatial Goals

Secondly the possible spatial goals of the human and the robot need to be
defined. In case the robot does not possess the full map of its surroundings,
the HMD can also provide that as we demonstrated in (Puljiz et al., 2020).
This can also include possible goals and regions of interest such as the table
or the conveyor belt. If goals are specific objects, here too the HMD can
provide types and positions of those objects. One such possibility is through
the use of bounding-box classifiers such as YOLOv4 (Bochkovskiy et al., 2020).
In Fig. 2, one can see the result of running YOLOv4 on the HoloLens camera
data. Having the egomotion data of the hololens, as well as data from the
HMD’s depth sensor, allows a full spatial definition of the objects and
therefore possible goals.

The user should also be able to add and remove goals manually during the
interaction step. Therefore the goal estimation algorithm was selected to
allow such a modality.

Figure 3. The HMM states of the human goal intention estimation system as
presented in (Petković et al., 2019). It consists of g goal states, a state of
unknown intention $\textbf{G}_{?}$ and a state of the human acting
irrationally $\textbf{G}_{x}$.

###  2.3. Goal Estimation

Finally, having a common coordinate system, mapped working environment and
possible goals, one can infer the goals using a human intention recognition
algorithm. We base our HIR algorithm on previous work by Petkovic et al.
(Petković et al., 2019) where a hidden Markov model framework was used to
estimate the goal of the human in an automated warehouse. The approach is
quite general and with minimal modification can be adapted to be used in our
use case. In this section we will present a brief overview of the calculation,
for more details please refer to the original paper.

Instead of the position of the human coworker as in the original paper, we
consider the position of the hand in relation to goal objects. To simplify the
calculations, we assume there is an almost straight line between the hand
position and each goal object. By doing that we can forego the complex path
planning step to determine the modulated distance and instead use the
euclidean distance to calculate the vector d that represents the distance of
the hand to each goal. As in (Petković et al., 2019), we define additional 32
points $p_{i}$ on a circle around the previous hand position $l^{\prime}$ and
a radius $r$ equal to the distance between the current $l$ and the previous
$l^{\prime}$ hand positions. We calculate the vector d for each point $p_{i}$
and append them to the modulated distance matrix D.

Additionally we consider the gaze validation s of the HMD. The motivation
being that the user is more likely to look approximately towards the goal of
the hand motion than towards other goals. The gaze validation is calculated
as:

(1)
$\textbf{s}_{i}=\begin{dcases}\textbf{g}\cdot\frac{\textbf{o}_{i}-\textbf{h}}{||\textbf{o}_{i}-\textbf{h}||},&\textbf{g}\cdot\frac{\textbf{o}_{i}-\textbf{h}}{||\textbf{o}_{i}-\textbf{h}||}\geq
0\\\ 0,&\text{otherwise.}\end{dcases}$

Where g is the HMD orientation in the world coordinate system,
$\textbf{o}_{i}$ is the position of object i and h is the position of the HMD.
We expand the motion validation vector v as follows:

(2) $\textbf{v}=\frac{\underset{1\leq i\leq
n}{\max}{\textbf{D}_{ij}}-\textbf{d}}{\underset{1\leq i\leq
n}{\max}{\textbf{D}_{ij}}-\underset{1\leq i\leq
n}{\min}{\textbf{D}_{ij}}}\cdot\textbf{s}$

The rest follows exactly the algorithm described in (Petković et al., 2019).
We use the same transition matrix with g goals $\textbf{T}^{g+2\times g+2}$
defined as:

(3) $\textbf{T}=\begin{bmatrix}1-\alpha&0&\dots&\alpha&0\\\
0&1-\alpha&\dots&\alpha&0\\\ \vdots&&\ddots&&\vdots\\\
\beta&\beta&\dots&1-g\beta-\gamma&\gamma\\\ 0&0&\dots&\delta&1-\delta\\\
\end{bmatrix},$

This transition matrix corresponds to the hidden Markov model (HMM)
architecture visible in Fig. 3. The parameter $\alpha$ captures the worker
tendency to change their mind, while the parameter couple $\beta$ and $\gamma$
set the threshold for estimating intention for each goal location. Increasing
$\beta$ leads to quicker inference of worker’s intentions and increasing
$\gamma$ speeds up the decision making process. Parameter $\delta$ captures
model’s reluctance to return to estimating the other goal probabilities once
it estimated that the worker is irrational. We performed several tests to
determine the optimal values of these parameters which will be described in
the ”Experiments” section.

The worker intention is estimated using the Viterbi algorithm (Forney, 1973),
which takes as inputs the hidden states set
$S=\\{G_{1},...G_{g},G_{?},G_{x}\\}$, hidden state transition matrix T,
initial state $\Pi$, sequence of observations O, and the emission matrix B.

The emission matrix B is calculated using the motion validation vector v.
Since the observation is the validation vector v with continuous element
values, the input to the Viterbi algorithm was modified by introducing an
expandable emission matrix $\textbf{B}^{k\times g}$, where $k$ is the recorded
number of observations, are functions of the observation value. Once a new
validation vector $v$ is calculated, the emission matrix is expanded with the
row $\textbf{B}^{\prime}$, where the element $\textbf{B}^{\prime}_{i}$ stores
the probability of observing v from hidden state $G_{i}$. The average of the
last $m$ vectors v is also calculated and the maximum average value $\phi$ is
selected. It is used as an indicator if the worker is behaving irrationally,
i.e., is not moving towards any predefined goal. The value of the
hyperparameter $m$ indicates how much evidence is to be collect before the
worker is declared irrational. If the worker has been moving towards at least
one goal in the last $m$ iterations ($\phi>0.5$), $\textbf{B}^{\prime}$ is
calculated as:

(4)
$B^{\prime}=\zeta\cdot\begin{bmatrix}\tanh(\textbf{v})&\tanh(1-\Delta)&0\end{bmatrix},$

and otherwise as:

(5) $B^{\prime}=\zeta\cdot\begin{bmatrix}\boldsymbol{0}_{1\times
g}&\tanh(0.1)&\tanh(1-\phi)\end{bmatrix},$

where $\zeta$ is a normalising constant and $\Delta$ is calculated as the
difference of the largest and second largest element of v.

Finally, the initial probabilities of worker’s intentions are set as:

(6) $\Pi=\begin{bmatrix}0&\dots&0&1&0\end{bmatrix},$

indicating that the initial state is $G_{?}$ and the model does not know which
goal the worker desires the most. The Viterbi algorithm outputs the most
probable hidden state sequence and the probabilities $P(G_{i})$ of each hidden
state in each step. These probabilities are the worker’s intention estimates.

Goals can be added and removed during runtime as well making such a intention
estimation framework quite flexible.

Figure 4. The result of hand detection presented in (Mueller et al., 2018) on
HoloLens RGB camera data. One can see robust performance even during object
handling.

###  2.4. Action Estimation

Though estimating the goal of the human motion is extremely important for
replanning robot motions to keep the interaction both safe and efficient,
estimating actions the human wishes to perform could also bring additional
information and flexibility to intention estimation systems.

Although the first generation of the HoloLens possesses inbuilt hand-tracking
capabilities, these are quite limited and only four gestures can be tracked
and classified. For a more robust hand following and classification we
expanded the hand tracking by using the work presented in (Mueller et al.,
2018) on the HoloLens’ RGB camera data. The algorithm tracks 21 hand joints
and works with occlusions, surface contacts and object handling.

The detected hand joints are to classify actions - intention to grasp an
object, grasped an object, pointing and stop. More actions can be classified
in the future. The stop and pointing gestures are used as simple cues to
control the robot. In addition, common gestures of fear or distress shall be
classified as stop gestures, allowing the system indirect reaction to stress.

We presently only detect and use the right hand, however with a slight
overhead the algorithm can detect both provided there is no significant
overlap between them.

Figure 5. An example of visualisation of the intentions of the robot to the
human coworker - virtual execution with a holographic robot and plan
visualisation.

###  2.5. Robot Intention visualisation

The benefits of HMDs extend also to visualising robot intention. Instead of
adapting robot motions to make them more legible (Dragan et al., 2015), one
can use holograms to signal the desired goal. In (Williams et al., 2019) it
was shown that holographic information is adequate to show the goal of the
robot, and even solve ambiguities if intention is expressed via synthesised
voice. General motion intent can also be effectively visualised using
holographic cues (Walker et al., 2018). In our work we chose to indicate the
goal via a hologram containing a 3D sound source (spatial sound), as well as
virtual execution - having a hologram of the robot execute the motion before
the real robot performs it, such as shown in Fig. 5.

Figure 6. The effect of increasing the parameter $\alpha$. The parameter
captures the worker tendency to change their mind. A low $\alpha$ will make
the algorithm ”too sure” about the intention, while a too high $\alpha$
produces a chaotic and unusable output.

##  3\. Experiments

The experiments were aimed at testing the performance of the goal intention
estimation. We used three goals in a circular pattern, from left to right - a
green cylinder, a red cube and a blue sphere, as shown in Fig. 7.

The first set of experiments was aimed to find the optimal set of parameters
$\alpha$, $\beta$, $\gamma$ and $\delta$ for our use case. Here we looked at
the goal states and the transitions between them. The path was a simple left
to right one, first going towards the green cylinder, then the red cube and
finally the blue sphere. Figure Fig. 6 shows the behaviour of the parameter
$\alpha$. A low value makes the algorithm too certain and almost does not
spend time in the unknown state, while a high value makes the estimated goals
jump too much. The optimal value of the parameter $\alpha$ was therefore set
to $\alpha=0.3$.

With the parameter $\alpha$ set, we tested the behaviour of changing the
parameter $\beta$. A low $\beta$maintains the unknown intention state too
long, while a high $\beta$ lowers the general certainty but eliminates the
insecurities between state transitions which is an unwanted behaviour. The
value of $\beta$ was set to $\beta=0.05$. In Fig. 8 one can see the behaviour
of changing the parameter $\beta$.

The parameters $\gamma$ and $\delta$ did not significantly influence the
outputs and were kept at the same value as in (Petković et al., 2019), namely
$\gamma=0.05$ and $\delta=0.1$.

With the parameters set we examined how the algorithm behaves with different
sequences of goals. The results are visible in Fig. 9.

The first test on the left is is a simple sequence of goals from left to
right. One can notice the algorithm goes into the state of unknown intention
during the transitions. The slower the transition the longer the unknown
state. One can also see the small drop near the end when the hand tracking was
lost and the gaze was not directly towards the sphere.

The second experiment starts with the middle goal, the cube, then goes left to
the cylinder, back to the cube and then towards the sphere. One can see that
the transition from cylinder to cube lasts slightly longer than from cube to
cylinder. This is due to the fact that the algorithm is reluctant to estimate
an already visited or skipped goal. One can also see the long transition
between the cube and the sphere, as the algorithm prefers the goal that has
already been visited two times. This shows that the estimation follows our
intuition.

Finally, the third experiment shows what happens when the user does a complete
rotation and faces away from all of the three goals. Again the algorithm
performs quite intuitively and proclaims the user ”irrational” as all the
possible goals were completely on the other side.

Additionally, we tested simple interactions between an industrial manipulator
and the human user. In the first one the robot was selecting goals and
randomly. Should the goal intention estimation detect that the human is moving
to the same goal the robot would stop and select a new goal. Additionally we
used the same framework to navigate the manipulator to the estimated goal of
the human, proving that the framework can also be beneficial in teleportation
scenarios.

Figure 7. The interaction setup with an industrial robot, the three spatial
goals are represented by colour and shape. In this experiment the user took
direct control of the robot and the human intention estimation is used to
detect to which object does the user wish the robot to move, illustrating
another use case of intention estimation. Figure 8. The effect of increasing
the parameter $\beta$. The parameter couple $\beta$ and $\gamma$ set the
threshold for estimating intentions for each goal location. A low $\beta$ will
make the algorithm estimate the unknown intention too much rendering it
unusable , while a too high $\beta$ lowers the general certainty but
eliminates the insecurities between transitions which is an unwanted
behaviour.

##  4\. Conclusion and Impact

In this work we presented a completely portable, robot agnostic system for
intention estimation and visualisation for human-robot collaboration
scenarios. Our system does not require any special set-up or sensors on or
around the robot and is capable of both estimating the human coworkers goals
and actions as well as visualise the goals and intentions of the robot
coworker.

Having an intention estimation system, in addition to explicit intention
declarations, can greatly reduce the mental and physical workload on the user,
while providing constant, information rich data to the robot, thereby
improving the safety and efficiency of robot motions.

We have shown that the goal prediction part of the HAIR system works as
intended and indeed the goal intentions estimated follow a reasoning that
humans might find intuitive and agree with.

Predicting the goal and motion of the human coworker can increase both safety
for the human and the efficiency of robot motions. The goal estimation system
(Petković et al., 2019) was integrated into a mobile robot fleet management
system of a simulated automated warehouse. In (Petković et al., 2020) it was
shown that the proposed system markedly improved warehouse efficiency compared
to no goal estimation or even a simplistic one. It is to be expected that such
an efficiency increase would also be observed in interactions with a robotic
manipulator. Further testing is going to be needed however, to support that
claim.

Likewise the action estimation component as well as the entire system needs to
be evaluated in user studies. More specifically the change in mental and
physical workloads between various intention sharing modalities is of great
interest and quite important in proving the claims that the intention
estimation algorithms presented here significantly decrease the workload
compared to explicitly stating the goals.

As HMDs become ever more common, and the amount of robot coworkers per human
coworker continues to increase, having intuitive HRI using systems that are
cheap, simple and portable becomes essential. Lowering the complexity and
price of each robot by exploiting wearables will lead to a wider use of robots
and increased human-robot collaboration. We hope that the research presented
here provides the first stepping stones towards such a system. HAIR should
raise the quality of interaction between robots and humans, instead of such
interactions raising the hair on the necks of the human coworkers.

Figure 9. Three tests with different goal order. On the left, the user was
selecting goals left to right - cylinder, cube then sphere. In the middle the
user starts with the cube, moves to the cylinder, back to the cube and finally
goes to the sphere. In the test on the right the starting goal is sphere, then
cube, then cylinder, after which the user turns around completely and ends
back on the sphere.

## References

  * (1)
  * Bascetta et al. (2011) L. Bascetta, G. Ferretti, P. Rocco, H. Ardö, H. Bruyninckx, E. Demeester, and E. Di Lello. 2011.  Towards safe human-robot interaction in robotic cells: An approach based on visual tracking and intention estimation. In _2011 IEEE/RSJ International Conference on Intelligent Robots and Systems_. 2971–2978.  https://doi.org/10.1109/IROS.2011.6094642
  * Bochkovskiy et al. (2020) Alexey Bochkovskiy, Chien-Yao Wang, and Hong-Yuan Mark Liao. 2020.  YOLOv4: Optimal Speed and Accuracy of Object Detection.  _arXiv preprint arXiv:2004.10934_ (2020).
  * Chakraborti et al. (2018) T. Chakraborti, S. Sreedharan, A. Kulkarni, and S. Kambhampati. 2018.  Projection-Aware Task Planning and Execution for Human-in-the-Loop Operation of Robots in a Mixed-Reality Workspace. In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. 4476–4482.  https://doi.org/10.1109/IROS.2018.8593830
  * Dragan et al. (2015) A. D. Dragan, S. Bauman, J. Forlizzi, and S. S. Srinivasa. 2015.  Effects of Robot Motion on Human-Robot Collaboration. In _2015 10th ACM/IEEE International Conference on Human-Robot Interaction (HRI)_. 51–58.
  * Forney (1973) G. D. Forney. 1973\.  The viterbi algorithm.  _Proc. IEEE_ 61, 3 (1973), 268–278.  https://doi.org/10.1109/PROC.1973.9030
  * Krupke et al. (2018) D. Krupke, F. Steinicke, P. Lubos, Y. Jonetzko, M. Görner, and J. Zhang. 2018\.  Comparison of Multimodal Heading and Pointing Gestures for Co-Located Mixed Reality Human-Robot Interaction. In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. 1–9.  https://doi.org/10.1109/IROS.2018.8594043
  * Mueller et al. (2018) Franziska Mueller, Florian Bernard, Oleksandr Sotnychenko, Dushyant Mehta, Srinath Sridhar, Dan Casas, and Christian Theobalt. 2018.  GANerated Hands for Real-Time 3D Hand Tracking from Monocular RGB. In _Proceedings of Computer Vision and Pattern Recognition (CVPR)_. 11 pages.  https://handtracker.mpi-inf.mpg.de/projects/GANeratedHands/
  * Ostanin et al. (2020) M. Ostanin, S. Mikhel, A. Evlampiev, V. Skvortsova, and A. Klimchik. 2020\.  Human-robot interaction for robotic manipulator programming in Mixed Reality. In _2020 IEEE International Conference on Robotics and Automation (ICRA)_. 2805–2811.  https://doi.org/10.1109/ICRA40945.2020.9196965
  * Petković et al. (2020) Tomislav Petković, Jakub Hvězda, Tomáš Rybecký, Ivan Marković, Miroslav Kulich, Libor Přeučil, and Ivan Petrović. 2020\.  Human Intention Recognition for Human Aware Planning in Integrated Warehouse Systems.  arXiv:2005.11202 [cs.RO]
  * Petković et al. (2019) Tomislav Petković, David Puljiz, Ivan Marković, and Björn Hein. 2019.  Human intention estimation based on hidden Markov model motion validation for safe flexible robotized warehouses.  _Robotics and Computer-Integrated Manufacturing_ 57 (2019), 182 – 196.  https://doi.org/10.1016/j.rcim.2018.11.004
  * Puljiz et al. (2020) David Puljiz, Franziska Krebs, Fabian Bösing, and Björn Hein. 2020.  What the HoloLens Maps Is Your Workspace: Fast Mapping and Set-up of Robot Cells via Head Mounted Displays and Augmented Reality.  _arXiv preprint arXiv:2005.12651_ (2020).
  * Puljiz et al. (2019) David Puljiz, Katharina S Riesterer, Björn Hein, and Torsten Kröger. 2019.  Referencing between a Head-Mounted Device and Robotic Manipulators. In _Proceedings of the 2nd Workshop on Virtual, Mixed and Augmented Reality Human.Robot Interaction, HRI 2019_.  http://arxiv.org/abs/1904.02480
  * Quintero et al. (2018) C. P. Quintero, S. Li, M. K. Pan, W. P. Chan, H. F. Machiel Van der Loos, and E. Croft. 2018.  Robot Programming Through Augmented Trajectories in Augmented Reality. In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. 1838–1844.  https://doi.org/10.1109/IROS.2018.8593700
  * Walker et al. (2018) Michael Walker, Hooman Hedayati, Jennifer Lee, and Daniel Szafir. 2018.  Communicating Robot Motion Intent with Augmented Reality. In _Proceedings of the 2018 ACM/IEEE International Conference on Human-Robot Interaction_ (Chicago, IL, USA) _(HRI ’18)_. Association for Computing Machinery, New York, NY, USA, 316–324.  https://doi.org/10.1145/3171221.3171253
  * Williams et al. (2019) T. Williams, M. Bussing, S. Cabrol, E. Boyle, and N. Tran. 2019.  Mixed Reality Deictic Gesture for Multi-Modal Robot Communication. In _2019 14th ACM/IEEE International Conference on Human-Robot Interaction (HRI)_. 191–201.  https://doi.org/10.1109/HRI.2019.8673275


#  CO Adsorption on the Surface of MgO from Periodic Coupled-Cluster Theory
with Local Natural Orbitals: Adding to the consensus

Hong-Zhou Ye<EMAIL_ADDRESS>Department of Chemistry, Columbia University,
New York, NY 10027 USA  Timothy C. Berkelbach<EMAIL_ADDRESS>Department of Chemistry, Columbia University, New York, NY 10027 USA

###### Abstract

Accurate determination of the adsorption energy of CO on the MgO (001) surface
has been a challenge for both computations and experiments over the past three
decades. A recent computational study by Shi and co-workers
(10.26434/chemrxiv-2023-h4czl) reported good agreement within $11$ meV ($1$
kJ/mol) between two popular theoretical methods: coupled-cluster with singles,
doubles, and perturbative triples [CCSD(T)] and diffusion Monte Carlo. In this
short note, we report results on the same problem from periodic Gaussian-based
MP2, CCSD, and CCSD(T), with the latter two performed using a recently
developed extension of the local natural orbital (LNO) approximation to
problems with periodic boundary conditions. Our final periodic LNO-CCSD(T)
adsorption energy ($-198\pm 11$ meV) is in quantitative agreement with the
embedded cluster-based LNO-CCSD(T) result ($-199\pm 11$ meV) by Shi and co-
workers. The computational cost of our periodic LNO-CCSD(T) calculations is
comparable to that of the embedded cluster-based LNO-CCSD(T) and is 10 times
less expensive than the plane-wave-based periodic canonical CCSD(T) or 50
times less expensive than the DMC calculations reported by Shi and co-workers.
Our findings highlight the accuracy and computational efficiency of the
periodic LNO-based approach for the simulation of surface chemistry with
correlated wavefunction methods.

##  I Introduction

A recent preprint by Shi and co-workers computationally studied the adsorption
of a single CO molecule on the MgO (001) surface Shi  _et al._ (2023). This is
an intriguing problem due to an enduring challenge of achieving a consensus
between theory and experiment regarding the adsorption energy,
$E_{\textrm{ads}}$, over the past three decades Sauer (2019). In ref 1, the
authors reported a remarkable agreement within $11$ meV or $1$ kJ/mol between
calculations using different theoretical methods, including coupled-cluster
theory with single, double, and perturbative triple excitations Raghavachari
_et al._ (1989) [CCSD(T)], commonly known as the “gold standard” of quantum
chemistry, and diffusion Monte Carlo Foulkes  _et al._ (2001) (DMC). Moreover,
these results were obtained using different computational frameworks,
including periodic calculations with plane wave basis functions Moerman  _et
al._ (2022) and an embedded cluster calculation with Gaussian basis functions
Shi  _et al._ (2022). Their best theoretical estimate of $E_{\textrm{ads}}$ is
$-199\pm 11$ meV, which was obtained by cluster-based local natural orbital
(LNO)-CCSD(T) Rolik and Kállay (2011); Rolik  _et al._ (2013); Nagy and Kállay
(2017); Shi  _et al._ (2022) and which agrees reasonably well with early
experiments Wichtendahl  _et al._ (1999); Dohnálek  _et al._ (2001).

In this short note, we show that results in quantitative agreement with those
in ref 1 can be obtained using a periodic Gaussian-based approach. Our results
demonstrate how recent developments in density fitting Ye and Berkelbach
(2021a, b); Bintrim, Berkelbach, and Ye (2022), correlation-consistent
Gaussian basis sets Ye and Berkelbach (2022), and periodic local correlation
theories Ye and Berkelbach ; Ye and Berkelbach (2023) enable fast and reliable
convergence along all computational axes necessary for simulating the
electronic ground state of condensed-phase systems with up to 100 atoms per
unit cell Bintrim, Berkelbach, and Ye (2022); Ye and Berkelbach . Our main
results are summarized in TABLE 1. The final adsorption energy from our
periodic LNO-CCSD(T) Ye and Berkelbach ; Ye and Berkelbach (2023) calculations
is $-198\pm 11$ meV and agrees almost perfectly with the cluster-based LNO-
CCSD(T) result from ref 1. Importantly, the computational cost of our periodic
LNO-CCSD(T) calculations is comparable to the cluster-based LNO-CCSD(T) but is
about 10 times less expensive than the plane-wave-based periodic canonical
CCSD(T) and about 50 times less expensive than the DMC calculations reported
in ref 1. Comparing our numbers with other independent estimates of
$E_{\textrm{ads}}$ at a similar level of theory from recent literature
highlights the challenges of fully converging the calculations even for
methods like MP2, which is often considered computationally inexpensive.

Table 1: Our final interaction energy ($E_{\textrm{int}}$) and adsorption energy ($E_{\textrm{ads}}$) for CO on MgO obtained through various wavefunction methods compared to recent results in literature. Like the cluster-based approach in ref 1, the main source of uncertainties in our $E_{\textrm{ads}}$ arises from the DFT-determined geometry relaxation energy $\Delta_{\textrm{geom}}$ [eqn 1]. Computational cost (measured in kCPU hours) is also shown when data are available.  Method | Comput. details | $E_{\textrm{int}}/\textrm{meV}$ | $E_{\textrm{ads}}/\textrm{meV}$ |  Cost111Our calculations were run on the AMD EPYC 7302 16-Core Processor (2.994 GHz). | Reference
---|---|---|---|---|---
HF | periodic | $+17\pm 2$ | $+28\pm 11$ | $<0.1$ | this work
MP2 | periodic | $-198\pm 3$ | $-187\pm 11$ | $0.1$
CCSD | periodic, LNO | $-171\pm 4$ | $-160\pm 11$ | $11$
CCSD(T) | periodic, LNO | $-209\pm 4$ | $-198\pm 11$ | $18$
MP2 | cluster embedding | $-200\pm 5$ | $-192\pm 11$ |  | ref 1 (2023)
CCSD(T) | cluster embedding, LNO | $-207\pm 6$ | $-199\pm 11$ | $20$
CCSD | periodic, FNO222Frozen natural orbital | $-153\pm 22$ | $-145\pm 24$ |
CCSD(T) | periodic, FNO | $-201\pm 22$ | $-193\pm 24$ | $200$
DMC | periodic | $-196\pm 24$ | $-188\pm 26$ | $1000$
MP2444Basis set superposition error (BSSE) is not corrected. | periodic, DMET333Density matrix embedding theory | $-431\pm\phantom{1}?$ |  |  | ref 18 (2022)
CCSD444Basis set superposition error (BSSE) is not corrected. | periodic, DMET | $-398\pm\phantom{1}?$ |  |
MP2 | cluster embedding | $+60\pm\phantom{1}?$ |  |  | ref 19 (2016)
CCSD | cluster embedding | $+90\pm\phantom{1}?$ |  |
CCSD(T) | cluster embedding | $+70\pm\phantom{1}?$ |  |

##  II Computational Details

Following ref 1, we calculate the adsorption energy,

$E_{\textrm{ads}}=E_{\textrm{int}}+\Delta_{\textrm{geom}}$ (1)

where $E_{\textrm{int}}$ is the (adiabatic) interaction energy between CO and
MgO calculated with their respective geometries fixed to be those in the
MgO+CO composite system, and $\Delta_{\textrm{geom}}$ is the geometry
relaxation energy. We obtained equilibrium geometries for CO, MgO, and MgO+CO
using density functional theory Hohenberg and Kohn (1964) (DFT) with the
Perdew-Burke-Ernzerhof (PBE) functional Perdew, Burke, and Ernzerhof (1996)
and the D3 dispersion correction Grimme  _et al._ (2010) as implemented in
Quantum Espresso Giannozzi  _et al._ (2009, 2017). As shown in TABLE 2, PBE+D3
gives a lattice constant of bulk MgO and a relaxation energy
$\Delta_{\textrm{geom}}$ that agree well with ref 1, which used a slightly
different DFT protocol, revPBE+D4 Zhang and Yang (1998); Caldeweyher,
Bannwarth, and Grimme (2017). However, PBE+D3 predicts an Mg-C distance that
is slightly shorter, and so we adjust it to 2.460 Å to be consistent with ref
1. The PBE+D3 interaction energy $E_{\textrm{int}}$ is too large by nearly
$100$ meV compared to both revPBE+D4 and other more accurate methods in ref 1.

Table 2: Comparison of DFT results obtained using PBE+D3 (this work) and revPBE+D4 (ref 1). | lat. const. 111Unit: $\mathrm{\AA}$ |  $d(\textrm{Mg-C})$ 111Unit: $\mathrm{\AA}$ |  $\Delta_{\textrm{geom}}$ 222Unit: meV |  $E_{\textrm{ads}}$ 222Unit: meV
---|---|---|---|---
PBE+D3 | $4.222$ | $2.377$ | $11.4$ | $-306$
revPBE+D4 | $4.220$ | $2.460$ | $8$ | $-207$

In what follows, we use our DFT-determined $\Delta_{\textrm{geom}}=11.4\pm 10$
meV, where the uncertainty is taken from ref 1 (which is based on variations
in $\Delta_{\textrm{geom}}$ calculated using different DFT methods), and we
calculate the interaction energy $E_{\textrm{int}}$ using a series of
wavefunction methods as implemented in PySCF Sun  _et al._ (2018, 2020). The
GTH-cc-pV$X$Z basis sets Ye and Berkelbach (2022), augmented by diffuse
functions for the CO molecule and nearby surface atoms, are employed with the
GTH pseudopotential optimized for HF Goedecker, Teter, and Hutter (1996);
Hartwigsen, Goedecker, and Hutter (1998); Hutter (2019). All interaction
energies are corrected for basis set superposition error. The two-electron
integrals are treated by the range-separated density fitting algorithm Ye and
Berkelbach (2021a, b). The integrable divergence of the HF exchange is treated
using a Madelung constant correction Paier  _et al._ (2006); Broqvist,
Alkauskas, and Pasquarello (2009); Sundararaman and Arias (2013). We first use
slab models of two atomic layers (2L) to establish the protocol for converging
$E_{\textrm{int}}(A,X)$ with both the surface size $A$ and the basis set size
$X$ to their respective limits: the infinite surface (IS) limit, $A_{\infty}$,
and the complete basis set (CBS) limit, $X_{\infty}$. We then apply the
protocol to thicker slab models to obtain our final estimate of
$E_{\textrm{int}}$ and hence $E_{\textrm{ads}}$ at both the CBS limit and the
thermodynamic limit (TDL). The results are summarized in TABLEs 1 and 3.

Table 3: The interaction energy (1) for CO on MgO in both the IS limit ($A_{\infty}$) and the CBS limit ($X_{\infty}$) obtained using HF, MP2, CCSD, and CCSD(T) with the respective protocols (2), (5), and (7) for 2L – 4L slabs. The IS corrections $\Delta_{\textrm{IS}}^{\textrm{HF}}$ and $\Delta_{\textrm{IS}}^{\textrm{MP2(FC),corr}}$ evaluated for 2L are used for all slab sizes. | $2\textrm{L}$ | $3\textrm{L}$ | $4\textrm{L}$
---|---|---|---
$E_{\textrm{int}}^{\textrm{HF}}(3\times 3,\textrm{TZ})$ | $+23.1$ | $+20.2$ | $+20.9$
$\Delta_{\textrm{CBS}}^{\textrm{HF}}$ | $-1.9$ | $-0.8$ | $-0.4$
$\Delta_{\textrm{IS}}^{\textrm{HF}}$ | $-3.6$ |  |
$E_{\textrm{int}}^{\textrm{HF}}(A_{\infty},X_{\infty})$ | $+17.6$ | $+15.8$ | $+16.9$
$E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{2\times 2,3\times 3},X_{\textrm{DZ},\textrm{TZ}})$ | $-206.5$ | $-206.7$ | $-207.4$
$\Delta_{\textrm{CBS}}^{\textrm{MP2(FC),corr}}$ | $-4.1$ | $-4.3$ | $-3.9$
$\Delta_{\textrm{IS}}^{\textrm{MP2(FC),corr}}$ | $+7.1$ |  |
$\Delta_{\textrm{FC}}^{\textrm{MP2(FC),corr}}$ | $-12.3$ | $-12.3$ | $-12.6$
$E_{\textrm{int}}^{\textrm{MP2}}(A_{\infty},X_{\infty})$ | $-198.2$ | $-200.5$ | $-199.9$
$\Delta_{\textrm{CCSD}}^{\textrm{LNO-CC(FC)}}$ | $+26.9$ |  |
$E_{\textrm{int}}^{\textrm{CCSD}}(A_{\infty},X_{\infty})$ | $-171.3$ |  |
$\Delta_{\textrm{CCSD(T)}}^{\textrm{LNO-CC(FC)}}$ | $-10.9$ |  |
$E_{\textrm{int}}^{\textrm{CCSD(T)}}(A_{\infty},X_{\infty})$ | $-209.1$ |  |

###  II.1 HF

Figure 1: Convergence of HF interaction energy,
$E_{\textrm{int}}^{\textrm{HF}}$, with surface and basis set sizes for CO on
2L MgO slabs.

We start with the HF interaction energy $E_{\textrm{int}}^{\textrm{HF}}$. FIG.
1 shows the fast convergence of $E_{\textrm{int}}^{\textrm{HF}}(A,X)$ with
both parameters for 2L slabs: using a $4\times 4$ surface with the QZ basis
set essentially reaches both the CBS limit and the IS limit, giving
$E_{\textrm{int}}^{\textrm{HF}}(A_{\infty},X_{\infty})\approx 18.6\pm 1$ meV
for the 2L slab, where half the difference between TZ and QZ results is taken
as a (conservative) estimate of the error bar. While this protocol ($4\times
4$/QZ) can be applied to thicker slabs to probe
$E_{\textrm{int}}^{\textrm{HF}}$ in the TDL, it is computationally expensive.
Given the small finite-size error of $A=3\times 3$ and the nearly perfect
parallelism between TZ and QZ data for all surface sizes, we propose a
computationally efficient alternative

$\displaystyle E_{\textrm{int}}^{\textrm{HF}}(A_{\infty},X_{\infty})\approx
E_{\textrm{int}}^{\textrm{HF}}(3\times
3,\textrm{TZ})+\Delta_{\textrm{CBS}}^{\textrm{HF}}+\Delta_{\textrm{FS}}^{\textrm{HF}}$
(2a)
$\displaystyle\Delta_{\textrm{CBS}}^{\textrm{HF}}=E_{\textrm{int}}^{\textrm{HF}}(2\times
2,\textrm{QZ})-E_{\textrm{int}}^{\textrm{HF}}(2\times 2,\textrm{TZ})$ (2b)
$\displaystyle\Delta_{\textrm{IS}}^{\textrm{HF}}=E_{\textrm{int}}^{\textrm{HF}}(4\times
4,\textrm{DZ})-E_{\textrm{int}}^{\textrm{HF}}(3\times 3,\textrm{DZ})$ (2c)

For the 2L slab, eqn 2 gives
$E_{\textrm{int}}^{\textrm{HF}}(A_{\infty},X_{\infty})=17.6$ meV, which
differs from the result above by only 1 meV. We take this difference to be the
error bar of the protocol (2). TABLE 3 lists
$E_{\textrm{int}}^{\textrm{HF}}(A_{\infty},X_{\infty})$ obtained using
protocol (2) for slabs of 2L to 4L, which are seen to all agree with each
other within our error bar of about $2$ meV. We thus obtain our final HF
interaction energy, $E_{\textrm{int}}^{\textrm{HF}}\approx 17\pm 2$ meV.

###  II.2 MP2

Figure 2: Convergence of the correlation part of the frozen-core MP2
interaction energy, $E_{\textrm{int}}^{\textrm{MP2(FC),corr}}$, with surface
and basis set sizes for CO on 2L MgO slabs. (D,T) denotes CBS extrapolation
using DZ and QZ results [similar for (T,Q)].

Transitioning to correlated methods, we first obtain the converged interaction
energy for MP2. FIG. 2 shows the correlation part of the MP2 interaction
energy evaluated with electrons in the [Ne] core of Mg being frozen (FC),
$E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A,X)$, for 2L slabs. The convergence
of $E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A,X)$ with both parameters is
slower than that of HF, but reliable extrapolations to both limits can be
performed based on the asymptotic behaviors

$\displaystyle E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A,X)$
$\displaystyle\approx
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X)+c_{1}A^{-1}$ (3a)
$\displaystyle E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A,X)$
$\displaystyle\approx
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A,X_{\infty})+c_{2}X^{-3}$ (3b)

A good estimate of
$E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X_{\infty})$ is obtained
by extrapolating the surface size with $A=3\times 3$ and $4\times 4$ and the
basis set size with TZ ($X=3$) and QZ ($X=4$), which gives
$E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{3\times 3,4\times
4},X_{\textrm{TZ},\textrm{QZ}})=-202.1$ meV. An error bar of $2$ meV is chosen
as half the difference between the (DZ,TZ) and the (TZ,QZ) extrapolated CBS
results.

Like for the HF case, the nearly perfect parallelism between results of
different basis sets justifies a computationally efficient alternative to the
protocol above

$\displaystyle
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X_{\infty})\approx
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{2\times 2,3\times
3},X_{\textrm{DZ},\textrm{TZ}})$
$\displaystyle\phantom{E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X_{\infty})}\qquad{}+\Delta^{\textrm{MP2(FC),corr}}_{\textrm{CBS}}+\Delta^{\textrm{MP2(FC),corr}}_{\textrm{IS}}$
(4a)
$\displaystyle\Delta^{\textrm{MP2(FC),corr}}_{\textrm{CBS}}=E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(2\times
2,X_{\textrm{TZ},\textrm{QZ}})$
$\displaystyle\phantom{\Delta^{\textrm{MP2(FC),corr}}_{\textrm{CBS}}}\qquad-
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(2\times
2,X_{\textrm{DZ},\textrm{TZ}})$ (4b)
$\displaystyle\Delta^{\textrm{MP2(FC),corr}}_{\textrm{IS}}=E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{3\times
3,4\times 4},X_{\textrm{DZ},\textrm{TZ}})$
$\displaystyle\phantom{\Delta^{\textrm{MP2(FC),corr}}_{\textrm{IS}}}\qquad-
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{2\times 2,3\times
3},X_{\textrm{DZ},\textrm{TZ}})$ (4c)

For the 2L slab, protocol (4) gives
$E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X_{\infty})\approx-203.5$
meV, which differs from the result above by less than $2$ meV. We take this
difference as the uncertainty of the protocol (4). Finally, we account for the
error of freezing the $[2s^{2}2p^{6}]$ semicore electrons of Mg by a composite
correction

$\displaystyle
E_{\textrm{int}}^{\textrm{MP2,corr}}(A_{\infty},X_{\infty})\approx
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{\infty},X_{\infty})+\Delta_{\textrm{FC}}^{\textrm{MP2(FC),corr}}$
(5a)
$\displaystyle\Delta_{\textrm{FC}}^{\textrm{MP2(FC),corr}}=E_{\textrm{int}}^{\textrm{MP2,corr}}(2\times
2,X_{\textrm{DZ},\textrm{TZ}})$
$\displaystyle\phantom{\Delta_{\textrm{FC}}^{\textrm{MP2(FC),corr}}=}\qquad-
E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(2\times
2,X_{\textrm{DZ},\textrm{TZ}})$ (5b)

TABLE 3 lists the final MP2 interaction energy,
$E^{\textrm{MP2}}_{\textrm{int}}$, calculated using protocol (5) for slabs
from 2L to 4L. Like in the HF case, the TDL of
$E^{\textrm{MP2}}_{\textrm{int}}$ is essentially reached by the 2L model, and
our final estimate of $E^{\textrm{MP2}}_{\textrm{int}}$ is $-200\pm 3$ meV
(the uncertainty from $E_{\textrm{int}}^{\textrm{HF}}$ is included), which
differs from many previous periodic or cluster-based MP2 calculations but
agrees quantitatively with the cluster embedding result from ref 1.

###  II.3 CCSD and CCSD(T)

Figure 3:  Convergence of the correlation part of the frozen-core LNO-CCSD and
LNO-CCSD(T) interaction energy, $E_{\textrm{int}}^{\textrm{LNO-CC(FC),corr}}$,
extrapolated to the CBS limit based on DZ and TZ results for CO on a 2L,
$2\times 2$ MgO slab with respect to the LNO truncation threshold. The hollow
and filled circles denote results obtained before and after applying the
$\Delta_{\textrm{cano}}$ correction (6). The canonical results are shown as
green horizontal lines, with the green shaded area indicating $\pm 2$ meV.

Finally, we calculate $E_{\textrm{int}}$ at the more expensive CCSD and
CCSD(T) levels employing our recently developed periodic extension of the LNO
approximation, which we have used elsewhere to study the dissociation of water
on the surface of Al2O3 and TiO2 Ye and Berkelbach , and whose details will be
described in a separate manuscript Ye and Berkelbach (2023); Clearly, our
periodic LNO-CC method parallels the cluster-based LNO-CCSD(T) method used in
ref 1, but is in some respects simpler because it does not require the
definition of a cluster and embedding protocol, which can be highly non-
trivial for systems like metal oxides Shi  _et al._ (2022). Moreover, periodic
LNO-CC is free of the extra errors introduced at the artificial cluster
boundary. The accuracy of the LNO approximation can be systematically improved
by adjusting a single parameter: the threshold $\eta$ used to truncate the
LNOs by their occupation numbers Rolik and Kállay (2011) (we use a threshold
for unoccupied orbitals that is fixed to be 10 times smaller than that of the
occupied orbitals; all reported thresholds henceforth are for unoccupied
orbitals). As $\eta\to 0$, the LNO-CCSD/CCSD(T) results converge to the result
of canonical CCSD/CCSD(T) calculations. To expedite this convergence, we use a
standard correction based on an MP2 calculation at the same value of $\eta$
Rolik and Kállay (2011).

FIG. 3 shows the $\eta$-convergence of the correlation part of the frozen-core
LNO-CCSD and LNO-CCSD(T) interaction energy, $E_{\textrm{int}}^{\textrm{LNO-
CC(FC),corr}}$, extrapolated to the CBS limit based on DZ and TZ results for
CO on a 2L, $2\times 2$ MgO slab. The system is small enough (about $650$
orbitals in the TZ basis set) to perform canonical CC calculations, results of
which are also shown in FIG. 3 for comparison. We see that the error of LNO-
CCSD and LNO-CCSD(T) (hollow circles) is about $12$ and $6$ meV with a modest
threshold of $10^{-6}$, but an error of $2$ meV requires a tight threshold of
$10^{-7}$ and $10^{-8}$. While LNO-CC calculations with a tight threshold of
$10^{-8}$ are feasible for larger systems, a more efficient alternative is to
correct the bare LNO-CC results by the difference between LNO-CC and canonical
CC in the DZ basis set, which we denoted as

$\begin{split}&\Delta_{\textrm{cano}}^{\textrm{LNO-
CC(FC)}}(\eta)=E_{\textrm{int}}^{\textrm{CC(FC),corr}}(2\times
2,\textrm{DZ})\\\ &\phantom{\Delta_{\textrm{cano}}^{\textrm{LNO-
CC(FC)}}(\eta)}\qquad{}-E_{\textrm{int}}^{\textrm{LNO-CC(FC),corr}}(2\times
2,\textrm{DZ},\eta)\end{split}$ (6)

The corrected LNO-CC results are shown as filled circles in FIG. 3 and seen to
converge much faster than the uncorrected ones.

Figure 4:  Convergence of the final estimate of the CCSD and CCSD(T)
interaction energy obtained using protocol (7) with respect to the LNO
truncation threshold. The best estimate is results from the tightest threshold
and highlighted by the solid horizontal lines, with the shaded area indicating
$\pm 2$ meV.

Our final estimate of the CC interaction energy can be obtained by first
extrapolating $E_{\textrm{int}}^{\textrm{LNO-CC(FC),corr}}(A,X)$ to the IS and
the CBS limits in an approximate manner, chosen to be $A_{2\times 2,3\times
3}$ and $X_{\textrm{DZ},\textrm{TZ}}$, and then correcting the remaining
finite-size and basis set incompleteness error (including the frozen-core
approximation) by a composite correction using the fully converged MP2 result
from above. This can be equivalently formulated as correcting the fully
converged MP2 energy with a CC correction term

$\displaystyle E_{\textrm{int}}^{\textrm{LNO-
CC,corr}}(A_{\infty},X_{\infty},\eta)\approx
E_{\textrm{int}}^{\textrm{MP2,corr}}(A_{\infty},X_{\infty})$
$\displaystyle\phantom{E_{\textrm{int}}^{\textrm{LNO-
CC,corr}}(A_{\infty},X_{\infty},\eta)}\qquad{}+\Delta^{\textrm{LNO-
CC(FC)}}_{\textrm{CC}}(\eta)$ (7a) $\displaystyle\Delta^{\textrm{LNO-
CC(FC)}}_{\textrm{CC}}(\eta)=E_{\textrm{int}}^{\textrm{LNO-
CC(FC),corr}}(A_{2\times 2,3\times 3},X_{\textrm{DZ},\textrm{TZ}},\eta)$
$\displaystyle\phantom{\Delta^{\textrm{LNO-
CC(FC)}}_{\textrm{CC}}(\eta)}\qquad{}-E_{\textrm{int}}^{\textrm{MP2(FC),corr}}(A_{2\times
2,3\times 3},X_{\textrm{DZ},\textrm{TZ}})$
$\displaystyle\phantom{\Delta^{\textrm{LNO-
CC(FC)}}_{\textrm{CC}}(\eta)}\qquad{}+\Delta^{\textrm{LNO-
CC(FC)}}_{\textrm{cano}}(\eta)$ (7b)

The $\eta$-convergence of protocol (7) is illustrated in FIG. 4. Our best
estimate of the CCSD and CCSD(T) interaction energy (including
$E_{\textrm{int}}^{\textrm{HF}}$), obtained from $\eta=10^{-7}$, is $-171.3$
meV and $-209.1$ meV, also listed in TABLE 3. An estimated error bar for these
numbers, estimated using the difference between the two calculations with the
tightest thresholds shown in FIG. 4, is about $2$ meV. Combined with the error
bar of $E_{\textrm{int}}^{\textrm{MP2}}$ ($3$ meV) and of
$\Delta^{\textrm{LNO-CC(FC)}}_{\textrm{cano}}$ ($2$ meV) used in protocol (7),
the final error bar for our CCSD and CCSD(T) numbers is about $4$ meV.

##  III Conclusion

To summarize, we have presented the converged calculation of the adsorption
energy of CO on the surface of MgO using periodic calculations with Gaussian
orbitals. CCSD and CCSD(T) results were obtained using our recently developed
periodic LNO-CC code, which is complementary to both the canonical periodic
plane-wave method and the cluster embedding LNO-CC method used in ref 1, and
the numerical agreement is excellent. Just as for molecules, LNO-CC
calculations for periodic solids can be performed at a fraction of the cost of
canonical periodic CC calculations, thus significantly expanding the scope of
applicability for CC-based methods in periodic simulations.

The results presented here were generated within a few days from scratch,
after learning of ref 1, suggesting that periodic Gaussian-based correlated
theories are on their way to becoming routine in the accurate study of surface
chemistry and other condensed-phase phenomena. We strongly echo the conclusion
of ref 1, that the ability to achieve quantitative agreement between different
methods using different codes and different basis sets, for a nontrivial
problem in surface science, marks an important milestone for the community.

## References

  * Shi  _et al._ (2023) B. Shi, A. Zen, V. Kapil, P. Nagy, A. Grüneis, and A. Michaelides, Chemrxiv preprint (2023), 10.26434/chemrxiv-2023-h4czl.
  * Sauer (2019) J. Sauer, Acc. Chem. Res. 52, 3502 (2019).
  * Raghavachari  _et al._ (1989) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 (1989).
  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Moerman  _et al._ (2022) E. Moerman, F. Hummel, A. Grüneis, A. Irmler, and M. Scheffler, J. Open Source Softw. 7, 4040 (2022).
  * Shi  _et al._ (2022) B. X. Shi, V. Kapil, A. Zen, J. Chen, A. Alavi, and A. Michaelides, J. Chem. Phys. 156, 124704 (2022).
  * Rolik and Kállay (2011) Z. Rolik and M. Kállay, J. Chem. Phys. 135, 104111 (2011).
  * Rolik  _et al._ (2013) Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, J. Chem. Phys. 139, 094105 (2013).
  * Nagy and Kállay (2017) P. R. Nagy and M. Kállay, J. Chem. Phys. 146, 214106 (2017).
  * Wichtendahl  _et al._ (1999) R. Wichtendahl, M. Rodriguez-Rodrigo, U. Härtel, H. Kuhlenbeck, and H.-J. Freund, Phys. Status Solidi A 173, 93 (1999).
  * Dohnálek  _et al._ (2001) Z. Dohnálek, G. A. Kimmel, S. A. Joyce, P. Ayotte, R. S. Smith, and B. D. Kay, J. Phys. Chem. B 105, 3747 (2001).
  * Ye and Berkelbach (2021a) H.-Z. Ye and T. C. Berkelbach, J. Chem. Phys. 154, 131104 (2021a).
  * Ye and Berkelbach (2021b) H.-Z. Ye and T. C. Berkelbach, J. Chem. Phys. 155, 124106 (2021b).
  * Bintrim, Berkelbach, and Ye (2022) S. J. Bintrim, T. C. Berkelbach, and H.-Z. Ye, J. Chem. Theory Comput. 18, 5374 (2022).
  * Ye and Berkelbach (2022) H.-Z. Ye and T. C. Berkelbach, J. Chem. Theory Comput. 18, 1595 (2022).
  * (16) H.-Z. Ye and T. C. Berkelbach, arXiv:2309.14640 (2023) .
  * Ye and Berkelbach (2023) H.-Z. Ye and T. C. Berkelbach, in preparation (2023).
  * Mitra  _et al._ (2022) A. Mitra, M. R. Hermes, M. Cho, V. Agarawal, and L. Gagliardi, J. Phys. Chem. Lett. 13, 7483 (2022).
  * Mazheika and Levchenko (2016) A. Mazheika and S. V. Levchenko, J. Phys. Chem. C 120, 26934 (2016).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys. Condens. Matter 29, 465901 (2017).
  * Zhang and Yang (1998) Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
  * Caldeweyher, Bannwarth, and Grimme (2017) E. Caldeweyher, C. Bannwarth, and S. Grimme, J. Chem. Phys. 147, 034112 (2017).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, Wiley Interdiscip. Rev. Comput. Mol. Sci 8, e1340 (2018).
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020).
  * Goedecker, Teter, and Hutter (1996) S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996).
  * Hartwigsen, Goedecker, and Hutter (1998) C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev. B 58, 3641 (1998).
  * Hutter (2019) J. Hutter, “New optimization of GTH pseudopotentials for PBE, SCAN, PBE0 functionals. GTH pseudopotentials for Hartree-Fock. NLCC pseudopotentials for PBE.” https://github.com/juerghutter/GTH (2019).
  * Paier  _et al._ (2006) J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber, and J. G. Ángyán, J. Chem. Phys. 124, 154709 (2006).
  * Broqvist, Alkauskas, and Pasquarello (2009) P. Broqvist, A. Alkauskas, and A. Pasquarello, Phys. Rev. B 80, 085114 (2009).
  * Sundararaman and Arias (2013) R. Sundararaman and T. A. Arias, Phys. Rev. B 87, 165122 (2013).