This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.

# Adaptive Surgical Robotic Training Using Real-Time Stylistic Behavior
Feedback Through Haptic Cues

Marzieh Ershad, Robert Rege, and Ann Majewicz Fey M. Ershad, the Department of
Electrical Engineering, University of Texas at Dallas, Richardson, TX, 75080.
e-mail<EMAIL_ADDRESS>Majewicz Fey, the Department of
Mechanical Engineering, University of Texas at Austin, Austin, TX 78712A.
Majewicz Fey and R. Rege, the Department of Surgery at UT Southwestern Medical
Center, Dallas, TX, 75390

###### Abstract

Surgical skill directly affects surgical procedure outcomes; thus, effective
training is needed to ensure satisfactory results. Many objective assessment
metrics have been developed and some are widely used in surgical training
simulators. These objective metrics provide the trainee with descriptive
feedback about their performance however, often lack feedback on how to
proceed to improve performance. The most effective training method is one that
is intuitive, easy to understand, personalized to the user and provided in a
timely manner. We propose a framework to enable user-adaptive training using
near-real-time detection of performance, based on intuitive styles of surgical
movements (e.g., fluidity, smoothness, crispness, etc.), and propose a haptic
feedback framework to assist with correcting styles of movement. We evaluate
the ability of three types of force feedback (spring, damping, and spring plus
damping feedback), computed based on prior user positions, to improve
different stylistic behaviors of the user during kinematically constrained
reaching movement tasks. The results indicate that four out of the six styles
studied here were statistically significantly improved (p$<$0.05) using spring
guidance force feedback and a significant reduction in task time was also
found using spring feedback. Path straightness and targeting error were other
task performance metrics studied which were improved significantly using the
spring-damping feedback. This study presents a groundwork for adaptive
training in robotic surgery based on near-real-time human-centric models of
surgical behavior.

###### Index Terms:

Surgical Robotics, Force Feedback, Adaptive and Intelligent Educational
Systems

##  I Introduction

Surgical outcomes are highly dependent on surgeon skill levels. Efficient
training that provides trainees with appropriate feedback and assists them
with achieving expert-like performance is critical for mastering technical
skills in surgery and achieving successful procedural outcomes [1].
Traditional methods in surgical training typically involve an expert observing
and evaluating a trainee’s performance in the operating room, and providing
feedback to the trainee on how to improve performance [2]. Automating skill
assessment can alleviate the time intensiveness and subjectiveness of these
methods; Furthermore, finding an effective and efficient feedback method,
which is intuitive and easy to understand is crucial [3].

For patient-free and more objective training environments, virtual reality
(VR) simulators have begun to find their way into surgical training [4, 5].
Simulators provide factual and quantitative data to the human user upon
completion of each simulated task, such as number of instrument collisions,
time to complete the task, and the number of missed targets. These metrics
indicate the success rate of the trainee but do not necessarily provide them
with meaningful feedback on how to modify their movements to improve
performance [6].

To address this issue, an ongoing development in surgical simulators is to
enable real-time feedback to users based on the calculated metrics. Errors are
computed from the user’s interaction with the virtual environment with respect
to a reference, and to correct this error, feedback is provided to the user
accordingly.

The error can be calculated based on deviating from an expected trajectory or
performance variable. Haptic feedback has been widely used for training
purposes in simulators to assist with following a specified trajectory or
providing a sense of touch in interaction between the tissue and the
instrument while performing a task. An example of haptic feedabck for
providing tissue-tool interaction is the work from Pezzementi et al. They
implemented a platform for interaction with soft tissue in a simulated
environment using the Phantom Omni haptic device by training a linear 2D mass-
spring-damper system which performs similar to a nonlinear finite element (FE)
model [7]. For trajectory guidance, Ko et al. developed a training simulator
to assist the trainee with following a desired catheter insertion path through
haptic feedback by calculating forces during catheter insertion [8]. These
methods have proved to be effective in improving performance however, do not
incorporate user’s behavior or movement style which provide rich information
regarding user’s proficiency, and can lead to more intuitive training methods.

An effective training method should be easily interpretable by the user. In an
earlier study [9], we showed that the quality of movement during task
performance which is intuitively perceived by a human observer can be used to
distinguish different expertise levels; thus the user’s style of movement
includes valuable information regarding his/her skill level, and deviations
from expert-like movements can be used to calculate relevant feedback for
training.

Another recently explored source for improvement in virtual reality surgical
simulators is adaptive training which provides relevant and customized
training feedback to trainees, based on individual strengths and weaknesses,
and could enhance learning outcomes. The large amount of data recorded and
stored by VR simulators enables data-driven analysis and automatic performance
evaluation. This enables adaptive training based on each individual’s
performance [10]. An example of an adaptive robotic surgical training
framework is presented in [11]. This study compares adaptive curriculum
training to self-managed training and shows significant improvement in
performance and learning skill using an adaptive framework. However, these
performance assessment and adaptive feedback methods are largely task-
dependent, which limit the generalizability of these approaches.

In the following, we will discuss previous studies in this field and describe
our proposed methodology which addresses the issues mentioned above to assist
with improving training in robotic surgery. The rest of the paper is
structured as follows. In section II we summarize related work in adaptive
training, force-reflective feedback, and guidance force feedback. In section
III our proposed stylistic assessment and feedback method is discussed in
detail. It includes a deficiency detection phase and a feedback applying
phase. A deficiency in style is detected from user hand position and velocity
data, by comparing to expert style for a variety of stylistic descriptors.
Then, the user is provided with either spring, damping, or spring-damping
force feedback, depending on their randomized assignment of our feedback
groups. We evaluate the effectiveness of our adaptive stylistic force feedback
using both performance metrics as well as stylistic changes over the duration
of the experimental study. Section IV describes the experiment design and
tools used to conduct the experiment. In section V, we present the results of
the proposed training method, and discuss the effect of the different types
force feedback on styles of movement. Section VI concludes the paper and
suggests the future research in this field.

##  II Related Work

###  II-A Adaptive Training

Adaptive technology can be introduced into training devices to develop user-
specific training that results in more effective learning. An adaptive system
can be known as a supervisor that instructs each trainee based on his/her
unique performance and provides specific instructions on how to proceed, or
adjusts the training task for each individual to ensure best results for each
trainee. These systems consist of a control loop that detect changes in the
output from a desired point, this can be done using machine learning
approaches that enable deficiency detection or performance classification .
Feedback is then applied to modify the response to move it towards the desired
performance levels [10]. Adaptive systems require three main elements
including constant monitoring and measurement of performance, an adaptive
variable, and a methodology to adjust the variable to enhance performance
[12].

In adaptive training, the user’s performance is evaluated based on specific
criteria (detection phase) and then in the next step training is adapted
accordingly (feedback or training phase). Task difficulty level is one element
of focus in the training phase. The difficulty of the task can be updated
based on user’s performance to adjust the level of challenge and enhance
learning. This has been studied in digital games [13]. The stimulus or the
type of feedback provided to the user is another element of focus in the
training phase in which the feedback is adapted based on the user’s
performance. Visual, audio, and haptic feedback are some of types of feedback
used in the training phase. Different types of haptic feedback used in
training systems, will be discussed in the following section.

###  II-B Haptic Feedback for Training

To study the effect of haptic feedback on user’s performance, two types of
haptic feedback are noticeable: reflective feedback and guidance feedback. The
former provides the user with a feeling of touch and force in interacting with
an object in environments where these senses are missing. In virtual
environments, this is done through haptic rendering. The latter provides the
user with haptic cues, and assists in correcting movements to improve
performance.

####  II-B1 Reflective Feedback

Haptic feedback in VR simulators improves training [14]. The lack of haptic
feedback (both force and tactile) causes an inappropriate level of force
applied to the tissue which can lead to safety issues [15]. Tactile feedack
decreases the force applied to the tissue and hence reduces tissue damage. A
study was conducted to show this effect on robotic manipulation by mounting
force feedback system onto a da Vinci surgical robotic system performing
multiple peg transfer tasks [16]. This study showed that higher force was
applied by all subjects in the absence of haptic feedback; thus, indicating
that tactile feedback assists surgeons with tissue handling by applying an
appropriate amount force to the tissue. In another study, the effect of
tactile force feedback was evaluated in vivo [17]. This study also showed a
significant reduction in grasping forces and thus, tissue damage in the
presence of an integrated tactile feedback. A study by Abiri et al. showed
that a multi-modal feedback including tactile, kinesthetic, and vibrotactile
feedback for providing a sense of touch in tissue grasping and manipulation
tasks resulted in an average of 50% reduction in force compared to a no
feedback scenario [18].

This type of reflective feedback though proving to be helpful in providing the
user with a feeling of touch and force in teleoperated environments where
these sense are missing, do not provide any cues to the user on how to modify
movement to improve performance. Guidance feedback which addresses this issue
is discussed in the following.

####  II-B2 Guidance Feedback

Haptic guidance can enhance learning new motor skills in robotic environments
where an instructor is not present to guide the user on how to modify his/her
movement. Different studies have shown the effectiveness of haptic feedback in
developing motor skills [19], as well as movement guidance [20]. A common type
of training motor skills using haptic feedback, is transferring expert skills
in which a an expert’s movements are recorded and played back to train a
novice [21]. However, Gibo et al. showed that haptic feedback can help
discover new movement strategy rather than following a specific trajectory or
enforcing a specific movement [22]. They provided the subjects with an
environment to explore different types of movement using haptic feedback and
adopt the best strategy. Haptic disturbances which suggests disturbing the
movemnet instead of guding the user can also improve motor skills [23].

while all these methods prove the effectiveness of haptic feedback in movement
guidance, they do not focus on performance feedback. Jantscher et al. designed
and implemented a framework that provides vibrotactile feedback method based
on movement smoothness. They proved that the smoothness-based feedback
improved accuracy compared to trajectory based feedback methods [24]. They
provided the subjects with vibrotactile cues with a degree of pleasantness
relative to their performance; however, results in the literature show
scenarios in which force feedback assists the user with performing a task, yet
is perceived negatively by the human user [25], and other scenarios in which
force feedback does not improve performance, yet is preferred [26]. These
results indicate that objective performance metrics and subjective user
response surveys may not be sufficient for understanding the intuitiveness of
a control interface.

Similar to [24], to further investigate performance based haptic feedback, we
examine the effect of haptic cues on six different stylistic performance
behaviors and study three different types of feedback to find the best type
which improves each movement style. We propose a framework to provide task-
independent stylistic feedback to the human user during movement-based
training tasks to provide the user with a more intuitive and global
understanding of their movement styles. We designed, implemented, and
evaluated an adaptive training method composed of the following elements: (1)
Our proposed framework first evaluates the user’s stylistic behavior
performance in near real time and detects deficiencies in some movement styles
[27]. (2) Next, it provides the user with haptic cues to modify their movement
to improve performance. We also evaluate the effectiveness of three common
types of haptic feedback, namely, spring, damping and spring-damping feedback
that is computed from prior user positions and velocities. The goal of our
study is to find intuitive ways to communicate with the user on how to modify
his/her movement to enhance performance. We evaluate the user performance,
based on the quality of movement through monitoring their styles of movement
(movement styles are is described in section III.A) while performing a task.
We then provide haptic feedback to the user to help correct their style in
near-real-time.

##  III Methods

Our goal is to improve robot assisted training to help achieve mastery in
surgical robotics. For this purpose, we aim to (1) introduce a customized
framework in which each individual is provided training based on his/her
performance, (2) provide the trainee with feedback in a timely manner and in
near-real-time, (3) introduce a generalizable and task-independent framework
which evaluates performance based on the user’s style of movement, and (4)
develop a more understandable and intuitive way to communicate with the user
on how to modify movement to improve performance.

A systematic framework for recognizing the quality of movement through
stylistic behavior and applying appropriate feedback for correcting the style
was developed using a human machine interface (i.e., a haptic device) and a
simulated task. Fig. 1 shows the block diagram of the proposed method.

Figure 1: System Block Diagram: The human user interacts with a haptic device
and the simulation environment (a). In advance of the experiment, training
movement data is used to learn a dictionary of stylistic features and a
classifier is trained to predict stylistic deficiencies in near-real-time [28]
(b). During the experiment, kinematic measurements from the haptic device is
represented into stylistic behaviors by projecting it on the learned
dictionary (c). The quality of the user’s style is detected using a classifier
which takes the coefficients of the new representation of the data as an input
(d). Finally, force feedback is provided to the user if negative performance
is detected. Three different types of feedback forces were evaluated in this
study for their effectiveness in improving user style. Feedback is computed
from prior user positions and velocities (e).

###  III-A Surgical Skill Assessment Using Stylistic Behavior

In a previous study, we presented a novel surgical skill assessment method
based on the user’s stylistic behavior [9]. These styles of movement
represents the appearance of movement in action described by common adjectives
used to describe the quality of movement such as smoothness, fluidity,
decisiveness, etc., which is easily distinguishable to a casual observer with
different expertise levels. The idea behind this method is that the quality of
movement holds fundamental information about a subject’s skill; thus,
quantifying these universally understandable movement descriptors enables the
development of effective and intuitive training strategies. We proposed a
lexicon of contrasting adjectives representing surgical styles through
collaboration with expert surgeons (Table I). To evaluate the ability of these
stylistic descriptors in differentiating among different expertise levels, we
used crowd-sourced assessment which has proven to show comparable results to
expertise evaluation. Paired videos of a subject performing a simulated
surgical task and the task being performed was posted to Amazon Mechanical
Turk and the crowd rated the videos based on the stylistic descriptors.

To quantify the qualitative assessment based on stylistic behavior, we found
data metrics associated with each stylistic behavioral adjective in the
lexicon through an extensive search among different calculated metrics. For
each adjective, we found the metric that correlated best with the crowd
ratings. These metrics were calculated from kinematic and physiological
measurements recorded from multiple sensors from the user’s hand movement
while performing a simulated task on the da Vinci surgical simulator.
Furthermore, we evaluated the ability of the stylistic descriptors to
differentiate between different expertise levels. For this purpose, the
metrics associated with the stylistic behavior were used to train a classifier
which was then used to distinguish among four levels of expertise (novice,
intermediate, expert, fellow) [27]. The results showed that these styles of
movement were able to distinguish among different expertise levels.

In the next step, to avoid the feature engineering required in the previous
study for identifying the stylistic behavior and to detect the deficiency in
the styles of movement during user’s performance, we proposed an automatic
method for extracting underlying structures that represent stylistic behavior
from raw kinematic data within 0.25 seconds of movement [28].

In this study, we design an experiment, to implement, and test the framework
for automatically detecting the deficiency in movement styles in near real-
time. In addition, to assist with correcting the style of movement as a ground
work for developing a training framework, we examine the effect of haptic
guidance using three different types of force feedback (spring, damping,
spring and damping) on the six different styles of movement in Table I.

Table I: Lexicon of Stylistic Behavior Positive Adjective | Negative Adjective
---|---
Fluid | Viscous
Smooth | Rough
Crisp | Jittery
Relaxed | Tense
Deliberate | Wavering
Coordinated | Uncoordinated

###  III-B Detecting Deficiencies in Stylistic Behavior

A framework for detecting the stylistic behavior performance is described in
[28]. A similar approach is used in this study, however, a different data set
is used and the model tuned to best fit the new data set. This approach is
discussed in the following.

####  III-B1 Crowd-Sourced Assessment for Positive and Negative Performance
for Each Style

To be able to train a model to recognize a deficiency in movement styles, we
first label the data based on a positive or negative performance of the
stylistic behavior. For this purpose, we use the JIGSAWS data set [29], which
is a publicly available data set that contains robotic surgical training
videos and kinematic recordings. JIGSAWS videos were uploaded to Amazon
Mechanical Turk and crowd workers rated the videos based on the quality of
performance in the six styles of movement mentioned in Table I. The crowd-
workers were asked to rate the video based on either a positive or negative
adjective for a given stylistic descriptor (e.g., smooth v.s. rough movement,
crisp vs. jittery movement). Each video was rated by 20 crowd-workers. The
trial was eventually assigned a positive label if it was rated positive by
more than or equal to 50 % of the crowd-workers and was otherwise assigned as
negative.

####  III-B2 Dictionary Training and Classifier Model Training

For each stylistic behavior, a dictionary containing the basis vectors for the
good and bad performance is created using the kinematic data from the right
hand manipulator of da Vinici skill simulator from the JIGSAWS data set. This
data set includes position, velocity, and angular velocity from the robot end
effectors. A separate dictionary was learned from the positive as well as the
negative performance and then the total dictionary was obtained from the
concatenation of the these two sets of dictionary such that the first half of
the basis vectors were dictionary learned from the good performance and the
second half were the dictionary learned from the bad performance. The positive
and negative labels regarding each stylistic behavior adjective used to train
the model were obtained from crowd-sourced assessment on the JIGSAWS video
data set (section III-B1). The input data is then represented as a linear
combination of the basis vectors in the dictionary. The dictionary and the
coefficients are calculated using an optimization algorithm that iterates
between two problems: 1) finding the basis vectors such that the reconstructed
signal is as similar as possible to the input signal, and 2) finding the
coefficients such that they are sparse. The sparseness reduces the
computational complexity and enables near-real-time implementation. These
sparse codes are then used to train a SVM classifier. Six separate codebooks
are learned for each of the six stylistic behavior adjectives, leading to six
trained classifiers.

####  III-B3 Coefficient Calculation

For a new set of input data (i.e. a new frame of 30 samples), dimensionality
reduction is done using principle component analysis (PCA) to remove
correlations in the data set, then this reduced dimension data set is
projected onto the learned codebook (described in section III-B2). The new
representation of the input signal is sparse. The sparse codes form the new
data frame at each point of time, which are then fed into the trained
classifier (described in section III-B2) for performance evaluation. Algorithm
1 shows the pseudo code for this method.

Algorithm 1 Style Performance Detection Algorithm

Input: new data
Output: stylistic behavior performance $S_{i}$

1: while trial not finished do

2: Get every data frame of 30 samples ($df$)

3: perform PCA on the new data frame ($df_{PCA}$)

4: Project the reduced dimension data set onto the dictionary
$D=[D_{p},D_{N}]$, $D_{p}\in R^{l}$ and $D_{N}\in R^{l}$, l: Number of basis
vectors in the dictionary,
$df_{PCA}=W_{p1}D_{p1},W_{N1}D_{N1}+W_{p2}D_{p2},W_{N2}D_{N2}+$
$...+W_{pl}D_{pl},W_{Nl}D_{Nl}$

5: use the sparse codes ($W_{p1},W_{N1},...,W_{pl},W_{Nl}$) as input to the
pre-trained classifier

6: $S_{i}$ = classifier output $S_{i}$ = 1 if poor performance is detected,
$S_{i}$ = 0 if good performance is detected)

7: end while

###  III-C Providing Feedback for Correcting Stylistic Behavior

To avoid confusing the operator with multiple, potentially competing feedback
cause, the experiment was divided into 6 blocks and only one stylistic
deficiency was detected within this set of movement trials. Based on which
style detection algorithm was activated for a given block in the experiment
protocol, when a poor performance was detected using the proposed near real-
time algorithm, one of the three type of force feedback was turned on. In the
following the three types of force feedback compared in this study (Fig. 2)
are discussed.

Figure 2: Three types of haptic feedback: spring, damping, and spring +
damping feedback were studied here for their ability to provide stylistic cues
to the human operator. A force feedback was generated based on the user’s
prior position in time.

  * •

Spring Feedback: This was calculated using the difference between the position
of the hand at time t $(D_{t})$, and the position at time t-1 $(D_{t-1})$.

$F_{s}=K_{s}(D_{t}-D_{t-1})$ (1)

The gain $K_{s}$ was obtained through trial and error and chosen to be 30. The
gain was chosen to be high enough so that the user would be able to feel the
feedback, but also maintain the stability of the system. This gain was fixed
throughout the experiment.

  * •

Damping Feedback: Was calculated using the difference between the velocity of
the hand at time t $(D_{t})$, and the velocity at time t-1 $(D_{t-1})$.

$F_{d}=B_{d}(V_{1}-V_{2})$ (2)

The gain, $B_{d}$ was chosen through trial and error and was set to be 15. A
lowpass filter with a cutoff frequency of 100 HZ was used to remove noise and
smooth the velocity signal and prevent the system from becoming unstable.

  * •

Spring + Damping Feedback: Was calculated using the difference between the
velocity of the hand at time t $(D_{t})$, and the velocity at time t-1
$(D_{t-1})$.

$F_{s}d=K_{sd}(D_{1}-D_{2})+B_{sd}(V_{1}-V_{2})$ (3)

The gains, $K_{sd}$ and $B_{sd}$ were chosen through trial and error and set
to be 10 and 5. A lowpass filter with a cutoff frequency of 100 HZ was used to
remove noise and smooth the velocity signal and prevent the system from
becoming unstable.

Algorithm 2 shows the pseudo code for the haptic feedback algorithm.

Algorithm 2 Feedback Generation Algorithm

Input: $S_{i}$: style performance , 0 otherwise
$S_{i}=1$ if good performance is detected,
$S_{i}=0$ if bad performance is detected,
Output: $f_{out}$: force feedback to be applied

1: define feedback type (F) $F=F_{s}$: Spring Feedback, $F=F_{d}$: Damping
Feedback, $F=F_{s}d$: Spring+Damping Feedback

2: while trial not finished do

3: if $S_{i}=0$ then

4: $f_{out}=F$

5: end if

6: end while

##  IV Experimental Setup

###  IV-A Data Acquisition and Simulated Task

The Geomagic Touch haptic device (3D Systems, Rock Hill, SC) was used in this
study. This device allows for 3-degree-of-freedom force feedback and 6-degree-
of-freedom sensing. It is used to both provide the user with the desired
movement tasks, as well as force feedback guidance cues based on stylistic
deficiencies. Position, and linear and angular velocity measurements were
recorded from the stylus of the haptic device at a frequency of 256 Hz. To
enable near-real-time performance, stylistic detection was performed on every
frame of 30 samples of incoming data (representing 0.12 seconds). The
simulated task consisted of reaching a set of targets under a kinematically
constrained environment, simulating the control of a steerable needle using
Cartesian Space teleoperation [30]. This task was chosen due to its complexity
as a single-handed movement and one that naturally hunders movement is a
straight-line path, which we felt would not be difficult enough to illicit
stylistic changes in the user movements. The movement tasks were developed
using C++ and the CHAI 3D haptic rendering library. Users were asked to reach
four 5 mm targets, mirrored vertically, at predefined locations which were
presented to the user at random. The user was instructed to initialize each
trial by moving the virtual stylus to the starting point. After reaching the
target the user would end the trail by pressing a button on the stylus. Data
was collected from the time the user initialized the haptic device until they
defined the end of the trial (Fig. 3).

###  IV-B Experimental Protocol

The experiment was divided into six blocks of kinematically constrained
movement trials (e.g., controlling a steerable needle under cartesian space
teleoperation), each block corresponding to one of the six stylistic
adjectives. Each block includes a baseline segment consisting of two
repetitions for each target (a total of 16 reaching trails) with no force
feedback, and a segment that contains an applied force feedback for five
repetitions of movements for each target (a total of 40 reaching trials). This
resulted in 336 trials (6 blocks x 56 trials per block) for each participant.
In each block, force feedback was provided when a stylistic deficiency was
detected for the given adjective corresponding to the block. A 20 sec break
was provided to the user between each block. Both target location and ordering
of stylistic blocks were randomized. Figure 4 shows an example of the
experiment protocol.

###  IV-C Participants

A total of 21 subjects participated in this study. The study protocol was
approved by UTD IRB office (UTD # 14-57). Participants had no previously
reported muscular-skeletal injuries or diseases, or neurological disorders.
The subjects were divided into 3 groups of 7 subjects each. Each group was
assigned the same randomized movement task, but only received either spring,
damping, or spring-damping feedback for each of the stylistic adjective
blocks. This parallel study design was chosen to allow us to evaluate the
effect of the type of haptic feedback on corresponding changes in stylistic
behavior.

(a) Simulated Task (b) Target Layout

Figure 3: (a) User interface: user interacting with simulated environment
using the Geomagic Touch haptic device. The task was initiated by moving the
virtual stylus to the red doughnut and would end by reaching the specified
target. (b) Target layout.
Figure 4: An example of an experiment protocol for one subject. The protocol
consists of six blocks, each related to one stylistic behavior detection
algorithm that was activated for that block. For each block, the user first
performed a set of reaching movements with no feedback to enable a baseline
computation of style, followed by a set of trials with feedback that was
provided, based on measured stylistic deficiencies. For each subject, a single
feedback method was provided throughout the whole experiment, but at different
points of time, depending on the style detection algorithm for that subject.
Hence, a unique feedback relevant to style was provided to each subject.

###  IV-D Stylistic Behavior Performance Detection

The kinematic data recorded from the haptic device includes user hand
position, velocity, and angular velocity all in X, Y, and Z directions,
resulting in 9 signal channels which are similar to the class of signals used
to train the performance detection model from the JIGSAWS dataset, described
in section III-B2. These set of basis vectors are used here for obtaining the
new spars representation of the input signal. The basis vector obtained from a
class of signal similar to the input signal, better enables capturing the
underlying information in the signal as opposed to using predefined
dictionaries (e.g. Fourier, Wavelet, etc.) Based on the style detection
algorithm activated in each block, the new frame of data was projected onto
the set of an over complete dictionary that was calculated as discussed in
III-B. The sparse codes for each incoming frame of data was calculated and
used as input to a classifier to detect the performance quality based on the
activated detection algorithm for the specific style. The classifier returns 0
if a poor performance is detected and returns 1 otherwise. The detection
algorithms were implemented in MATLAB.

###  IV-E Providing Feedback to the User

For each frame of incoming data if a poor performance was detected (output of
the classifier was 0), one type of force feedback (spring force feedback,
damping force feedback, or spring-damping force feedback), was activated and
applied to the user’s hand. A custom C++ code was developed to apply the force
through the Geomagic Touch device. Robot Operating System (ROS) was used to
build the connection between detection algorithm in MATLAB and applying the
force to the user through the Geomagic Touch haptic device in C++. Three types
of forces, as discussed in section III-C, were studied in this experiment.
Each group of subjects was provided with one type of force feedback throughout
the whole experiment for all different blocks of style detection.

###  IV-F User Performance Evaluation Metrics

To quantify the quality of performance in each trial in which a feedback was
applied, the performance quantity $P$ was calculated. For each style (i.e.,
each block in the protocol (Fig 4)), the first section of the block where no
feedback is applied is used as a baseline for that style. For each trial, the
performance of the user was evaluated by the sum of number of times a one was
detected (good performance), divided by the total number of detections in that
trial. This was done for the baseline trials for each style and averaged over
all force-feedback trials for the same style.

$p=\dfrac{\sum_{i}^{N=40}{(num\\_positive\\_WF/num\\_total\\_WF)/N}}{\sum_{j}^{M=16}{(num\\_positive\\_NF/num\\_total\\_NF)/M}}$
(4)

Where: num_positive_WF is the number of good performance detected in a trial
with feedback, num_total_WF is the number of total detections in a trial with
feedback, i is the trial index for feedback trials, and N is the total number
of trials with feedback for one style. In the denominator, we
defined:num_positive_NF as the number of good performance detected in the
baseline trial (no feedback), num_total_NF as the number of total detections
in the baseline trial (no feedback), j as the trial index for baseline trials,
and M as the total number of baseline trials.

###  IV-G Task Performance Evaluation Metrics

To compare the effect of the three types of feedback on the task performance,
three metrics were calculated including: (1) time taken to reach the target,
(2) needle trajectory straightness (the distance traveled by the needle
divided by a straight line to the target), and (3) the needle position error
(the distance between the needle and the target at the end of the trial).

##  V Results and Discussion

We collected a total of 7056 trials (21 subjects, 336 each). Data analysis was
carried out for all trials. The results include the evaluation of stylistic
behavior improvement, as well as an evaluation of task performance as a
function of the different types of haptic force feedback. A NASA Task Load
index survey was also conducted to show how users percieved the feedback
provided to them in terms of workload.

###  V-A Effect of Force Feedback on Styles

The effect of each type of force feedback on each style of movement is shown
in Figure 5. The mean and standard deviation of the quantity associated with
good performances ($P$) for the three different types of force feedback
(spring, damping, spring-damping) are plotted. This is the average number of
good performances detected in the feedback segment of one block (i.e., one
style detection algorithm, normalized to the average of the number of good
performances detected in the baseline, no feedback, segment of the same
block). The values above the horizontal line crossing at 1 show the
improvement of the movement style when applying feedback with respect to the
no feedback condition and the values below this line indicate that receiving
feedback did not improve the movement style compared to not receiving any
feedback.

This plot indicates that the spring force feedback was able to improve the
average performance of the styles “crisp”, “deliberate”, and “relaxed”. The
damping force feedback improved the “crisp” and “deliberate” styles on
average, and the spring+damping force feedback was able to improve the
“smooth”, “calm”, “deliberate” styles on average.

Overall, all styles except for “fluid”, showed an average improvement by
applying one or more types of force feedback. The “fluid” style however showed
the best performance in the absence of the forces studied here. This can be
due to the fact that other kinematic metrics, rather than the position and
velocity, contribute to the fluidness of movement. In this study only force
feedback associated with position and velocity were studied. According to our
previous study [31], the angular velocity of the hand movement is related to
the fluidity of the movement. Thus, in future work, applying other types of
force feedback which incorporate the effect of angular velocity might help to
improve the fluidity of movement. This study was limited by the fact that the
haptic device used was not able to provide rotational feedback cues. The style
“deliberate” was improved by all types of forces when compared to the no
feedback condition; however the most improvement occurred when applying spring
force feedback.

Figure 5: Comparing the effects of three different types of haptic feedback on
each style. For each group of subjects receiving the same type of feedback,
the mean and standard deviation is shown for the number of positive
performance normalized to the total number of detections and divided by the
baseline stylistic positive performance, for each style. The values above 1
show an improvement in the performance when applying feedback compared to the
no feedback condition.

A post hoc statistical analysis was done to determine significant differences
in the three types of force feedback, different targets, and task repetitions
for each stylistic behavior adjective. The normality test to identify a normal
distribution in the data was rejected and thus, the Kruskal Wallis test was
used to identify significantly different groups. Effect significance is
identified for p-values less than 0.05.

The results from the statistical analysis on different styles of movement
(Table II) indicate that for the styles, Fluid/Viscous, Relaxed/Tense,
Deliberate/Wavering, the spring force feedback showed significant difference
in improving the user performance compared to the other two types of feedback.
For the Crisp/Jittery style both the spring feedback and damping feedback
showed significant improvement in performance. For the Calm/Anxious styles,
the spring+damping force feedback showed significant improvement in
performance.

The statistical analysis indicate that task repetition shows no statistically
significant effect on the different types of stylistic behavior as opposed to
the target location which show visible statistically significant effects on
different styles due to different target locations. This is shown in the third
and fourth column of Table II.

Table II: Statistical analysis summary of the effect of different force
feedback types, targets, and task repetitions on the stylistic behavior

| Force Feedback | Target | Repetition
---|---|---|---
Style | $p$ | Significance | $p$ | Significance | $p$ | Significance
Fluid/Viscous |  $<$0.0035 |  S$>$SD, D | 0.64259 | N/A | 0.9916 | N/A
Smooth/Rough | 0.2008 | N/A | 0.0045 |  1$>$2 | 0.3854 | N/A
Crisp/Jittery |  $<$0.0035 |  S, D $>$ SD | 0.0045 |  1$>$2,3,4 | 0.6953 | N/A
Calm/Anxious |  $<$0.0035 |  SD$>$D$>$S | 0.2987 | N/A | 0.7711 | N/A
Deliberate/Wavering |  $<$0.0035 |  S$>$D, SD |  $<$0.0035 |  1$>$2$>$3$>$4 | 0.4111 | N/A
Relaxed/Tense |  $<$0.0035 |  S$>$D, SD |  $<$0.0035 |  1$>$2,3,4, 2$>$3 | 0.6892 | N/A

S – Spring, D – Damping, SD – Spring+Damping

###  V-B Effect of Force Feedback on Task Performance

Target configuration and needle trajectory for all trials are shown in Figure
6, grouped by the force feedback and color-coded by target error. This figure
shows the traces from all trials receiving each type of force feedback
regardless of the style. The plot visually demonstrates that in general, the
needle trajectory is more confined when applying spring force feedback
compared to the other two types of feedback.

Figure 6: Target layout and resulting needle paths for all subjects. Green
paths indicating smallest error.

(a) Time to Complete the Task (b) Target Positioning Error (c) Needle
Trajectory Straightness

Figure 7: Three metrics including (a) time to complete the task, (b) target
positioning error, (c) and needle trajectory straightness were used to
evaluate the effect of the haptic cues on task performance. For each group
(who received the same type of force feedback), the mean and standard
deviation of each task performance metric is calculated and compared for 4
target locations.

For each trial, task-specific metrics including the time taken to complete the
task, needle trajectory straightness, and target error, were calculated to
evaluate the effects of different types of feedback, as well as the no
feedback condition, on task performance. The mean and standard deviation of
each of these metrics were calculated. Figure 7 (a) shows that the time taken
to reach the target is improved using all three types of force feedback
compared to the no feedback condition. For targets 1 and 2, the group with
spring force feedback showed the least time taken to reach the target, and for
targets 3 and 4, the group with the spring-damping force feedback showed the
least time taken to complete the task.

The target positioning error in all four targets was improved using all three
groups of receiving force feedback compared to the no feedback condition;
however, the group which received the spring+damping force feedback showed the
least error (Fig. 7 (b)). This indicates that applying force feedback
increases the accuracy reaching task regardless of the target.

The straightness of the trajectory traveled by the needle is compared in Fig.
7 (c). This figure indicates that for all four targets, the group which
received the damping feedback showed a straighter needle path compared to the
other two feedback groups and the no feedback condition; however, spring
feedback and spring+damping force feedback caused less straightness in the
needle trajectory compared to the no feedback condition.

In general, the time to complete the task and the target error are improved by
applying at least one of the forces. The straightness of the trajectory
however is only improved compared to absence of feedback only when the damping
feedback is applied. A post hoc statistical analysis was also done for task
performance metric, to determine significant differences in the three types of
force feedback, different targets, and task repetitions. The normality test
was rejected and the Kruskal Wallis test was used to identify significantly
different groups. Effect significance is identified for p-values less than
0.05.

Table III: Statistical analysis summary of the effects of force feedback
types, targets, task repetition on performance metrics

| Force Feedback | Target | Repetition
---|---|---|---
Performance Metric | $p$ | Significance | $p$ | Significance | $p$ | Significance
Target Error | 0.0389 |  D$>$SD |  $<$0.0035 |  2$>$3$>$1$>$4 | .02033 | N/A
Needle Trajectory Straightness | 0.0398 |  SD$>$D |  $<$0.0035 |  1 $>$ 2,3,4, 3 $>$ 4 | 0.914 | N/A
Time taken to complete the task |  $<$0.0035 |  SD, D$>$S |  $<$0.0035 |  2,3,4 $>$ 1 | 0.0498 | N/A

S – Spring, D – Damping, SD – Spring+Damping

For the task performance metrics, the statistical analysis indicates that for
the target positioning error, spring+damping feedback shows significant
difference in reducing the error compared to damping feedback but is not
significant compared to spring feedback; however, it results in a
statistically significant less straighter path traveled by the needle compared
to damping feedback. The spring feedback results in statistically significant
less task completion time compared to damping feedback and spring+damping
feedback (Table III). No statistically significant effect in task performance
metrics, were found due to task repetition; however the target location shows
significant importance in task-specific metrics.

###  V-C Subject Survey

The results from user survey (NASA Task Load Index) indicate that subjects who
received the spring force feedback found the feedback unpleasant and the tasks
more demanding compared to subjects who received the other two types of
feedback. This indicates that a haptic feedback can improve task and user
performance but still be unpleasant to the user.

Figure 8: NASA Task Load Index

##  VI Conclusion

In this study we proposed an automatic training framework in a simulated
environment which detects a poor behavioral performance in the user’s style of
movement in near-real-time and applies force feedback using a haptic device to
help correct the style of movement. We conducted a human subject study to
evaluate the effect of three different types of force feedback: spring force
feedback, damping force feedback, and spring+damping force feedback on six
different behavioral styles. The relation between the quality/style of
movement and user’s skill level was investigated in a previous study [27].
This study builds the groundwork for using haptic as “performance feedback”
(as opposed to reflective or guidance feedback) for improving stylistic
behavior and hence, quality of movement.

The results indicated that “Spring” force feedback resulted in less time to
complete the task, hence faster performance speed. It also helped demonstrate
a more fluid, crisp, calm, and deliberate behavior in the user’s movement.

“Spring+Damping” feedback reduced the target error resulting in a more
accurate performance; however, it resulted in a less straight needle path and
more time to complete the task. It helped demonstrate a more relaxed
performance.

“Damping” force feedback resulted in a straighter line traveled by the needle
in the simulated task towards the target; however, it lead to an increased
target error and a slower speed resulting in more time taken to complete the
task. It helped the user to demonstrate a more crisp performance.

In this study we considered only three type of force feedback related to
position and velocity; however, this might not be sufficient for all styles,
since other kinematic metrics can also be associated with some styles. In
future studies, considering other types of force feedback can help improve the
performance of the styles that were not improved using only a position or
velocity force feedback.

Another area of focus in the future will be to evaluate the effectiveness of
the proposed haptic training method for improving stylistic behavior in long
term. For this purpose, subjects will be trained in three groups, each group
will be receiving a different type of haptic feedback and the their
performance will be monitored and measured over time.

Furthermore this study was carried out in a simulation environment which has
limitations in representing a real surgical task. Future studies will focus on
addressing this issue by implementing the proposed method on the da Vinici
research kit (dvrk) along with a training task.

This study provides the groundwork for continued research on user performance
based feedback for adaptive training.

## Acknowledgment

The authors would like to thank Ziheng Wang for his contributions develop a
bridge between MATLAB and ROS.

The research reported in this publication was supported by the National Center
for Advancing Translational Sciences of the National Institutes of Health
under award Number UL1TR001105. The content is solely the responsibility of
the authors and does not necessarily represent the official views of the NIH.

## References

  * [1] M. Curry, A. Malpani, R. Li, T. Tantillo, A. Jog, R. Blanco, P. K. Ha, J. Califano, R. Kumar, and J. Richmon, “Objective assessment in residency-based training for transoral robotic surgery,” _The Laryngoscope_ , vol. 122, no. 10, pp. 2184–2192, 2012.
  * [2] J. L. Cameron, “William stewart halsted. our surgical heritage.” _Annals of surgery_ , vol. 225, no. 5, p. 445, 1997.
  * [3] R. L. Hoffman, J. Petrosky, M. Eskander, L. Selby, and A. Kulaylat, “Feedback fundamentals in surgical education: Tips for success,” _Bull Am Coll Surg_ , vol. 100, no. 8, pp. 35–39, 2015.
  * [4] A. G. Gallagher, E. M. Ritter, H. Champion, G. Higgins, M. P. Fried, G. Moses, C. D. Smith, and R. M. Satava, “Virtual reality simulation for the operating room: proficiency-based training as a paradigm shift in surgical skills training,” _Annals of Surgery_ , vol. 241, no. 2, pp. 364–372, 2005.
  * [5] I. Badash, K. Burtt, C. A. Solorzano, and J. N. Carey, “Innovations in surgery simulation: a review of past, current and future techniques,” _Annals of translational medicine_ , vol. 4, no. 23, 2016.
  * [6] C. Sewell, D. Morris, N. H. Blevins, S. Dutta, S. Agrawal, F. Barbagli, and K. Salisbury, “Providing metrics and performance feedback in a surgical simulator,” _Computer Aided Surgery_ , vol. 13, no. 2, pp. 63–81, 2008.
  * [7] Z. Pezzementi, D. Ursu, S. Misra, and A. M. Okamura, “Modeling realistic tool-tissue interactions with haptic feedback: A learning-based method,” in _2008 Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems_. IEEE, 2008, pp. 209–215.
  * [8] J. Ko, S.-w. Jang, and Y. S. Kim, “Development of epiduroscopy training simulator using haptic master device,” in _2017 14th International Conference on Ubiquitous Robots and Ambient Intelligence (URAI)_. IEEE, 2017, pp. 542–543.
  * [9] M. Ershad, Z. Koesters, R. Rege, and A. Majewicz, “Meaningful assessment of surgical expertise: Semantic labeling with data and crowds,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI)_. Springer, 2016, pp. 508–515.
  * [10] N. Vaughan, B. Gabrys, and V. N. Dubey, “An overview of self-adaptive technologies within virtual reality training,” _Computer Science Review_ , vol. 22, pp. 65–87, 2016.
  * [11] A. Mariani, E. Pellegrini, N. Enayati, P. Kazanzides, M. Vidotto, and E. De Momi, “Design and evaluation of a performance-based adaptive curriculum for robotic surgical training: a pilot study,” in _2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_. IEEE, 2018, pp. 2162–2165.
  * [12] C. R. Kelley, “What is adaptive training?” _Human Factors_ , vol. 11, no. 6, pp. 547–556, 1969.
  * [13] D. Charles, A. Kerr, M. McNeill, M. McAlister, M. Black, J. Kcklich, A. Moore, and K. Stringer, “Player-centred game design: Player modelling and adaptive digital games,” in _Proceedings of the Digital Games Research Conference_ , vol. 285, 2005, p. 00100.
  * [14] C. Basdogan, S. De, J. Kim, M. Muniyandi, H. Kim, and M. A. Srinivasan, “Haptics in minimally invasive surgical simulation and training,” _IEEE Computer Graphics and Applications_ , vol. 24, no. 2, pp. 56–64, 2004\.
  * [15] N. Enayati, E. De Momi, and G. Ferrigno, “Haptics in robot-assisted surgery: Challenges and benefits,” _IEEE Reviews in Biomedical Engineering_ , vol. 9, pp. 49–65, 2016.
  * [16] C.-H. King, M. O. Culjat, M. L. Franco, C. E. Lewis, E. P. Dutson, W. S. Grundfest, and J. W. Bisley, “Tactile feedback induces reduced grasping force in robot-assisted surgery,” _IEEE Transactions on Haptics_ , vol. 2, no. 2, pp. 103–110, 2009.
  * [17] C. R. Wottawa, B. Genovese, B. N. Nowroozi, S. D. Hart, J. W. Bisley, W. S. Grundfest, and E. P. Dutson, “Evaluating tactile feedback in robotic surgery for potential clinical application using an animal model,” _Surgical Endoscopy_ , vol. 30, no. 8, pp. 3198–3209, 2016.
  * [18] A. Abiri, J. Pensa, A. Tao, J. Ma, Y.-Y. Juo, S. J. Askari, J. Bisley, J. Rosen, E. P. Dutson, and W. S. Grundfest, “Multi-modal haptic feedback for grip force reduction in robotic surgery,” _Scientific reports_ , vol. 9, no. 1, pp. 1–10, 2019.
  * [19] P. Boulanger, G. Wu, W. F. Bischof, and X. D. Yang, “Hapto-audio-visual environments for collaborative training of ophthalmic surgery over optical network,” in _2006 IEEE International Workshop on Haptic Audio Visual Environments and their Applications (HAVE 2006)_ , Nov 2006, pp. 21–26.
  * [20] A. A. Stanley and K. J. Kuchenbecker, “Evaluation of tactile feedback methods for wrist rotation guidance,” _IEEE Transactions on Haptics_ , vol. 5, no. 3, pp. 240–251, Third 2012.
  * [21] X. Yang, W. F. Bischof, and P. Boulanger, “Validating the performance of haptic motor skill training,” in _2008 Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems_ , March 2008, pp. 129–135.
  * [22] T. L. Gibo and D. A. Abbink, “Movement strategy discovery during training via haptic guidance,” _IEEE Transactions on Haptics_ , vol. 9, no. 2, pp. 243–254, 2016.
  * [23] J. Lee and S. Choi, “Effects of haptic guidance and disturbance on motor learning: Potential advantage of haptic disturbance,” in _2010 IEEE Haptics Symposium_. IEEE, 2010, pp. 335–342.
  * [24] W. H. Jantscher, S. Pandey, P. Agarwal, S. H. Richardson, B. R. Lin, M. D. Byrne, and M. K. O’Malley, “Toward improved surgical training: Delivering smoothness feedback using haptic cues,” in _2018 IEEE Haptics Symposium (HAPTICS)_ , March 2018, pp. 241–246.
  * [25] J. C. Gwilliam, M. Mahvash, B. Vagvolgyi, A. Vacharat, D. D. Yuh, and A. M. Okamura, “Effects of haptic and graphical force feedback on teleoperated palpation,” in _2009 IEEE International Conference on Robotics and Automation_. IEEE, 2009, pp. 677–682.
  * [26] W. McMahan, J. Gewirtz, D. Standish, P. Martin, J. A. Kunkel, M. Lilavois, A. Wedmid, D. I. Lee, and K. J. Kuchenbecker, “Tool contact acceleration feedback for telerobotic surgery,” _IEEE Transactions on Haptics_ , vol. 4, no. 3, pp. 210–220, 2011.
  * [27] M. Ershad, R. Rege, and A. M. Fey, “Meaningful assessment of robotic surgical style using the wisdom of crowds,” _International Journal of Computer Assisted Radiology and Surgery_ , pp. 1–12, 2018.
  * [28] ——, “Automatic and near real-time stylistic behavior assessment in robotic surgery,” _International Journal of Computer Assisted Radiology and Surgery_ , vol. 14, no. 4, pp. 635–643, 2019.
  * [29] Y. Gao, S. S. Vedula, C. E. Reiley, N. Ahmidi, B. Varadarajan, H. C. Lin, L. Tao, L. Zappella, B. Béjar, D. D. Yuh _et al._ , “Jhu-isi gesture and skill assessment working set (jigsaws): A surgical activity dataset for human motion modeling,” in _MICCAI Workshop: M2CAI_ , vol. 3, 2014, p. 3.
  * [30] A. Majewicz and A. M. Okamura, “Cartesian and joint space teleoperation for nonholonomic steerable needles,” in _2013 World Haptics Conference (WHC)_. IEEE, 2013, pp. 395–400.
  * [31] M. Ershad, R. Rege, and A. Majewicz, “Surgical skill level assessment using automatic feature extraction methods,” in _Medical Imaging: Image-Guided Procedures, Robotic Interventions, and Modeling_. International Society for Optics and Photonics, 2018.


# Special Glass Structures for First Principles Studies
of Bulk Metallic Glasses

Siya Zhu<EMAIL_ADDRESS>School of Engineering, Brown University,
Providence, RI, 02912, USA  Jan Schroers  Department of Mechanical Engineering
and Materials Science, Yale University, New Haven, CT, 06511, USA  Stefano
Curtarolo  Department of Mechanical Engineering and Materials Science and
Center for Autonomous Materials Design, Duke University, Durham, NC, 27708,
USA  Hagen Eckert  Department of Mechanical Engineering and Materials Science
and Center for Autonomous Materials Design, Duke University, Durham, NC,
27708, USA  Axel van de Walle  School of Engineering, Brown University,
Providence, RI 02912, USA

###### Abstract

The atomic-level structure of bulk metallic glasses is a key determinant of
their properties. An accurate representation of amorphous systems in
computational studies has traditionally required large supercells that are
unfortunately computationally demanding to handle using the most accurate ab
initio calculations. To address this, we propose to specifically design small-
cell structures that best reproduce the local geometric descriptors (e.g.,
pairwise distances or bond angle distributions) of a large-cell simulation. We
rely on molecular dynamics (MD) driven by empirical potentials to generate the
target descriptors, while we use reverse Monte Carlo (RMC) methods to optimize
the small-cell structure. The latter can then be used to determine mechanical
and electronic properties using more accurate electronic structure
calculations. The method is implemented in the Metallic Amorphous Structures
Toolkit (MAST) software package.

††preprint: APS/123-QED

##  I Introduction

For decades, bulk metallic glasses (BMGs) have been investigated for their
excellent strength and fracture toughness [1, 2, 3, 4, 5, 6, 7]. Traditional
amorphous alloys require high quenching rates to avoid crystallization, which
restricts the size of powders and films [8]. In contrast, the critical cooling
rates of BMGs are much lower, allowing the size to be greater than
$1\text{\,}\mathrm{mm}$ in the smallest dimension, and promising the
application as structural materials. [9, 10, 11]

The atomic structures of BMGs are of prime importance for predicting and
explaining their macroscopic properties [12, 13, 14]. Experimental
characterizations with X-ray diffraction (XRD), small-angle neutron
scattering(SANS), and transmission electron microscopy (TEM) help us
understand some structural properties [15], yet they are insufficient to fully
reconstruct the atomic packing of BMGs. Molecular dynamics (MD) is often used
in the investigation of atomic structures of BMGs [16]. Traditional embedded
atom method (EAM) MD can efficiently simulate the cooling process over a time
scale of $100\text{\,}\mathrm{ns}$ for systems of over 10,000 atoms [13].
While this is typically sufficient to reproduce the general atomic structure
of BMG, EAM-MD exhibits two main shortcomings: First, due to the limited
accuracy of EAM, the finer details of the structure, such as specific bond
lengths or bond angles, may not be accurately reproduced. As a result, the
accuracy of some physical properties, particularly viscosity and plastic
mechanical properties, as well as normal mode frequencies, or even material
density, may suffer. Second, the structures obtained from EAM-MD cannot, due
to their sizes, be directly used for further, more accurate analysis using
first-principles calculations. This, therefore, precludes the calculations of
electronic properties as well as more accurate calculations of mechanical
properties.

To address this, it has been proposed [12] to use first-principles methods for
both the simulation of the quenching process and the calculation of the
properties of interest. While first-principles methods guarantee a higher
accuracy in calculating the energies and forces, the system size is typically
limited to $<200$ atoms and the time scale to picoseconds. As a result, the
large statistical fluctuation and periodic boundary conditions in such a small
system size as well as the high quenching rate due to the short simulation
time can significantly affect the structure of the simulated system.

In an effort to leverage the advantages of both EAM-MD and fully first-
principles methods, we propose to combine these methods with the help of
Reverse Monte Carlo (RMC). The idea of RMC has often been used to infer the
atomic structures of BMG and other disordered systems [17, 18, 19] from
indirect experimental characterization techniques. For instance, X-ray
diffraction, neutron diffraction, extended X-ray absorption fine structure
(EXAFS), etc., can provide structure factors, radial distribution functions,
or atomic coordinations. RMC modeling, a variation of the standard Metropolis
Monte Carlo method, can then be used to generate structural model (a cell or
an ensemble of atoms) by minimizing an error function (also called loss
function) quantifying the consistency of the generated structure with the
experimental data, while subjecting to a set of constraints like
concentrations or cell shape. Various programs implementing this algorithm,
such as fullrmc [20], RMCProfile [21], and RMC++ [22] are widely used to study
the atomic modeling of liquids, glasses, polymers and other amorphous
materials.

Here we use RMC in differently. Instead of using experimental data as input,
we perform RMC based on structures from EAM-MD results in order to build
representative atomic structures that are small enough for DFT calculations,
while maintaining the statistical geometrical properties of the large
structure from EAM-MD.

We have developed a C++ software package — Metallic Amorphous Structures
Toolkit (MAST), as a part of Alloy Theoretic Automated Toolkit (ATAT) [23,
24]. MAST includes RMC codes for generating random structures suitable for MD
and MC run, as well as scripts to automate the entire process, including MD,
MC and ab initio calculations. It also contains scripts for calculating
mechanical properties such as elastic constants.

We first generate a typical BMG atomic configuration (with over 10,000 atoms)
using EAM-MD. Then, we use the RMC algorithm to construct a small periodic
system of 32 atoms consistent with the MD result in terms of statistical
properties such as pair distribution functions and bond angles. We give the
name Special Glass Structure (SGS) to such atomic structure, paraphrasing the
concept of Special Quasirandom Structures (SQS) [25, 26, 27] for crystalline
structures. Once an optimal small-cell representation has been determined, ab
initio calculations for structural optimization and mechanical properties
computations are performed.

In this work, we study, as an example, Zr-Cu and Zr-Cu-Al based BMGs over a
range of compositions, since these systems represent the transition metal BMG
former and are the basis of various important BMGs with attractive properties
[28, 29]. Atomic structures and mechanical properties are calculated and
compared with experimental data and other theoretical predictions.

##  II Methods

###  II.1 Error Functions

To measure the structural differences between a large target structure with
over 10,000 atoms and an SGS with 32 atoms, we devise an error function
containing various kinds of commonly used structural information.

This information includes the pair distribution functions (also known as the
partial radial distribution functions, PDF), defined as

$g_{\alpha\beta}\left(r\right)=\frac{1}{4\pi r^{2}\rho
c_{\beta}}\frac{\text{d}n_{\alpha\beta}\left(r\right)}{\text{d}r}$ (1)

where $\rho$ is the atomic number density, $n_{\alpha\beta}\left(r\right)$ is
the average number of atoms type $\beta$ within a distance of $r$ from an atom
of type $\alpha$, and $c_{\beta}$ is the concentration of atom type $\beta$ in
the structure. To make it symmetric to $\alpha$ and $\beta$, it can also be
written as

$g_{\alpha\beta}\left(r\right)=\frac{N}{4\pi r^{2}\rho
N_{\alpha}N_{\beta}}\frac{\text{d}N_{\alpha\beta}\left(r\right)}{\text{d}r}=\frac{V}{4\pi
r^{2}N_{\alpha}N_{\beta}}\frac{\text{d}N_{\alpha\beta}\left(r\right)}{\text{d}r}$
(2)

where $V$ is the volume, $N$ is the total number of atoms,
$\text{d}N_{\alpha\beta}\left(r\right)$ is the total number of $\alpha$ and
$\beta$ atoms pairs within distance $r$ to $r+\text{d}r$, $N_{\alpha}$ and
$N_{\beta}$ are the total numbers of atoms type $\alpha$ and $\beta$,
respectively. A traditional way to compare the PDFs is to obtain the structure
factors with a Fourier transformation on the PDFs, sum over all partial
structure factors for the total structure factors, and perform a $\chi^{2}$
test with the one from experiments like XRD. Since we only compare between
both atomic structures instead of experimental data, we can compute the
difference between two PDFs directly. In order to get a smooth PDF without
sharp peaks for small structures or noises for large systems, we perform
kernel smoothing [30]:

$\bar{g}_{\alpha\beta}\left(r\right)\equiv\frac{1}{s}\int_{-\infty}^{\infty}K\left(\frac{r^{\prime}-r}{s}\right)g_{\alpha\beta}\left(r^{\prime}\right)\text{d}r^{\prime}$
(3)

where $K$ is a symmetric kernel function integrating to one and $s$ is a user-
specified smoothing parameter.

Then, the error function for the PDFs is defined as

$E_{\text{PDF}}=\sum_{\alpha,\beta}\int_{0}^{r_{\text{MAX}}}\frac{1}{r^{2}}\left[\bar{g}^{\text{o}}_{\alpha\beta}\left(r\right)-\bar{g}_{\alpha\beta}\left(r\right)\right]^{2}\text{d}r$
(4)

where $r_{\text{MAX}}$ is a cutoff distance (usually set as the distance of
the 2nd or 3rd nearest neighbors),
$\bar{g}^{\text{o}}_{\alpha\beta}\left(r\right)$ is the PDF of a target
structure. We include a factor $1/r^{2}$ to give higher weights to shorter
distances, as the chemical environment closer to an atom is more important
than farther away. Notice that, since the PDF is symmetric in $\alpha$ and
$\beta$, we do not double count the $g_{\beta\alpha}$ term in the sum. For
example, for a Zr-Cu system, we only sum over three terms of partial PDFs: Zr-
Zr, Zr-Cu, and Cu-Cu, but not Cu-Zr. For a ternary system like Zr-Cu-Al, we
sum over 6 terms in $E_{\text{PDF}}$.

For bond angles, we define a statistical quantity for the angles between two
bonds connecting an atom $\beta$ to two nearest neighbors site occupied by
species $\alpha$ and $\gamma$:

$h_{\alpha\beta\gamma}\left(\theta\right)=\frac{VN}{N_{\alpha}N_{\beta}N_{\gamma}}\frac{\text{d}N_{\alpha\beta\gamma}\left(\theta\right)}{\text{d}\theta}$
(5)

where $\text{d}N_{\alpha\beta\gamma}\left(\theta\right)$ is the number of
angles forming with bond $\alpha$-$\beta$ and $\beta$-$\gamma$, within the
range of $\theta$ to $\theta+\text{d}\theta$. The bond length should be less
than a cutoff length, usually set between the distance of the 1st and 2nd
nearest neighbours. We also kernel smooth the function
$h_{\alpha\beta\gamma}\left(\theta\right)$. The error function for angles is
defined as

$E_{\text{angles}}=\sum_{\alpha,\beta,\gamma}\int_{0}^{\pi}\left[\bar{h}^{\text{o}}_{\alpha\beta\gamma}\left(\theta\right)-\bar{h}_{\alpha\beta\gamma}\left(\theta\right)\right]^{2}\text{d}\theta$
(6)

Notice that $\alpha$ and $\gamma$ are equivalent and $\beta$ is not,
therefore, there are 6 terms in the sum for a binary system and 18 terms for a
ternary system.

For a total error function, we weigh the errors from each aspect:

$E=w_{\text{PDF}}E_{\text{PDF}}+w_{\text{angles}}E_{\text{angles}}$ (7)

A basic principle for choosing $w_{i}$ is to make every term $w_{i}E_{i}$ of
the same order of magnitude, so that they are all taken into account in the
Monte Carlo simulations.

###  II.2 RMC process

With an initial guess, which can be either a MD simulated structure, a
previous MC result, or a random structure, we do an RMC process:

  * i.

Randomly pick an operation from: (1) move an atom by a random amount; (2) swap
positions of two atoms; (3) reshape the supercell by a random amount. The
probabilities of each operation are chosen for better convergence. A suggested
improvement is to set the initial structure as a cubic cell with volume based
on the density of target structure, so that we do not need to reshape the cell
(which is a default option in our code). Reshaping cells is sometimes
detrimental to the convergence of the method and is not recommended;

  * ii.

Calculate the error function $E_{\text{new}}$ for the new structure with the
target structure;

  * iii.

Compare with the error function $E$ of the old structure. If
$E_{\text{new}}<E$, we always accept the operation. Otherwise, we accept the
operation with a probability
$\text{exp}\left[\left(E-E_{\text{new}}\right)/{T}\right]$, where $T$ is a
“temperature” factor greater than 0. The larger $T$ is, the more likely to get
over barriers and find the global minimum, while a smaller $T$ speeds up the
convergence to a local minimum.

  * iv.

Go back to the first step, until we reach the $E$ we set, or the maximum
number of Monte Carlo steps.

###  II.3 Generation of atomic modeling

As a specific example, to generate a 32-atom SGS with concentration Zr50Cu50,
we follow this process:

  * i.

Use MAST to generate a random structure with 8000 Zr atoms and 8000 Cu atoms;

  * ii.

Use LAMMPS [31, 32] to run an EAM-MD simulation. The EAM potentials in this
work are from the group of Sheng [13]. We first do a NVT and a NPT run at
$2000\text{\,}\mathrm{K}$, and quenching to $300\text{\,}\mathrm{K}$ at
$10^{12}$ K/s with NPT ensemble. Then a NPT at $300\text{\,}\mathrm{K}$ is
performed for a stable structure;

  * iii.

Use MAST for an RMC process from a cell with 16 Zr atoms and 16 Cu atoms
randomly arranged to generate an SGS;

  * iv.

Use LAMMPS to do a structural optimization for the SGS from step(iii) with
command “minimize”, then a NPT MD at $300\text{\,}\mathrm{K}$. This is because
RMC may sometimes give unreasonable structures (with two atoms too close to
each other). An optimization and MD run with EAM help us avoid this situation;

  * v.

Another RMC run on the result of step (iv). More loops of MD and RMC (step
(iv) and step (v))can be performed for a better convergence. In this work, we
do not perform extra loops of MD and RMC since it already gives very good
fitting.

  * vi.

Use Vienna ab initio Simulation Package (VASP) [33, 34, 35] to do an ab initio
optimization and get the final structure for Zr16Cu16. We use the Perdew-
Burke-Ernzerhof (PBE) exchange and correlation functional at level of the
generalized gradient approximation(GGA) [36] . A cutoff energy of 600eV is set
to eliminate the effect of Pulay stress in structural optimization with
ISIF=3. A $7\times 7\times 7$ K-points mesh is set.

###  II.4 Calculation of Elastic Moduli

As an isotropic material, a BMG only has two independent elastic constants:
$C_{11}$ and $C_{12}$. To calculate them, we apply strains on the relaxed
structure, optimize with VASP (ISIF = 2), calculate the energy and fit the
strain-energy relationship. A planewave energy cutoff of 300eV is used in the
mechanical properties calculations and other settings are the same as
previously stated. We apply uniaxial strain from $-0.03$ to $+0.03$ on $x$,
$y$ and $z$ directions (the strain-energy curves on each direction are almost
the same, which will be shown in the next section, as an evidence of isotropy
of our structure). $C_{11}$ can be derived from

$C_{11}=\frac{2}{V_{0}}\frac{\text{d}^{2}E_{\text{uniaxial}}}{\text{d}\epsilon^{2}}$
(8)

Then we add a bulk strain $\epsilon$ from $-0.03$ to $+0.03$ on each direction
and plot the strain-energy curve. Then $C_{12}$ can be calculated as

$C_{12}=\frac{1}{3V_{0}}\frac{\text{d}^{2}E_{\text{triaxial}}}{\text{d}\epsilon^{2}}-\frac{C_{11}}{2}$
(9)

The bulk modulus is

$K=\frac{C_{11}}{3}+\frac{2C_{12}}{3}=\frac{2}{9V_{0}}\frac{\text{d}^{2}E_{\text{triaxial}}}{\text{d}\epsilon^{2}}$
(10)

and shear modulus is

$G=C_{44}=\frac{C_{11}-C_{12}}{2}$ (11)

##  III Results

In this section, we first use an example of Zr50Cu50 BMG to illustrate our
results. First, we perform a MD simulation to generate a structure with 8000
Zr and 8000 Cu atoms, as shown in Figure 1a. Pair distribution functions and
angle distribution functions are calculated and shown in Figure 1b and 1c.

Figure 1: (a) An 16000-atom atomic structure of Zr50Cu50 constructed with MD.
Green atoms are Zr and blue atoms are Cu. (b) Pair distribution functions of
the structure.(c) Angle distribution functions of the structure, with cutoff
distance 4A.

Following the steps in Section II.3, we generate the atomic structures with
steps (iii) to (v), calculate pair and angle distribution functions as shown
in Figure 2. The RMC steps sometimes generate structures with nonzero PDF at
very short distances, thus implying that some atoms are nonphysically close to
each other. This can be easily remedied either by performing relaxations and
MD of the atoms using the EAM model, as illustrated in Figure 2d-2f. Another
RMC process reducing the error function gives an SGS Figure 2g with good
fitting, as shown in Figure 2h and 2i. The SGS is now physically reasonable
and successfully represents the structural information of the structure in
Figure 1a.

Figure 2: Atomic structures, PDFs and angle distributions at different steps
in Section II.3. (a) SGS after step (iii); (b) PDFs for structure(a); (c)
Angle distributions for structure(a); (d) SGS after step (iv); (e) PDFs for
structure(b); (f) Angle distributions for structure(b); (g) SGS after step
(v); (h) PDFs for structure(c); (i) Angle distributions for structure(c).

After the structural optimization with VASP, we apply uniaxial strain and
triaxial strain to our SGS, optimize the structure under strain, and calculate
the energy, as shown in Figure 3. Two things need to be noted here: (i) we add
the strain from $-0.03$ to $+0.03$, which might exceed the range of elastic
deformation of BMG. However, in this case, the SGS is too small to perform any
plastic deformation or fracture, and the energies still fit well with large
strain. Therefore we keep the data of large strain up to 0.03 in this example.
The strain range can be set manually in MAST; (ii) all ab initio calculations
of energy here are at the temperature of $0\text{\,}\mathrm{K}$. That would
bring about some error in the elastic modulus compared to experimental data at
room temperature. Nevertheless, the differences between mechanical properties
at $0\text{\,}\mathrm{K}$ and $300\text{\,}\mathrm{K}$ are insignificant.
Figure 3b shows the energy under uniaxial strain applied to $x$, $y$ and $z$
directions, respectively. The strain-energy curves of three directions almost
overlap, manifesting that our SGS is “isotropic” in some sense, which provides
additional evidence of its validity.

Figure 3: (a) Strain-energy curve of Zr50Cu50 SGS with uniaxial and triaxial
strain; (b) Strain-energy curve of Zr50Cu50 SGS with uniaxial on $x$, $y$ and
$z$ directions, respectively

In Table 1, we perform a robustness test of our method on 32 atoms structure
of Zr-Cu BMG at concentration $x=0.5$. With four different initial structures
in step (iii) in Section II.3, we get four different SGSs and calculate their
mechanical properties. The results show our Monte Carlo process converges and
is robust and reliable.

Table 1: Robustness test of MAST on Zr50Cu50 (32 atoms) | C11(GPa) | C12(GPa) | K(GPa) | G(GPa) | E (GPa) | $\nu$
---|---|---|---|---|---|---
1 | 141.34 | 89.77 | 106.96 | 25.78 | 71.60 | 0.388
2 | 142.11 | 89.23 | 106.86 | 26.44 | 73.27 | 0.386
3 | 142.30 | 92.63 | 109.19 | 24.83 | 69.25 | 0.394
4 | 142.13 | 92.88 | 109.30 | 24.63 | 68.72 | 0.395

In Table 2, we compare the mechanical properties of Zr50Cu50 from experiments,
our method, and some other methods to show the advantage of our method. For
the SGS in Figure 2(c) with 32 atoms, the differences in bulk modulus, shear
modulus, Young’s modulus and Poisson’s ratios are less than 15% (labeled MAST
(32 atoms) in Table 2). For comparison, we apply strain on the structure with
16,000 atoms in Figure 1a, run MD with EAM potentials for energies, and
calculate mechanical properties(labeled as EAM-MD (16,000 atoms) in Table 2).
The results from our process MAST agree with experimental data a little bit
better — both methods overestimate the bulk modulus K, while MAST
underestimates shear and Young’s modulus and EAM-MD overestimates them,
compared with experiments. To show the beneficial contribution of our Monte
Carlo process in MAST, we also show computations for the glass structure with
32 atoms, generated with EAM-MD only. We calculate the mechanical properties
with EAM-MD itself(labeled as EAM-MD (32 atoms) in Table 2) and ab initio
(labeled as EAM-MD + ab initio (32 atoms) in Table 2). Both of them show poor
agreement with experiments. An SGS with only 12 atoms (Zr6Cu6) is also
generated with MAST(labeled as MAST(12 atoms) in Table 2). The mechanical
properties of it are close to MAST 32 atoms SGS, and are even better than EAM-
MD structures with 32 atoms, which shows the robustness and accuracy of our
method.

Table 2: Mechanical properties of Zr50Cu50 by different methods | K(GPa) | G(GPa) | E (GPa) | $\nu$ | $R^{2}$ | $R_{\text{adj}}^{2}$
---|---|---|---|---|---|---
Experiment [37] | 101.2 | 31.3 | 84.0 | 0.35 |  |
MAST(32 atoms, average) | 108.08 | 25.42 | 70.7 | 0.391 |  $>$0.998 |  $>$0.998
MAST(12 atoms) | 111.9 | 24.5 | 68.6 | 0.398 | 0.987 | 0.986
EAM-MD (16,000 atoms) | 112.3 | 38.9 | 104.5 | 0.344 | 0.999 | 0.998
EAM-MD + ab initio (32 atoms) | 107.3 | 15.9 | 45.3 | 0.430 | 0.991 | 0.990
EAM-MD (32 atoms) | 110.6 | 40.7 | 108.7 | 0.336 | 0.926 | 0.911

Similarly, our method can be applied to BMG systems with different
concentrations, or even more than two elements. We show some calculation
results for ZrCu and ZrCuAl BMG in Table 3, together with some experimental
data as reference. For most of the concentrations, our results agree well with
experiments on the bulk moduli, but the levels of agreement decrease for shear
moduli and Young’s moduli. This originates, in large part, from the fact that
ab initio bulk modulus calculations exhibit a higher signal-to-noise ratio,
because a triaxial strain induces a larger change in energy than a uniaxial
strain and only the bulk modulus depends solely on triaxial strain data (c.f.
Eq. 10). As seen in Figure 3, the overall energy changes associated with
triaxial strains are about 10 times bigger than those associated with uniaxial
strains. In spite of this, the SGS approach still provides sufficient accuracy
to investigate the compositional-dependence of elastic properties and, in this
case, provides corroboration of the relatively small sensitivity to
composition found in experiments.

Table 3: Mechanical properties of ZrCu and ZrCuAl BMGs with different
concentrations

| K(GPa,cal) | K(GPa,exp) | G(GPa,cal) | G(GPa,exp) | E (GPa,cal) | E(GPa,exp) |  $\nu$(cal) |  $\nu$(exp)
---|---|---|---|---|---|---|---|---
Zr36Cu64 [28, 38] |  | 104.3 |  | 34 |  | 92.3 |  | 0.352
Zr43.7Cu56.3 | 112.9 |  | 25.3 |  | 70.6 |  | 0.396 |
Zr46.9Cu53.1 | 112.8 |  | 20.2 |  | 57.1 |  | 0.416 |
Zr50Cu${}_{5}0$ [37] | 108.8 | 101.2 | 25.4 | 31.3 | 70.7 | 84.0 | 0.391 | 0.35
Zr53.1Cu46.9 | 106.9 |  | 18.5 |  | 52.6 |  | 0.418 |
Zr56.3Cu43.7 | 107.8 |  | 21.1 |  | 59.5 |  | 0.408 |
Zr54Cu46 [39, 38] |  | 128.5 |  | 30.0 |  | 83.5 |  | 0.391
Zr47.5Cu47.5Al5 [37] |  | 113.7 |  | 33.0 |  | 90.1 |  | 0.365
Zr47.5Cu47.5Al5 [40] |  | - |  | - |  |  $\sim$60† |  | -
Zr47.5Cu47.5Al5 [41] |  | 117$\pm$1 |  | 32$\pm$1 |  | 88$\pm$1 |  | -
Zr47.5Cu47.5Al5 [42] |  | - |  | - |  | 87$\pm$5 |  | -
Zr47Cu47Al6 [42] |  | - |  | - |  | 84$\pm$5 |  | -
Zr46.9Cu46.9Al6.2 | 108.2 |  | 22.3 |  | 62.6 |  | 0.40 |
Zr45Cu45Al10 [43] |  | - |  | - |  |  $\sim$60† |  | -
Zr45Cu45Al10 [40] |  | - |  | - |  |  $\sim$80† |  | -

  * •

† estimated from reported graph;

The usefulness of SGS is not limited to mechanical properties. For instance,
they make it possible to analyze detailed electronic structure of BMGs, which
would be impossible with EAM alone and intractable with large supercells using
ab initio methods. In Figure 4a we show the calculated DOS for different
compositions of Zr-Cu glasses. This type of analysis enables us, for instance,
to assess whether this BMG is well described by a rigid-band approximation. In
the case of Zr-Cu, we can actually determine that shifts in the DOS as a
function of composition, relative to the Fermi level, are not uniform over all
energies. For instance, around 3eV, there is a larger shift between Zr14Cu18
and Zr15Cu17 composition than between Zr15Cu17 and Zr18Cu14 composition, while
to opposite is true around $-2$ eV.

These effects would be difficult to isolate using other methods. For instance,
in Figure 4b, we make a comparison between the DOS obtained via an SGS from
MAST and the one obtained via EAM-MD only, showing that the electronic
properties would likely be incorrect without the Monte Carlo optimization
process. The two DOS differ by shifts of up to 0.2 eV on the energy scale,
which is large enough to mask the effects identified in Figure 4a.

Figure 4: (a) DOS for ZrCu SGS generated by MAST at different concentrations;
(b) DOS for Zr16Cu16 SGS by MAST and EAM.

##  IV Conclusion

In this work, we develop a method to generate small-cell approximations to
amorphous structures suitable for efficient ab initio calculations. These so-
called Special Glass Structures (SGS) are constructed by matching geometric
descriptors from large-cell empirical potential MD simulations using reverse
Monte Carlo. The method is benchmarked by comparing the predicted mechanical
properties of common Bulk Metallic Glasses to corresponding experimental
results.

More generally, our method can be used for the study of any amorphous system,
and is applicable to the calculation of properties besides the elastic
constants, such as thermodynamics properties like enthalpy and vibrational
entropy with a higher accuracy than classical potential-based MD. Larger SGS
($\sim 100$ atoms) should also apply to the study of kinetic processes, such
as viscous flow or plastic deformation processes via shear band evolution in
BMGs [44]. While, we leave the demonstration of these capabilities for
subsequent work, it should be clear that the SGS method opens the way to the
systematic study of BMG properties via high-throughput ab initio methods over
a broad range of chemistries [24, 45].

###### Acknowledgements.

This work is supported by Office of Naval Research grant N00014-20-1-2225.
Computational resources were provided by (i) the Center for Computation and
Visualization at Brown University, (ii) the Extreme Science and Engineering
Discovery Environment (XSEDE) through allocation TG-DMR050013N, which is
supported by National Science Foundation Grant No. ACI-1548562 and (iii) the
Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support
(ACCESS) program through allocation DMR010001, which is supported by National
Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296.

## References

  * [1] Christopher A Schuh, Todd C Hufnagel, and Upadrasta Ramamurty.  Mechanical behavior of amorphous alloys.  Acta Materialia, 55(12):4067–4109, 2007.
  * [2] A Reza Yavari, JJ Lewandowski, and J Eckert.  Mechanical properties of bulk metallic glasses.  Mrs Bulletin, 32(8):635–638, 2007.
  * [3] J Eckert, J Das, S Pauly, and C Duhamel.  Mechanical properties of bulk metallic glasses and composites.  Journal of materials research, 22(2):285–301, 2007.
  * [4] A Lindsay Greer and Evan Ma.  Bulk metallic glasses: at the cutting edge of metals research.  MRS bulletin, 32(8):611–619, 2007.
  * [5] MF Ashby and A Lindsay Greer.  Metallic glasses as structural materials.  Scripta Materialia, 54(3):321–326, 2006.
  * [6] Ling Shao, Jittisa Ketkaew, Pan Gong, Shaofan Zhao, Sungwoo Sohn, Punnathat Bordeenithikasem, Amit Datye, Rodrigo Miguel Ojeda Mota, Naijia Liu, Sebastian Alexander Kube, et al.  Effect of chemical composition on the fracture toughness of bulk metallic glasses.  Materialia, 12:100828, 2020.
  * [7] Henry J Neilson, Alex S Petersen, Andrew M Cheung, S Joseph Poon, Gary J Shiflet, Mike Widom, and John J Lewandowski.  Weibull modulus of hardness, bend strength, and tensile strength of ni- ta- co- x metallic glass ribbons.  Materials Science and Engineering: A, 634:176–182, 2015.
  * [8] Jan Schroers.  Processing of bulk metallic glass.  Advanced materials, 22(14):1566–1597, 2010.
  * [9] Akihisa Inoue.  Stabilization of metallic supercooled liquid and bulk amorphous alloys.  Acta materialia, 48(1):279–306, 2000.
  * [10] XH Lin and WL Johnson.  Formation of ti–zr–cu–ni bulk metallic glasses.  Journal of Applied Physics, 78(11):6514–6519, 1995.
  * [11] Donghua Xu, Gang Duan, and William L Johnson.  Unusual glass-forming ability of bulk amorphous alloys based on ordinary metal copper.  Physical review letters, 92(24):245504, 2004.
  * [12] H.W. Sheng, W.K. Luo, F.M. Alamgir, J.M. Bai, and E Ma.  Atomic packing and short-to-medium-range order in metallic glasses.  Nature, 439(7075):419–425, 2006.
  * [13] Y.Q. Cheng, E Ma, and H.W. Sheng.  Atomic level structure in multicomponent bulk metallic glass.  Physical review letters, 102(24):245501, 2009.
  * [14] W.K. Luo, H.W. Sheng, F.M. Alamgir, J.M. Bai, J.H. He, and E Ma.  Icosahedral short-range order in amorphous alloys.  Physical Review Letters, 92(14):145502, 2004.
  * [15] YQ Cheng and E Ma.  Atomic-level structure and structure–property relationship in metallic glasses.  Progress in materials science, 56(4):379–473, 2011.
  * [16] Gang Duan, Donghua Xu, Qing Zhang, Guoyun Zhang, Tahir Cagin, William L Johnson, and William A Goddard III.  Molecular dynamics study of the binary cu 46 zr 54 metallic glass motivated by experiments: Glass formation and atomic-level structure.  Physical Review B, 71(22):224208, 2005.
  * [17] László Pusztai and Erzsébet Sváb.  Structure study of ni62nb38 metallic glass using reverse monte carlo simulation.  Journal of non-crystalline solids, 156:973–977, 1993.
  * [18] P Lamparter.  Reverse monte carlo simulation of amorphous ni80p20 and ni81b19.  Physica Scripta, 1995(T57):72, 1995.
  * [19] Robert L McGreevy.  Reverse monte carlo modelling.  Journal of Physics: Condensed Matter, 13(46):R877, 2001.
  * [20] Bachir Aoun.  Fullrmc, a rigid body reverse monte carlo modeling package enabled with machine learning and artificial intelligence.  Journal of computational chemistry, 37(12):1102–1111, 2016.
  * [21] Matthew G Tucker, Martin T Dove, and David A Keen.  Simultaneous analysis of changes in long-range and short-range structural order at the displacive phase transition in quartz.  Journal of Physics: Condensed Matter, 12(48):L723, 2000.
  * [22] Guillaume Evrard and László Pusztai.  Reverse monte carlo modelling of the structure of disordered materials with rmc++: a new implementation of the algorithm in c++.  Journal of Physics: Condensed Matter, 17(5):S1, 2005.
  * [23] Axel van de Walle, M Asta and G Ceder.  The alloy theoretic automated toolkit: A user guide.  Calphad, 26(4):539–553, 2002.
  * [24] Axel van de Walle, Ruoshi Sun, Qi-Jun Hong, and Sara Kadkhodaei.  Software tools for high-throughput calphad from first-principles data.  Calphad, 58:70–81, 2017.
  * [25] Alex Zunger, S-H Wei, LG Ferreira, and James E Bernard.  Special quasirandom structures.  Physical review letters, 65(3):353, 1990.
  * [26] A Van de Walle, P Tiwary, M De Jong, DL Olmsted, M Asta, A Dick, D Shin, Yi Wang, L-Q Chen, and Z-K Liu.  Efficient stochastic generation of special quasirandom structures.  Calphad, 42:13–18, 2013.
  * [27] James E Saal and Chris Wolverton.  Thermodynamic stability of co–al–w l12 $\gamma$.  Acta materialia, 61(7):2330–2338, 2013.
  * [28] Donghua Xu, Boonrat Lohwongwatana, Gang Duan, William L Johnson, and Carol Garland.  Bulk metallic glass formation in binary cu-rich alloy series–cu100- xzrx (x= 34, 36, 38.2, 40 at.%) and mechanical properties of bulk cu64zr36 glass.  Acta Materialia, 52(9):2621–2624, 2004.
  * [29] D Wang, Y Li, BB Sun, ML Sui, K Lu, and E Ma.  Bulk metallic glass formation in the binary cu–zr system.  Applied Physics Letters, 84(20):4029–4031, 2004.
  * [30] MP Wand and MC Jones.  Monographs on statistics and applied probability.  Kernel smoothing, 1995.
  * [31] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of computational physics, 117(1):1–19, 1995.
  * [32] Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan S Bolintineanu, W Michael Brown, Paul S Crozier, Pieter J in’t Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, et al.  Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Computer Physics Communications, 271:108171, 2022.
  * [33] Georg Kresse and Jürgen Hafner.  Ab initio molecular dynamics for liquid metals.  Physical review B, 47(1):558, 1993.
  * [34] Georg Kresse and Jürgen Furthmüller.  Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set.  Computational materials science, 6(1):15–50, 1996.
  * [35] Georg Kresse and Jürgen Furthmüller.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  Physical review B, 54(16):11169, 1996.
  * [36] John P Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Physical review letters, 77(18):3865, 1996.
  * [37] GJ Fan, M Freels, H Choo, PK Liaw, JJZ Li, Won-Kyu Rhim, WL Johnson, P Yu, and WH Wang.  Thermophysical and elastic properties of cu 50 zr 50 and (cu 50 zr 50) 95 al 5 bulk-metallic-glass-forming alloys.  Applied physics letters, 89(24):241917, 2006.
  * [38] WL Johnson and Konrad Samwer.  A universal criterion for plastic yielding of metallic glasses with a (t/t g) 2/3 temperature dependence.  Physical review letters, 95(19):195501, 2005.
  * [39] YQ Cheng, AJ Cao, and E Ma.  Correlation between the elastic modulus and the intrinsic plastic behavior of metallic glasses: The roles of atomic configuration and alloy composition.  Acta Materialia, 57(11):3253–3267, 2009.
  * [40] G Kumar, T Ohkubo, T Mukai, and K Hono.  Plasticity and microstructure of zr–cu–al bulk metallic glasses.  Scripta materialia, 57(2):173–176, 2007.
  * [41] NS Barekar, S Pauly, RB Kumar, U Kühn, BK Dhindaw, and J Eckert.  Structure–property relations in bulk metallic cu–zr–al alloys.  Materials Science and Engineering: A, 527(21-22):5867–5872, 2010\.
  * [42] S Pauly, G Liu, S Gorantla, G Wang, U Kühn, DH Kim, and J Eckert.  Criteria for tensile plasticity in cu–zr–al bulk metallic glasses.  Acta materialia, 58(14):4883–4890, 2010.
  * [43] Peter Thurnheer, Robert Maass, Stefan Pogatscher, and Jörg F Löffler.  Compositional dependence of shear-band dynamics in the zr–cu–al bulk metallic glass system.  Applied Physics Letters, 104(10):101910, 2014.
  * [44] S Xie and Easo P George.  Hardness and shear band evolution in bulk metallic glasses after plastic deformation and annealing.  Acta Materialia, 56(18):5202–5213, 2008.
  * [45] Stefano Curtarolo, Gus L. W. Hart, Natalio Nardelli, Marco Buongiorno an d Mingo, Stefano Sanvito, and Ohad Levy.  The high-throughput highway to computational materials design.  Nat. Mater., 12:191, 2013.