# LAMPAT: Low-Rank Adaption for Multilingual Paraphrasing Using Adversarial
Training

Khoi M. Le1, 2, Trinh Pham2, Tho Quan2, Anh Tuan Luu3 Corresponding author.

###### Abstract

Paraphrases are texts that convey the same meaning while using different words
or sentence structures. It can be used as an automatic data augmentation tool
for many Natural Language Processing tasks, especially when dealing with low-
resource languages, where data shortage is a significant problem. To generate
a paraphrase in multilingual settings, previous studies have leveraged the
knowledge from the machine translation field, i.e., forming a paraphrase
through zero-shot machine translation in the same language. Despite good
performance on human evaluation, those methods still require parallel
translation datasets, thus making them inapplicable to languages that do not
have parallel corpora. To mitigate that problem, we proposed the first
unsupervised multilingual paraphrasing model, LAMPAT (Low-rank Adaptation for
Multilingual Paraphrasing using Adversarial Training), by which monolingual
dataset is sufficient enough to generate a human-like and diverse sentence.
Throughout the experiments, we found out that our method not only works well
for English but can generalize on unseen languages as well. Data and code are
available at https://github.com/phkhanhtrinh23/LAMPAT.

## Introduction

Paraphrase generation involves the transformation of a given sentence or
phrase into its equivalent form while preserving its semantic content. By
leveraging the power of NLP techniques, paraphrase generation can enhance
several applications, such as machine translation (Freitag et al. 2020),
information retrieval (Lewis et al. 2020), question-answering systems (Gan and
Ng 2019), and text summarization (Cao et al. 2017). However, most existing
approaches in paraphrase generation focus on a single language such as
English, limiting their effectiveness in multilingual scenarios where accurate
and contextually appropriate paraphrases are crucial.

Figure 1: The training process of LAMPAT consists of multiple stages. Firstly,
we create a synthetic parallel corpus using unsupervised monolingual data.
Next, we utilize LoRA to effectively fine-tune our model. Finally, we obtain
the self-supervised model LAMPAT through the utilization of Virtual
Adversarial Training.

Most of the current approaches for multilingual paraphrasing are built around
the mechanism of machine translation. Thompson and Post (2020b) make use of
the multilingual neural machine translation (MNMT) model from Thompson and
Post (2020a) to translate between the same source and target languages with a
special decoding algorithm to reduce the lexical overlaps, creating a
paraphrase of the input text. Meanwhile, round-trip translation (Federmann,
Elachqar, and Quirk 2019), which creates a pivot language for source language
to translate forth and back to create a different wording output, is another
approach for paraphrase generation.

While the machine translation approach is widely recognized for its
effectiveness in multilingual paraphrasing, it encounters various limitations.
One primary hurdle is the necessity of obtaining reliable parallel corpora of
high quality for machine translation, which can be challenging to obtain in
real-world scenarios. Another drawback lies in the lack of diversity in the
generated output, often resulting in output sequences that closely resemble
the input and fail to preserve crucial information present in the original
input, as depicted in Table 1. Many MNMT models utilize heuristic blocking
techniques during inference to avoid generating output sequences that are
identical to the input. However, this approach limits the model’s ability to
alter word order or employ diverse syntactic structures, such as inversion or
active-to-passive transformations. Consequently, the generation lacks
diversity. Furthermore, these blocking algorithms heavily rely on the
distribution of the vocabulary. For instance, Thompson and Post (2020b) reduce
the probability of selecting subsequent subwords in $n$-grams to encourage the
model to choose different subwords. Nevertheless, these alternative subwords
may include antonyms or unrelated words, potentially leading to paraphrases
with semantic meanings deviating from the intended direction and producing
inappropriate results.

Type | Sentence
---|---
Input | I like to eat pasta.
Human reference |  My favourite food is pasta.
Lack of diversity | I like eating pasta.
Incorrect meaning | I like to eat paste.
Table 1: While human can change the structure of the sentence to create
paraphrase (in bold), paraphrasing model usually tends to replace words or
slightly modify the syntax (in italic). In the worst case, paraphrasing model
even changes the meaning of the sentence by using inappropriate word
replacement (in underline).

In this work, we propose LAMPAT (Low-rank Adaptation for Multilingual
Paraphrasing using Adversarial Training) as an approach to mitigate the strict
requirements of parallel corpora by learning in an unsupervised manner and
alleviate the problem of duplicate generation by using adversarial training
objectives. According to Figure 1, to eliminate the need for parallel corpora,
we use the monolingual dataset and apply a series of processes (i) identify
stop words (ii) remove stop words (iii) randomly shuffle the words to create
the corrupted version of the original input. The training of the model focuses
on the objective of reconstructing the original sentence from a corrupted
version, aiming to recreate the initial sentence. However, by this learning
objective, the model is drawn to generate the same sentence compared to the
original sentence, which does not create a syntactically diverse paraphrase.
To cope with this problem, we further propose using Virtual Adversarial
Training (VAT) (Zhang et al. 2019) and noise perturbation added directly to
the input embedding to steer the model towards a more diverse paraphrase
generation, as in Figure 2. In addition, Large Language Models (LLMs), are
known to experience the catastrophic forgetting (Kaushik et al. 2021) during
full fine-tuning; therefore, we adapt LoRA (Hu et al. 2021) as a parameter-
efficient fine-tuning method to partially update the model’s prior knowledge
and preserve all the linguistic knowledge on which the model has been pre-
trained. In general, LAMPAT can effectively generate human-like paraphrases
while preserving the original semantic meaning and employing different
syntactic structures to promote the diversity of the predictions.

In summary, the key contributions of this paper are as followed:

  * •

To resolve the requirement of parallel corpora for machine translation, we
propose the unsupervised learning method for multilingual paraphrasing.

  * •

To address the issue of predominantly generating identical or highly lexical-
similar outputs, we incorporate noise perturbation and a virtual labeling
strategy into the adversarial training process, aiming to alleviate this
limitation.

  * •

We expand the multilingual paraphrasing evaluation dataset to include more
languages and leverage future research in multilingual paraphrase generation.

## Methodology

Figure 2: LAMPAT is showcased using actual inputs. Initially, an input text
undergoes corruption by removing stopwords and shuffling. Then, a noise
perturbation, denoted as $\delta$, is introduced into the text embedding to
generate a paraphrase that exhibits lexical diversity. The transformer block
is replicated $N$ times, with the Multi-Head Attention component decomposed
into low-rank matrices for efficient fine-tuning. Lastly, LAMPAT is trained
using virtual adversarial training, incorporating a two-component loss
function: the reconstruction loss $\mathcal{L}_{rec}$ and the virtual
adversarial regularizer $\mathcal{L}_{vadv}$.

The training process of our paraphrasing model is illustrated in Figure 1,
comprising three essential elements: Synthetic Parallel Corpora, Parameter-
Efficient Fine-Tuning (PEFT), and Virtual Adversarial Training (VAT). We
employ the Self-supervised model to generate paraphrases, which not only
addresses the data shortage issue using unsupervised learning but also
maintains semantic similarity and enhances lexical diversity as well.

### Problem Definition

Given 2 sentences $x$ and $y$, where $x$ is the original sentence, and $y$ is
the paraphrase reference of $x$. Let $\mathcal{M}(x)$ be the meaning of $x$
and $\mathcal{S}(x,y)$ be the lexical or syntactic similarity between $x$ and
$y$. $\hat{y}=\underset{y}{argmax}[p(y|M(x))-S(x,y)]$. The term $p(y|M(x))$ is
the probability of generating $y$ that conveys the same meaning as $x$ (i.e.
$M(x)$). $S(x,y)\in[0,1]$ measures the lexical similarity of $x$ and $y$,
where $S(z,z)=1$ for every $z$. Based on this formulation, paraphrasing should
be a method which not only allows us to convey the same meaning but also
enhances lexical diversity in the generated text.

### Synthetic Parallel Corpora

Synthetic Parallel Corpora is a significant component of LAMPAT. First, we
corrupt the input by removing the stop words in the sentence, and then
randomly shuffling the words in the remaining text. The corrupted sentence is
referred to as the source sequence $S$, while the original uncorrupted
sentence is the target sequence $T$. We have a set of stop words $A$, which
are removed from the sentences. Our goal is to generate the paraphrase by
reconstructing $T$ from the keywords or the corrupted sentence $S$, where
$S=\text{Shuffle}(T-A)$. When we fine-tune the model, we create the input
sequence $X$ by combining the source and target sequences with a special
symbol in between. The input sequence $X$ is represented as
$(x_{1},x_{2},x_{3},...,x_{k}$, \n, $x_{k+1},x_{k+2},...x_{h})$, where the
source sequence is denoted as $S=(x_{0},x_{1},...,x_{k})$, and the target
sequence is denoted as $T=(x_{k+1},x_{k+2},...,x_{h})$. The special character
\n is included to differentiate between the source and target tokens, and it
also serves as a prompt during the inference process.

### Parameter-Efficient Fine-Tuning

The method chosen for parameter-efficient fine-tuning is Low-rank Adaptation
(LoRA) (Hu et al. 2021). A major drawback of fine-tuning is that the resulting
model contains the same number of parameters as the original model. LoRA
overcomes this by indirectly training some dense layers in a neural network by
optimizing rank decomposition matrices of the dense layers’ changes during
adaptation while keeping the pre-trained weights frozen. This also allows for
efficient task-switching by simply replacing the matrices, resulting in
reduced storage requirements and task-switching overhead.

LoRA (Hu et al. 2021) introduces a constraint on the update of a pre-trained
matrix $W_{x}\in\mathbb{R}^{d\times k}$ using a low-rank matrix
$W_{\beta}\in\mathbb{R}^{d\times k}$. Instead of directly updating $W_{x}$,
LoRA modifies it as $W^{\prime}_{x}=W_{x}+W_{\beta}$ and focuses on updating
the parameters involved in the construction of $W_{\beta}$. The construction
of $W_{\beta}$ involves the multiplication of two matrices
$B\in\mathbb{R}^{d\times r}$ and $A\in\mathbb{R}^{r\times k}$, where $r\ll
min(d,k)$, resulting in a low-rank matrix $W^{\prime}_{x}=W_{x}+BA$. By using
LoRA, we reduce the number of parameters to tune from $d\times k$ to
$r\times(d+k)$. Specifically, we follow the application of LoRA to the query
and value transformation matrices in the multi-head attention sublayers, as in
Hu et al. (2021).

### Virtual Adversarial Training

Consider a standard classification task with an underlying data distribution
$D$ over examples $x\in\mathbb{R}^{d}$ and corresponding labels $y$. Assume
that we are given a suitable loss function $\mathcal{L}_{rec}$, for example,
the cross-entropy loss for a neural network. As usual,
$\theta\in\mathbb{R}^{p}$ is the set of model parameters. Our goal is then to
find model parameters $\theta$ that satisfy:

$\min_{\theta}E_{(x,y)\sim\mathcal{D}}[\mathcal{L}_{rec}(f(x,\theta),y)]$ (1)

Once we have created the synthetic parallel data, we define a set of
permissible perturbations $\mathcal{C}\subseteq\mathbb{R}^{d}$ for each data
point $x$, which represents the manipulative ability of the adversary. Instead
of directly using samples from the distribution $\mathcal{D}$ in the loss
$\mathcal{L}$, we allow the adversary $\delta\in\mathcal{C}$ to first perturb
the input embedding. Specifically, $x$ denotes the sub-word embedding in
$f(x,\theta)$. We notice that by perturbing the embedding space $x+\delta$,
rather than the input space, adversarial training may unintentionally favour
on-manifold perturbations over regular perturbations, leading to improved
generalization. Hence, we apply perturbations to the embedding space. On the
other hand, complete label information may not always be available, and
especially in this unsupervised manner, we aim to output virtual labels other
than the label $y$. Consequently, we adopt a strategy to replace the label $y$
with its current approximation, $f(x,\theta)$. This approximation is not
necessarily naive, as $f(x,\theta)$ tends to be close to $y$ when the number
of labelled training samples is large. This rationale also explains the use of
the term “virtual” in Miyato et al. (2018). Essentially, we employ virtual
labels generated from $f(x,\theta)$ in place of paraphrasing labels, and
compute the adversarial direction based on these virtual labels. As a result,
we replace the Equation 1 by:

$\begin{split}\min_{\theta}E_{(x,y)\sim\mathcal{D}}\left[\max_{\delta}\left[\mathcal{L}_{rec}(f(x+\delta,\theta),y)\right.\right.\\\
\left.\left.+\alpha\mathcal{L}_{vadv}(f(x+\delta,\theta),f(x,\theta))\right]\vphantom{\max_{\delta}}\right]\end{split}$
(2)

Algorithm 1 Low-rank Adaptation Multilingual Paraphrasing using Adversarial
Training.

Input: $X$: Training samples, $f(x;\theta)$: the machine learning model
parameterized by $\theta$, $\epsilon$: the perturbation bound, $\tau$: the
global learning rate, $\alpha$: the smoothing proportion of the adversarial
training, $\eta$: the ascent step size, $\mathcal{H}$: the Hessian gradient
matrix, $N$: the number of epochs, $K$: the number of ascent steps, $e^{*}$:
the number of epochs trained with PGD algorithm, $\gamma$: the scaling factor.
Output: $\theta$.

1: for epoch = 1…$N$ do

2: for minibatch $B$ $\in$ X do

3: $\delta\sim\gamma\cdot\mathcal{N}(0,\sigma^{2}I)$

4: for m = 1…K do

5: Accumulate gradient of parameter $\theta$:

6: $g_{m}\leftarrow g_{m-1}+\frac{1}{K}E_{(x,y)\in
B}[\nabla_{\theta}\mathcal{L}_{rec}(f(x+\delta,\theta),y)+\alpha\nabla_{\theta}\mathcal{L}_{vadv}(f(x+\delta,\theta),f(x,\theta))]$

7: Calculate the gradient of the perturbation $\delta$:

8:
$g_{adv}\leftarrow\nabla_{\delta}\mathcal{L}_{vadv}(f(x+\delta,\theta),f(x,\theta))$

9: Update the perturbation $\delta$ through gradient ascent:

10: if epoch $\leq$ $e^{*}$ then

11:
$\delta\leftarrow\prod_{||\delta||\leq\epsilon}(\delta+\eta\frac{g_{adv}}{||g_{adv}||_{F}})$

12: else

13:
$\delta\leftarrow\prod_{\mathcal{H}}(\delta+\eta[\mathcal{H}]^{-1}\frac{g_{adv}}{||g_{adv}||_{F}})$

14: end if

15: end for

16: Update the parameter $\theta$ through gradient descent:

17: $\theta\leftarrow\theta-\tau g_{K}$

18: end for

19: end for

In reference to Algorithm 1, our method can reduce the computational cost of
adversarial training using projection over constraints algorithms. It achieves
similar levels while conducting fewer sweeps of forward and backward
propagation, making it faster and less computationally expensive. Our method
takes advantage of every propagation to update weights and allows multiple
updates per iteration, potentially leading to faster convergence. According to
Zhang et al. (2019), by combining these factors, the $K$ ascent steps
significantly accelerate standard adversarial training. After that, the
model’s parameter $\theta$ is updated all at once with the accumulated
gradients. By taking a descent step along the $K$ gradients, we can
approximately calculate the following objective function:

$\begin{split}\min_{\theta}E_{x,y\sim\mathcal{D}}\left[\frac{1}{K}\sum_{t=0}^{K-1}\max_{\delta}[\mathcal{L}_{rec}(f(x+\delta,\theta),y)\right.\\\
\left.+\alpha\mathcal{L}_{vadv}(f(x+\delta,\theta),f(x,\theta))]\vphantom{\frac{1}{K}}\right]\end{split}$
(3)

The rationale behind computing $g(\delta)$ with respect to the virtual
adversarial regularizer $\nabla_{\delta}\mathcal{L}_{vadv}$ instead of
$\nabla_{\delta}\mathcal{L}_{rec}$ in Algorithm 1 is due to the unsupervised
nature of the model training and the objective to guide the model towards the
virtual labels rather than reconstructing the original sentence, which could
result in duplication. Equation 3 is essentially replacing the original batch
$x$ with a virtual batch that is $K$ times larger, comprising samples with
embeddings of $X+\delta_{0},...,X+\delta_{K-1}$. While the original
adversarial training Equation 2 minimizes the maximum risk at a single
estimated point in the vicinity of each training sample, Equation 3 minimizes
the maximum risk at each ascent step and guides the model towards the virtual
labels with minimal additional overhead. Moreover, we use both Projected
Gradient Descent (PGD) and Projected-Newton Method (PNM) in Algorithm 1. The
Gradient Descent step involves descending along the linear estimate of the
function, while Newton’s step involves moving the point towards the minimum of
the parabola that approximates the function, which can lead to faster
convergence.

## Experimental Setup

### Primary Model

Pre-trained model: We utilize the mGPT model, a pre-trained multilingual GPT-
like model, with a 1.3B checkpoint consisting of 100K tokens and supporting 61
languages. The LoRA implementation from PEFT
111https://github.com/huggingface/peft.

Dataset: To assess the fine-tuning, we choose to use the latest version WMT19
(Foundation 2019) to train the model. This dataset covers a wide range of 15
languages including Arabic, Czech, German, English, Spanish, French, Hindi,
Indonesian, Italian, Japanese, Kazakh, Dutch, Portuguese, Russian, and
Chinese. The WMT19 dataset we use is in its latest version, which is just
released in 2023. To balance language resources, we employ a uniform
distribution to sample sentences, creating a training set of nearly 600k
sentences and a validation set of around 100k sentences. For training, we use
the monolingual version of WMT19 and corrupt the input sentence by removing
all of the stop words 222https://github.com/stopwords-iso/stopwords-iso,
further we randomly shuffle the words 33% of the time. The goal is to
reconstruct a sentence from its keywords, or its corrupted sentence. For the
machine translation approach, we use the available bilingual version available
in WMT19, we sample with the same strategy as LAMPAT’s training with 600k
samples for training and 100k for validation.

### Baseline Model

We compare our method with other approaches such as multilingual machine
translation (MMT) proposed in Thompson and Post (2020b) and denoising auto-
encoder (DAE) in Guo et al. (2019). Initially, these methods are trained on
different datasets, thus, we re-train them, following the procedures proposed
in each paper, on the WMT19 dataset.

### Evaluation Dataset

We follow Guo et al. (2019), randomly select 10k sentences respectively from
each language of English, Spanish, Russian, and Chinese to construct the test
set. However, the number of languages covered in Guo et al. (2019) is
relatively small compared to the number of languages around the world, over
7000 languages 333https://www.ethnologue.com/statistics. Therefore, we expand
the languages covered in the task of multilingual paraphrasing to 13
languages, including some figurative languages such as Japanese, and Chinese
or accented languages such as Vietnamese. The proposed evaluation dataset has
two types to assess different aspects of the model: Input-only and Input-
reference.

Input-only: We follow Thompson and Post (2020b) and use the validation set
from the WMT19 dataset, which is realeased in 2019 and available on
HuggingFace 444https://huggingface.co/datasets/wmt19. Since the WMT19 dataset
is used for the task Machine Translation, we thus, extract one side of the
dataset in order to build the Input-only evaluation dataset. The languages we
extracted from the WMT19 are Czech (cs), German (de), English (en), Finish
(fi), French (fr) and Chinese (zh). Most of the samples in this evaluation
dataset are of the news domains, which is used to test the model’s ability on
producing a paraphrase that conveys the same meaning.

Input-reference: We use in total 3 datasets to construct this evaluation set:

  * •

PAWS-X (Yang et al. 2019) is the cross-lingual paraphrase identification
dataset, thus, we extract only the sentence pairs with label 1 (indicating
paraphrase) for 6 languages: Japanese (ja), Chinese (zh), German (de), French
(fr), Spanish (es) and English (en).

  * •

Opusparcus (Creutz 2019): is a paraphrase corpus for six European languages:
German (de), English (en), Finnish (fi), French (fr), Russian (ru), and
Swedish (sv). We extract the test set of Opusparcus with a score of 4 to
ensure high-quality sentence pairs.

  * •

STAPLE (Duolingo 2020): is a multi-reference machine translation dataset in
which each reference could be viewed as the paraphrase. Since STAPLE does not
have the validation or test set, we randomly extract 1000 samples, each with 5
reference texts, covering 3 languages: Vietnamese (vi), Portuguese (pt) and
Hungarian (hu), to construct the first multi-reference paraphrase generation
corpus with 1000 samples and 4 references each.

### Automatic Evaluation

We follow Chowdhury, Zhuang, and Wang (2022) to report the result on BLEU
(Papineni et al. 2002), Self-BLEU, Self-TER, which adapted TER (Snover et al.
2006) to the input instead of the reference, BERTScore (Zhang et al. 2020)
with bert-base-multilingual-cased checkpoint, iBLEU with $\alpha=0.7$
following Hosking and Lapata (2021). In addition, we further use two latest
paraphrase metrics: ParaScore (Shen et al. 2022) and BERT-iBLEU (Niu et al.
2021).

### Human Evaluation

In addition to machine evaluation, we also conduct the human evaluation of the
paraphrase generated by our model. For each of the following languages:
English, Vietnamese, German, French and Japanese, we randomly extract 200
sentence triples of the input sentence, our model prediction and the output
from the model of Thompson and Post (2020b). For each mentioned language, we
ask 5 annotators to score 200 sentence pairs independently. Each annotator is
instructed to rate on the 1-5 scale (with 5 being the highest) based on 3
criteria: (i) Semantic preservation, evaluating how much information is
preserved in the output; (ii) Lexical similarity, evaluating how much similar
in term of syntax or word choices of the output compared to input; and (iii)
Fluency, assessing the fluency and coherence of the generated output. The
annotators’ agreement is measured using Krippendorff’s alpha (Krippendorff
1970), which provides a measure of inter-annotator reliability.

## Results

### Main Results

Method | en | es | zh | ru
---|---|---|---|---
BERTScore $\uparrow$
DAE | 79.25 | 80.56 | 78.91 | 77.83
MMT | 84.93 | 82.68 | 84.79 | 81.32
LAMPAT | 86.86 | 84.35 | 83.26 | 84.01
Self-BLEU $\downarrow$
DAE | 20.35 | 30.49 | 20.38 | 20.61
MMT | 19.89 | 28.57 | 10.32 | 15.19
LAMPAT | 19.46 | 20.16 | 14.95 | 12.46
Self-TER $\uparrow$
DAE | 50.48 | 45.19 | 45.92 | 48.31
MMT | 52.45 | 43.16 | 41.25 | 45.68
LAMPAT | 61.32 | 55.43 | 55.28 | 50.67
BERT-iBLEU $\uparrow$
DAE | 79.33 | 78.08 | 79.05 | 78.14
MMT | 83.92 | 80.16 | 85.72 | 81.99
LAMPAT | 85.52 | 83.41 | 83.61 | 84.69
ParaScore $\uparrow$
DAE | 88.75 | 90.46 | 89.37 | 88.50
MMT | 89.45 | 91.56 | 90.05 | 88.47
LAMPAT | 92.95 | 92.96 | 90.64 | 91.92
Table 2: Multilingual paraphrase generation test results over 4 languages English, Spanish, Chinese and Russian from the work of DAE (Guo et al. 2019). Method | fr | cs | fi | de | en | zh
---|---|---|---|---|---|---
BERTScore $\uparrow$
DAE | 70.39 | 78.45 | 80.96 | 69.32 | 65.68 | 81.35
MMT | 74.40 | 80.20 | 81.20 | 71.20 | 65.48 | 80.29
LAMPAT | 85.42 | 83.92 | 84.91 | 86.16 | 82.59 | 79.18
Self-TER $\uparrow$
DAE | 45.53 | 39.97 | 54.50 | 46.18 | 30.42 | 28.49
MMT | 48.40 | 37.12 | 56.74 | 48.42 | 33.76 | 31.56
LAMPAT | 49.53 | 40.49 | 63.92 | 53.76 | 39.31 | 40.80
BERT-iBLEU $\uparrow$
DAE | 71.10 | 79.56 | 81.38 | 71.24 | 67.36 | 81.05
MMT | 75.48 | 78.81 | 82.15 | 74.16 | 66.08 | 69.70
LAMPAT | 84.63 | 85.22 | 86.03 | 85.38 | 84.08 | 79.90
ParaScore $\uparrow$
DAE | 89.91 | 88.46 | 75.93 | 87.43 | 88.56 | 88.12
MMT | 89.95 | 90.10 | 81.15 | 89.47 | 91.49 | 90.35
LAMPAT | 93.24 | 92.45 | 86.82 | 93.21 | 94.84 | 92.24
Table 3: Multilingual paraphrase generation test results on our input-only evaluation dataset. Method | sv | fi | en | fr | de | ru | ja | zh | es | hu | pt | vi
---|---|---|---|---|---|---|---|---|---|---|---|---
BERTScore $\uparrow$
DAE | 80.20 | 79.50 | 84.00 | 85.20 | 89.44 | 89.90 | 74.50 | 83.00 | 88.80 | 75.30 | 78.40 | 80.50
MMT | 83.48 | 84.92 | 75.86 | 77.01 | 82.53 | 81.49 | 72.61 | 71.93 | 82.45 | 73.11 | 83.44 | 75.76
LAMPAT | 85.47 | 85.47 | 94.87 | 89.94 | 90.10 | 92.16 | 90.77 | 81.00 | 92.87 | 86.16 | 91.99 | 87.70
BLEU $\uparrow$
DAE | 5.22 | 4.20 | 10.50 | 15.51 | 14.44 | 5.46 | 20.62 | 47.80 | 18.58 | 8.10 | 16.01 | 15.05
MMT | 1.22 | 0.39 | 11.79 | 6.10 | 11.98 | 0.77 | 9.55 | 16.69 | 20.49 | 0.84 | 11.09 | 3.58
LAMPAT | 6.04 | 4.90 | 23.07 | 20.55 | 21.65 | 6.48 | 29.13 | 30.52 | 26.95 | 8.16 | 16.65 | 20.24
Self-BLEU $\downarrow$
DAE | 24.58 | 17.56 | 45.50 | 30.20 | 25.75 | 36.82 | 50.78 | 50.22 | 40.69 | 25.82 | 30.44 | 27.26
MMT | 25.72 | 16.76 | 50.33 | 23.47 | 35.65 | 39.67 | 45.62 | 63.32 | 36.65 | 27.82 | 34.70 | 27.24
LAMPAT | 23.68 | 14.53 | 43.47 | 19.55 | 24.54 | 30.25 | 40.56 | 47.58 | 29.00 | 22.94 | 31.43 | 25.51
iBLEU $\uparrow$
DAE | 0.05 | 0.01 | -0.04 | -0.05 | -0.01 | 0.01 | 0.08 | 0.10 | 0.04 | 0.05 | 0.01 | -0.01
MMT | 0.01 | 0.02 | 0.02 | 0.04 | 0.04 | -0.03 | 0.06 | 0.12 | 0.09 | 0.06 | 0.11 | 0.14
LAMPAT | 0.08 | 0.04 | 0.04 | 0.02 | 0.05 | 0.04 | 0.09 | 0.15 | 0.13 | 0.15 | 0.15 | 0.05
Self-TER $\uparrow$
DAE | 50.18 | 60.57 | 24.55 | 55.34 | 57.44 | 33.14 | 40.27 | 30.75 | 31.55 | 40.57 | 47.81 | 50.55
MMT | 47.16 | 53.59 | 23.90 | 47.07 | 43.48 | 29.27 | 40.88 | 26.85 | 30.97 | 36.88 | 27.10 | 39.95
LAMPAT | 57.20 | 59.45 | 56.56 | 56.49 | 61.33 | 47.80 | 42.62 | 35.34 | 38.55 | 41.51 | 55.11 | 57.91
BERT-iBLEU $\uparrow$
DAE | 70.20 | 78.20 | 66.60 | 56.20 | 62.70 | 68.20 | 52.00 | 67.20 | 57.20 | 67.50 | 68.80 | 77.20
MMT | 61.43 | 68.44 | 65.25 | 66.33 | 63.07 | 69.06 | 65.33 | 62.67 | 63.95 | 72.39 | 72.34 | 69.91
LAMPAT | 73.30 | 78.60 | 67.60 | 82.33 | 80.25 | 69.98 | 78.04 | 67.72 | 86.46 | 76.32 | 78.92 | 81.39
ParaScore $\uparrow$
DAE | 82.00 | 72.80 | 85.00 | 85.30 | 88.90 | 90.20 | 84.50 | 87.60 | 90.10 | 80.50 | 80.52 | 80.02
MMT | 83.17 | 81.72 | 76.83 | 77.62 | 82.59 | 81.69 | 73.92 | 73.07 | 82.70 | 74.54 | 83.84 | 77.01
LAMPAT | 86.42 | 86.00 | 94.51 | 90.73 | 90.29 | 91.47 | 91.33 | 93.01 | 92.96 | 86.89 | 91.74 | 89.78
Table 4: Multilingual paraphrase generation test results on our input-
reference evaluation dataset.

For the test dataset from Guo et al. (2019) and our input-only evaluation, we
evaluate using BERTScore, Self-BLEU, Self-TER, BERT-iBLEU and ParaScore, as
these evaluation metrics do not require the reference text. The test results
are depicted by Table 2 and 3. For input-reference, we report BLEU and iBLEU,
in addition, as depicted by Table 4. LAMPAT can generate diverse output
compared to the input, which is indicated by the low score in Self-BLEU and
high score in Self-TER. LAMPAT also preserves better information from the
input as demonstrated by the BERTScore results. We have manually examined the
text generated by all four methods by randomly selecting 100 samples per
language. Although these sentences convey the same meaning, the word choices
and syntax structures are largely different from both input and reference,
leading to the low score of iBLEU. Overall, BERT-iBLEU and ParaScore metrics,
which are the metrics that grade both lexical and semantic aspects, show that
all three methods could generate comprehensive sentences. However, our method,
LAMPAT, still achieves the highest score over 13 languages we have tested.
Even though our model learned in an unsupervised manner, it can still
outperform the supervised counterpart. Overall, LAMPAT has demonstrated that
integrating adversarial training into unsupervised learning could improve the
performance of multilingual paraphrase generation.

### Human Evaluation Results

The results of our human evaluation can be found in Table 5, where we present
the assessments provided by human evaluators. The evaluations conducted by
human experts provide valuable insights into the performance of our model.

Method | SP $\uparrow$ | LS $\downarrow$ | F $\uparrow$
---|---|---|---
MMT | 3.2 | 4.8 | 4.8
LAMPAT | 4.2 | 3.2 | 4.8
Human generated paraphrase | 4.4 | 2.4 | 4.8
Krippendorff’s alpha | 0.68 | 0.7 | 0.78
Table 5: Human evaluation results. SP: Semantic Preservation; LS: Lexical
Similarity; F: Fluency. All the scores reported are the average value of 5
chosen languages.

## Ablation Study

### Parameter-Efficient Fine-Tuning

In order to examine which ingredients help improve the performance of LAMPAT,
we experiment with Adapter (Houlsby et al. 2019), Prefix Tuning (Li and Liang
2021), Prompt Tuning (Lester, Al-Rfou, and Constant 2021), LoRA (Hu et al.
2021), P-Tuning (Liu et al. 2022) and Full fine-tuning to find out which PEFT
methods result in a more stable and better result. For each of the methods, we
train on the same number of epochs with 3 random seeds and report the mean and
standard deviation of ParaScore of all languages.

Method | ParaScore
---|---
Prefix Tuning | 77.45$\pm$8.5
Prompt Tuning | 82.25$\pm$7.9
P-Tuning | 82.94$\pm$5.5
Adapter | 89.56$\pm$3.2
LoRA | 90.67$\pm$2.5
Full | 89.46$\pm$1.5
Table 6: ParaScore of different fine-tuning methods on 13 languages. The mean
and standard deviation are the weighted mean and standard deviation in all 13
languages.

According to Table 6, LoRA experiences to be the most stable method and
achieves the highest scores, especially for generation tasks.

### Adversarial Optimization

Figure 3: The average ParaScore of each technique over 13 languages.

Since LoRA is a stable method for partially fine-tuning LLMs, we adapt LoRA as
the PEFT method for fine-tuning mGPT. In order to study the effect of
adversarial training in unsupervised learning, we employ two settings which
are Projected Gradient Descent (PGD) and Projected Newton Method (PNM)
(Bertsekas 1982) together with LoRA. According to Figure 3, we hypothesized
that the main driving force that makes LoRA + PGD + PNM better compared to PGD
is because it uses the Taylor expansion, which has a better approximation of
the objective function. In general, when we are near the local optima (for
example, at the last 2 epochs), we can take a few more Newton’s steps to reach
the optimum point instead of taking many small Gradient Descent steps.

## Related Works

### Multilingual Paraphrasing

Numerous techniques for paraphrasing in multiple languages employ Machine
Translation-based models (MNMT). To illustrate, Thompson and Post (2020b)
applied a pretrained MNMT model introduced by Thompson and Post (2020a), along
with a customized decoding algorithm aimed at reducing repetition of words and
encouraging diverse vocabulary usage. Another instance is the work by Guo et
al. (2019), which leveraged a language model pretraining task adapted from
Conneau and Lample (2019). During inference, the same language code is
provided to the model, and the sequence is generated sequentially in an
autoregressive manner. Despite the fact that the translation-based approach
produces high-quality and fluent paraphrases, it faces certain inherent
challenges. Firstly, there is a potential for bias stemming from the dominance
of certain languages used for training the MNMT model. Secondly, there may be
instances of incorrect synthetic paraphrases due to the inherent ambiguity in
the pivot language, as pointed out by Thompson and Post (2020b). These issues
need careful consideration in the development of multilingual paraphrasing
methods.

### Adversarial Training

Adversarial training is a powerful technique employed in the development of
resilient neural networks. While the computer vision community, as highlighted
by Goodfellow, Shlens, and Szegedy (2015), has generally accepted that
adversarial training can have a detrimental impact on model generalization,
the scenario appears to be quite different for language models, as evidenced
by studies such as Pereira et al. (2020) and Dong et al. (2021a).
Incorporating adversarial training into large language models (LLMs), as
explored by Miyato et al. (2018) and further elaborated upon by Dong et al.
(2021b), has been found to yield improvements in both model generalization and
robustness. An innovative training algorithm, denoted as YOPO (You Only
Propagate Once), was proposed by Zhang et al. (2019). YOPO takes advantage of
the “free” training strategies advocated by Shafahi et al. (2019) to diversify
the training data by incorporating various adversarial samples while imposing
different norm constraints. Also, Miyato et al. (2018) proposes a new training
method that regularizes the training objective by using virtual labels in
adversarial training. These approaches collectively showcase the effectiveness
of adversarial training in enhancing both the robustness and generalization
capabilities of neural networks.

## Conclusion and Future Work

In this research, we introduce an efficient method for generating paraphrases
in multiple languages using Low-Rank Adaptation combined with virtual labeling
during adversarial training. Importantly, our approach delivers satisfactory
results without relying on supervised learning. Additionally, we contribute to
the field by creating a novel multilingual multi-domain evaluation dataset.
While LAMPAT has demonstrated proficiency in generating human-like paraphrases
across various languages, it still requires improvements in handling idiomatic
expressions. Furthermore, our evaluation dataset covers only 13 languages,
leaving out many, especially those with limited resources. This highlights the
ongoing need for research to enhance and expand the capabilities of
multilingual paraphrase generation models.

## Acknowledgement

This research/project is supported by the National Research Foundation,
Singapore under its AI Singapore Programme, AISG Award No: AISG2-TC-2022-005.

## References

  * Bertsekas (1982) Bertsekas, D. P. 1982.  Projected Newton Methods for Optimization Problems with Simple Constraints.  _SIAM Journal on Control and Optimization_ , 20(2): 221–246.
  * Cao et al. (2017) Cao, Z.; Luo, C.; Li, W.; and Li, S. 2017.  Joint Copying and Restricted Generation for Paraphrase.  In _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence_ , AAAI’17, 3152–3158. AAAI Press.
  * Chowdhury, Zhuang, and Wang (2022) Chowdhury, J. R.; Zhuang, Y.; and Wang, S. 2022.  Novelty Controlled Paraphrase Generation with Retrieval Augmented Conditional Prompt Tuning.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(10): 10535–10544.
  * Conneau and Lample (2019) Conneau, A.; and Lample, G. 2019.  Cross-lingual Language Model Pretraining.  In Wallach, H.; Larochelle, H.; Beygelzimer, A.; d'Alché-Buc, F.; Fox, E.; and Garnett, R., eds., _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc.
  * Creutz (2019) Creutz, M. 2019.  Open subtitles paraphrase corpus for six languages.  _LREC 2018 - 11th International Conference on Language Resources and Evaluation_ , (2005): 1364–1369.
  * Dong et al. (2021a) Dong, X.; Luu, A. T.; Ji, R.; and Liu, H. 2021a.  Towards robustness against natural language word substitutions.  _arXiv preprint arXiv:2107.13541_.
  * Dong et al. (2021b) Dong, X.; Luu, A. T.; Lin, M.; Yan, S.; and Zhang, H. 2021b.  How should pre-trained language models be fine-tuned towards adversarial robustness?  _Advances in Neural Information Processing Systems_ , 34: 4356–4369.
  * Duolingo (2020) Duolingo. 2020.  Data for the 2020 Duolingo Shared Task on Simultaneous Translation And Paraphrase for Language Education (STAPLE).
  * Federmann, Elachqar, and Quirk (2019) Federmann, C.; Elachqar, O.; and Quirk, C. 2019.  Multilingual Whispers: Generating Paraphrases with Translation.  In _Proceedings of the 5th Workshop on Noisy User-generated Text (W-NUT 2019)_ , 17–26. Hong Kong, China: Association for Computational Linguistics.
  * Foundation (2019) Foundation, W. 2019.  ACL 2019 Fourth Conference on Machine Translation (WMT19), Shared Task: Machine Translation of News.
  * Freitag et al. (2020) Freitag, M.; Foster, G.; Grangier, D.; and Cherry, C. 2020.  Human-Paraphrased References Improve Neural Machine Translation.  In _Proceedings of the Fifth Conference on Machine Translation_ , 1183–1192. Online: Association for Computational Linguistics.
  * Gan and Ng (2019) Gan, W. C.; and Ng, H. T. 2019.  Improving the Robustness of Question Answering Systems to Question Paraphrasing.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , 6065–6075. Florence, Italy: Association for Computational Linguistics.
  * Goodfellow, Shlens, and Szegedy (2015) Goodfellow, I. J.; Shlens, J.; and Szegedy, C. 2015.  Explaining and Harnessing Adversarial Examples.  arXiv:1412.6572.
  * Guo et al. (2019) Guo, Y.; Liao, Y.; Jiang, X.; Zhang, Q.; Zhang, Y.; and Liu, Q. 2019.  Zero-Shot Paraphrase Generation with Multilingual Language Models.  arXiv:1911.03597.
  * Hosking and Lapata (2021) Hosking, T.; and Lapata, M. 2021.  Factorising Meaning and Form for Intent-Preserving Paraphrasing.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , 1405–1418. Online: Association for Computational Linguistics.
  * Houlsby et al. (2019) Houlsby, N.; Giurgiu, A.; Jastrzebski, S.; Morrone, B.; De Laroussilhe, Q.; Gesmundo, A.; Attariyan, M.; and Gelly, S. 2019.  Parameter-Efficient Transfer Learning for NLP.  In Chaudhuri, K.; and Salakhutdinov, R., eds., _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , 2790–2799. PMLR.
  * Hu et al. (2021) Hu, E. J.; Shen, Y.; Wallis, P.; Allen-Zhu, Z.; Li, Y.; Wang, S.; Wang, L.; and Chen, W. 2021.  LoRA: Low-Rank Adaptation of Large Language Models.  arXiv:2106.09685.
  * Kaushik et al. (2021) Kaushik, P.; Gain, A.; Kortylewski, A.; and Yuille, A. 2021.  Understanding Catastrophic Forgetting and Remembering in Continual Learning with Optimal Relevance Mapping.  arXiv:2102.11343.
  * Krippendorff (1970) Krippendorff, K. 1970.  Estimating the Reliability, Systematic Error and Random Error of Interval Data.  _Educational and Psychological Measurement_ , 30: 61 – 70.
  * Lester, Al-Rfou, and Constant (2021) Lester, B.; Al-Rfou, R.; and Constant, N. 2021.  The Power of Scale for Parameter-Efficient Prompt Tuning.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , 3045–3059. Online and Punta Cana, Dominican Republic: Association for Computational Linguistics.
  * Lewis et al. (2020) Lewis, M.; Ghazvininejad, M.; Ghosh, G.; Aghajanyan, A.; Wang, S.; and Zettlemoyer, L. 2020.  Pre-training via Paraphrasing.  In Larochelle, H.; Ranzato, M.; Hadsell, R.; Balcan, M.; and Lin, H., eds., _Advances in Neural Information Processing Systems_ , volume 33, 18470–18481. Curran Associates, Inc.
  * Li and Liang (2021) Li, X. L.; and Liang, P. 2021.  Prefix-Tuning: Optimizing Continuous Prompts for Generation.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , 4582–4597. Online: Association for Computational Linguistics.
  * Liu et al. (2022) Liu, X.; Ji, K.; Fu, Y.; Tam, W.; Du, Z.; Yang, Z.; and Tang, J. 2022.  P-Tuning: Prompt Tuning Can Be Comparable to Fine-tuning Across Scales and Tasks.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , 61–68. Dublin, Ireland: Association for Computational Linguistics.
  * Miyato et al. (2018) Miyato, T.; ichi Maeda, S.; Koyama, M.; and Ishii, S. 2018.  Virtual Adversarial Training: A Regularization Method for Supervised and Semi-Supervised Learning.  arXiv:1704.03976.
  * Niu et al. (2021) Niu, T.; Yavuz, S.; Zhou, Y.; Keskar, N. S.; Wang, H.; and Xiong, C. 2021.  Unsupervised Paraphrasing with Pretrained Language Models.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , 5136–5150. Online and Punta Cana, Dominican Republic: Association for Computational Linguistics.
  * Papineni et al. (2002) Papineni, K.; Roukos, S.; Ward, T.; and Zhu, W.-J. 2002.  Bleu: a Method for Automatic Evaluation of Machine Translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , 311–318. Philadelphia, Pennsylvania, USA: Association for Computational Linguistics.
  * Pereira et al. (2020) Pereira, L.; Liu, X.; Cheng, F.; Asahara, M.; and Kobayashi, I. 2020.  Adversarial Training for Commonsense Inference.  In _Proceedings of the 5th Workshop on Representation Learning for NLP_ , 55–60. Online: Association for Computational Linguistics.
  * Shafahi et al. (2019) Shafahi, A.; Najibi, M.; Ghiasi, M. A.; Xu, Z.; Dickerson, J.; Studer, C.; Davis, L. S.; Taylor, G.; and Goldstein, T. 2019.  Adversarial training for free!  In Wallach, H.; Larochelle, H.; Beygelzimer, A.; d'Alché-Buc, F.; Fox, E.; and Garnett, R., eds., _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc.
  * Shen et al. (2022) Shen, L.; Liu, L.; Jiang, H.; and Shi, S. 2022.  On the Evaluation Metrics for Paraphrase Generation.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , 3178–3190. Abu Dhabi, United Arab Emirates: Association for Computational Linguistics.
  * Snover et al. (2006) Snover, M.; Dorr, B.; Schwartz, R.; Micciulla, L.; and Makhoul, J. 2006.  A Study of Translation Edit Rate with Targeted Human Annotation.  In _Proceedings of the 7th Conference of the Association for Machine Translation in the Americas: Technical Papers_ , 223–231. Cambridge, Massachusetts, USA: Association for Machine Translation in the Americas.
  * Thompson and Post (2020a) Thompson, B.; and Post, M. 2020a.  Automatic Machine Translation Evaluation in Many Languages via Zero-Shot Paraphrasing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , 90–121. Online: Association for Computational Linguistics.
  * Thompson and Post (2020b) Thompson, B.; and Post, M. 2020b.  Paraphrase Generation as Zero-Shot Multilingual Translation: Disentangling Semantic Similarity from Lexical and Syntactic Diversity.  In _Proceedings of the Fifth Conference on Machine Translation_ , 561–570. Online: Association for Computational Linguistics.
  * Yang et al. (2019) Yang, Y.; Zhang, Y.; Tar, C.; and Baldridge, J. 2019.  PAWS-X: A Cross-lingual Adversarial Dataset for Paraphrase Identification.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 3687–3692. Hong Kong, China: Association for Computational Linguistics.
  * Zhang et al. (2019) Zhang, D.; Zhang, T.; Lu, Y.; Zhu, Z.; and Dong, B. 2019.  You Only Propagate Once: Accelerating Adversarial Training via Maximal Principle.  In Wallach, H.; Larochelle, H.; Beygelzimer, A.; d'Alché-Buc, F.; Fox, E.; and Garnett, R., eds., _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc.
  * Zhang et al. (2020) Zhang, T.; Kishore, V.; Wu, F.; Weinberger, K. Q.; and Artzi, Y. 2020.  BERTScore: Evaluating Text Generation with BERT.  In _International Conference on Learning Representations_.


# Simulations with machine learning potentials identify the ion conduction
mechanism mediating non-Arrhenius behavior in LGPS

Gavin Winter  Massachusetts Institute of Technology, 77 Massachusetts Ave,
Cambridge, MA 02139  Rafael Gómez-Bombarelli  Massachusetts Institute of
Technology, 77 Massachusetts Ave, Cambridge, MA 02139<EMAIL_ADDRESS>
###### Abstract

Li10Ge(PS6)2 (LGPS) is a highly concentrated solid electrolyte, in which
Coulombic repulsion between neighboring cations is hypothesized as the
underlying reason for concerted ion hopping, a mechanism common among
superionic conductors such as Li7La3Zr2O12 (LLZO) and Li1.3Al0.3Ti1.7(PO4)3
(LATP). While first principles simulations using molecular dynamics (MD)
provide insight into the Li+ transport mechanism, historically, there has been
a gap in the temperature ranges studied in simulations and experiments. Here,
we used a neural network (NN) potential trained on density functional theory
(DFT) simulations, to run up to 40-nanosecond long MD simulations at DFT-like
accuracy to characterize the ion conduction mechanisms across a range of
temperatures that includes previous simulations and experimental studies. We
have confirmed a Li+ sublattice phase transition in LGPS around 400 K, below
which the ab-plane diffusivity $D^{*}_{ab}$ is drastically reduced.
Concomitant with the sublattice phase transition near 400 K, there is less
cation-cation (cross) correlation, as characterized by Haven ratios closer to
1, and the vibrations in the system are more harmonic at lower temperature.
Intuitively, at high temperature, the collection of vibrational modes may be
sufficient to drive concerted ion hops. However, near room temperature, the
vibrational modes available may be insufficient to overcome electrostatic
repulsion, thus resulting in less correlated ion motion and comparatively
slower ion conduction. Such phenomena of a sublattice phase transition, below
which concerted hopping plays a less significant role, may be extended to
other highly concentrated solid electrolytes such as LLZO and LATP.

††: J. Phys. Energy

  * November 2022

##  1 Introduction

Materials in the LGPS family exhibit high ionic conductivity, with values
between 5 and 20 mS/cm [1, 2] that are on par with organic liquid electrolytes
in commercial lithium-ion batteries. Since a fully occupied Li+ sublattice
leads to significant electrostatic repulsion between neighboring cations, the
high ionic conductivity of LGPS is thought to arise from inter-cationic
repulsion among lithium ions. In particular, this Coulomb repulsion has been
postulated originally by He et al. to drive the concerted ion hop that boosts
ionic conductivity in these superionic materials systems, namely LGPS, LLZO,
and LATP [3, 4]. LGPS features a BCC-like structure with alternating polyhedra
of PS4 and GeS4 tetrahedra. Unique to crystalline superionic conducting
solids, LGPS exhibits its exceptional ionic conductivity at a stoichiometric
composition, unlike other crystalline solid electrolytes, which rely on
defects or non-stoichiometry to attain high ionic conductivity [2].

The temperature dependence of ion conduction in LGPS has been reported
experimentally, where non-Arrhenius behavior and temperature-dependent
activation energies have been observed. [5, 6, 2, 7]. However, the origin of
the change in activation energy remain a matter of debate. Some groups have
proposed grain boundary effects to explain non-Arrhenius or shifts in
Arrhenius behavior [8], where grain boundaries may dictate ion transport at
low temperatures. However, grain boundary effects have only been shown to be
significant below 250 K in LGPS [8], and at higher temperatures, bulk ion
conduction dominates.

Other groups have hypothesized the existence of a continuous (diffuse) phase
transition [6, 2, 7, 9], characterized as Type II by Boyce and Huberman [10],
only alters the Li+ sublattice. In particular, a transition from quasi one-
dimensional ion conduction at low temperature to three-dimensional ion
conduction at elevated temperatures has been proposed [7, 2, 11, 6, 12, 8,
13]. Nevertheless, ionic conductivity measurements are not available for the
full range of temperatures over which the transition is thought to occur.
Furthermore, in condensed Li+ systems such as LGPS, where the ion conduction
mechanism critically relies on Li+-Li+ interactions, the Haven ratio and other
metrics that quantify the Li+-Li+ interactions (correlation) have not yet been
studied over the range of temperatures from 300 to 800 K to further
characterize the non-Arrhenius behavior. The Haven ratio, as detailed in
Equation 10, can be measured experimentally by back-calculating the total
diffusivity $D_{\sigma}$ from the ionic conductivity and then comparing that
with the tracer diffusivity $D^{tr}$ (equivalent to $D^{*}$). In a review
paper, Kato et al. reported total diffusivities $D_{\sigma}$ and calculated
Haven ratios at room temperature to be less than 1 [2], but this all relied on
data from a single sample by Kuhn et al., the only working group to report
tracer diffusivity $D^{*}$ in LGPS, as characterized by nuclear magnetic
resonance (NMR) [5]. As a result, the values for $D_{\sigma}$ and $D^{*}$ did
not come from the same sample of LGPS. Furthermore, Kuhn et al. themselves
characterized the room-temperature Haven ratio for LGPS as on the order of 1,
by measuring the total diffusivity $D_{\sigma}$ and the tracer diffusivity
$D^{*}$ by NMR for the same sample [5].

Historically, ab initio molecular dynamics (AIMD) has been used to evaluate
the diffusivity and ionic conductivity of crystalline electrolytes. These
simulations require a costly quantum mechanical calculation of the atomic
forces with every femtosecond of simulated time, and consume hundreds of
thousands of CPU hours. Consequently, AIMD is constrained in both the length
and time scales that can be accessed. Thus, simulation of solid ion conductors
with AIMD is limited to high temperatures, where ion hops occur with very high
frequency and diffusivity statistics reach convergence with short simulation
times on the order of 100 ps. The high-temperature regimes in which
diffusivity is evaluated by AIMD often do not overlap with the near-room-
temperature regime of proposed commercial application, where diffusivities and
ionic conductivities are measured experimentally. Comparing the experiment and
simulation requires extrapolating high-temperature AIMD results, under the
assumption that there is no non-Arrhenius behavior.

In the case of LGPS, AIMD and NMR have been used to determine the tracer
diffusivity (self-diffusivity $D^{*}$) in the high- ($>$500 K) and low-
temperature regimes ($<$500 K), respectively. Kato et al. has acknowledged
that no single method has been able to cover the entire range due to the
aforementioned limitations of simulations over the temperature ranges with
respect to the non-Arrhenius behavior [2].

Machine learning provides a way to learn a surrogate for the expensive quantum
mechanical forces, preserving the accuracy but increasing speed by 4 or more
orders of magnitude [14]. As early as 2017, NN potentials have been used to
study ion conduction in crystalline solid electrolytes [15]. Recently, in the
composition space of sulfides, NN potential-driven MD was used by Huang et al.
[16], as well as by Miwa et al. [9], to study LGPS. Both NNMD studies
successfully benchmarked against literature values for high-temperature
diffusivity with AIMD. However, the NNMD simulations by Huang et al.
overestimated the experimental diffusivity near room temperature, and there
were only two datapoints for NNMD simulations near room temperature reported
by Miwa et al., which did not allow for a proper Arrhenius fit. These
discrepancies between MD simulations and experiment warranted a more
comprehensive study of the diffusion and correlation mechanism in LGPS at room
temperature. A NN potential is required to drive MD in order to efficiently
run simulations on the timescales required to study the effect of correlation
on diffusion. Ultimately, since these MD studies are currently limited to bulk
ion conduction (i.e. perfect crystals), the effects of grain boundaries cannot
be assessed —the one limitation of this method.

Using NN potentials, which predict the forces on each atom according to the
local environment at each timestep, we have been able to access timescales on
the order of tens of nanoseconds with MD and achieved quantum chemical
accuracy (comparable to ab initio methods), allowing us to evaluate
diffusivities for each composition of a crystalline material in a faster time
than AIMD. This has allowed us to bridge the gap between studies of ion
diffusion in the high temperature regime with AIMD and experimental studies of
ion diffusion near room temperature.

##  2 Methods

###  2.1 Neural Network Interatomic Potential

The underlying potential used to drive NNMD simulations in this work is a
message-passing neural network (MPNN), which initializes atom-wise features
and builds a graph within a local neighborhood of interacting atoms that
repeatedly exchanges messages through updates of node features. The
polarizable atom interaction neural network (PaiNN) is one example of a
message-passing NN that has rotationally equivariant representations [17] and
is the architecture that we employ as a NN potential for driving MD. Unlike
rotation-invariant representations, which are limited in their ability to
propagate directional information, the rotational equivariant representations
of vector features employed by PaiNN allow directional information to be
preserved, allowing for higher accuracy NN potentials with less DFT training
data. As a message-passing NN, if the cutoff is 5 Å, and there are 3
convolutional layers (the hyperparameters chosen for this study), due to the
updates with message-passing, the actual local environment is (3x5=)15 Å, and
effectively 4-body interactions are learned [18]. Intuitively, there is a
tradeoff between accuracy and computational speed in considering the cutoff
size and how large a local environment one considers an atom/ion interacting
with. Upon inference, the total energy can then be predicted for a specific
input geometry. On top of the output energy from PaiNN, a repulsive energy
term with a $1/r_{ij}^{12}$ dependence similar to the repulsive term
accounting for electron cloud overlap in the Lennard-Jones potential, was
added onto the output energy from PaiNN in order to prevent collisions of
atoms within their own excluded volume and improve the stability of the NN
potential. Auto-differentiating the total energy with respect to atomic
coordinates yields the negative of the per-atom forces as in Equation 1, where
$i$ is the atom index.

$\hat{\vec{F}}_{i}=-\frac{\partial{\hat{E}}}{\partial{\vec{r}_{i}}}$ (1)

In contrast, auto-differentiating the energy with respect to pairwise
distances between all atoms allows the calculation of stress as in Equation 3,
where $a$ and $b$ are running along $x$, $y$, $z$, and $i$ and $j$ are running
along all the atom indices.

$\hat{\vec{f}}_{ij}=\frac{\partial{\hat{E}}}{\partial{\vec{r}_{ij}}}$ (2)
$\hat{\tau}_{ab}=\frac{1}{2\Omega}\sum_{j=1,i\neq
j}(\vec{r}_{ij})_{a}(\hat{\vec{f}}_{ij})_{b}$ (3)

Each timestep in MD requires a call to evaluate forces by this method.

###  2.2 Generation of training data

Generating a high-fidelity NN potential for MD requires training data for
energy, forces, and stress, corresponding to structures in phase space
relevant to the phenomenon being investigated —in this case ion hopping.
Density functional theory (DFT), using the pseudopotential projector-augmented
wave method in the Vienna Ab Initio Simulation Package (VASP) [19, 20, 21],
was used to generate the training data with: optimizations (relaxations) of
periodic structures, nudged elastic band (NEB) calculations, and subsequent
single-point calculations of structures. Each datum in the training data
consists of a structure (atomic positions and lattice), the energy of that
structure, per-atom forces, and a stress tensor as evaluated by DFT. (For
details regarding DFT calculations with VASP, see A.2. For details regarding
dataset availability, see A.3.) While previous sampling strategies have
involved using AIMD as training data [22] or random perturbations to atomic
positions to generate off-equilibrium structures [16], we chose to primarily
rely on NEB as a means of creating an initial set of DFT training data that
spans the trajectory of structures representative of the ion conduction
pathways in a crystal. Using a transition state method such as NEB for
sampling helps acquire thousands of off-equilibrium geometries: (1) with
greater efficiency per ab initio calculation (as compared with sampling by
AIMD, where hops may occur once every several thousand timesteps), (2) in
physically relevant phase space (unlike sampling by random perturbations).

The equilibrium structure that was used for the NEB calculations was the
50-atom unit cell with stoichiometric LGPS mp-696128 [23], and two different
concerted ion hops were sampled with NEB. When there are a multitude of
possible sites in the crystal framework’s void space for the mobile species,
tools that employ Voronoi decomposition [24] or CAVD [25, 26] were used to re-
populate the interstitial sites with Li, and bond valence site energy methods
[27, 28] provided a guess ion migration path. This starting equilibrium
structure has tetragonal symmetry and the zig-zag (Z) arrangement of PS4 /
GeS4 polyhedra. Since only the stoichiometric LGPS compound was used in all
instances, no vacancy concentration was assumed. Training data exclusively
consisted of the 50-atom unit cell and is a sufficient representation for
sampling structures from the ion hopping trajectory. However, the MD was all
run with a 3,200-atom 4x4x4 supercell to mitigate any finite size effects and
provide a wealth of statistics for ion diffusion.

###  2.3 Active learning loops

In order to supplement DFT training data gathered using NEB, a two-pronged
active learning approach [29, 30, 14, 31] was used to further sample frames
from unseen regions in phase space. (1) One approach relied on sampling of
frames from an MD trajectory driven by a first-generation neural network (NN)
potential. First, non-physical frames from the MD trajectory were filtered out
of the pool of structures to sample from. These included: (i) frames with
atoms that collided; (ii) frames with a molten framework since they do not
represent physically-relevant phase space. Though the Li+ sublattice is
expected to be mobile, the criteria for a molten framework was based on the
root-mean-squared-displacement for framework atoms (in reference to the
initial structure) with respect to a threshold value. Subsequently, the
remaining frames in physically-relevant phase space were sampled by their
uncertainty and by their diversity relative to the structures comprising the
existing training dataset. To quantify uncertainty, the variance in forces was
evaluated by an ensemble of NNs. To ensure diversity in sampling, low-
dimensional embedding of the similarity matrix, quantified by root-mean-
squared displacement, between all filtered frames was performed to do k-means
clustering. From each of the 4 clusters, 100 frames were selected for
singlepoint DFT evaluation of forces and energies. (2) The other approach,
developed by Schwalbe-Koda et al. [31], relied on an ensemble of NNs to
identify trajectories that maximize uncertainty (adversarial loss) in the
forces or stress, where uncertainty was again measured as the variance in the
forces or stress evaluated by the ensemble of NNs. In addition to the variance
in forces or stress, the adversarial loss also includes a Boltzmann term
according to the structure’s predicted energy $\hat{E}$ to determine the
frames that were physically accessible in phase space.

With the collection of more training data through these two approaches and
through subsequent generations of training and simulations/sampling, this
active learning process continued until the mean average error of the NN on
the test set with respect to ground-truth DFT fell below 25 meV/Åfor forces
and below 50 meV for energies (as can be seen in Figure 1). While the test
error is a quantitative metric, nanosecond-long stability is the paramount
pre-requisite for MD driven by a NN potential. With this final generation NN
potential (details regarding model availability can be found in A.3),
production MD simulations were then run. MD was run on a 3,200-atom 4x4x4
supercell, using a timestep of 1 fs with the Nosé-Hoover thermostat in the NVT
ensemble and fixed lattice parameters from relaxation at 0 K, as implemented
in the atomic simulation environment (ASE) [32] and adapted for NN potentials.

Figure 1: Comparison between PaiNN-predicted properties and DFT-computed
properties, as evaluated on the test dataset: (1) energy gradient and (2)
total energy. Mean absolute errors (MAE) are below 25 meV/Åand 50 meV for
energies.

##  3 Results

###  3.1 Sublattice phase transition and distinct ion conduction mechanisms in
LGPS: bridging the gap between computation and experiment

After running an MD simulation at each of the fixed temperatures with the NVT
ensemble, the self-diffusivity $D^{*}$ was evaluated using Equation 6.

Figure 2: Top is Arrhenius plot of self-diffusivity $D^{*}$ for Li10Ge(PS6)2,
including AIMD diffusivity measured by Mo et al. [33] and NMR diffusivity
measured experimentally by Kuhn et al. [5] Bottom is the Arrhenius plot
decomposed by 1D ion conduction along the c-channel and 2D ion conduction in
the ab-plane. Lattice parameters used for this simulation are the 0 K lattice
parameters, as relaxed using DFT.

As seen in the top Arrhenius plot in Figure 2, the self-diffusivity $D^{*}$
values from NNMD match closely with both the AIMD values in the high
temperature regime (above 500 K) and with the experimental NMR values in the
near-room temperature regime (below 350 K). While the slope of the lines
(activation energy) is very similar in the two regimes, in between the two
regimes, there is a change in the intercept on the diffusivity axis of nearly
an order of magnitude. This can be attributed to a change in the $\Delta
S_{migration}$ term included in the pre-factor when plotting $\ln(D^{*})$
versus $1/T$. A continuous sublattice phase transition between 350 and 500 K
is hypothesized to be responsible for this change in pre-factor, leaving the
activation energy (migration barrier) essentially unaltered before and after
the transition. From the decomposed ab-plane (2D) diffusivity and c-channel
(1D) diffusivity in the bottom of Figure 2, this sublattice phase transition
may be attributed to a change in the Li+ sublattice, where the ab-plane and
c-channel diffusivity is drastically reduced below 500 K. Interestingly, below
350 K, the decomposed ab-plane diffusivity $D^{*}_{ab}$ (as denoted by the
upside down triangles) continues to follow the same Arrhenius regime as from
350-500 K, while the c-channel diffusivity $D^{*}_{c}$ (as denoted by the
triangles) deviates from the Arrhenius behavior established in the 350-500 K
regime and remains comparatively high.

We also observed a change in the vibrational density of states that supports
the hypothesis of a sublattice phase transition. The vibrational density of
states $D(\nu)$, decomposed by ionic species $i$, can be calculated by
integrating over the Fourier transform of the velocity $\vec{v}_{i}$
autocorrelation function and weighting by mass $m_{i}$, as in Equation 4 [34,
35].

Figure 3: Li+ vibrational density of states spectra, as given by Equation 4,
at each of the temperatures across the sublattice phase transition, from 20 ps
of simulation with a 5 fs sampling frequency. Errors depicted by thickness of
the density of states curve at a given frequency were obtained from variance
over ensemble of MD trajectories, using the same NN potential.

$D(\nu)=\frac{2}{k_{B}T}m_{i}\int<\sum\limits^{N_{i}}_{i=1}\vec{v}_{i}(\tau)\cdot\vec{v}_{i}(t+\tau)>_{\tau}e^{-i2\pi\nu
t}dt$ (4)

Generally, from the temperature-resolved vibrational density of states spectra
in Figure 3, the vibrations in LGPS becomes more harmonic at lower
temperature. For reference, the Li dimer harmonic frequency is around 10 THz
[36], so the peak in the vibrational density of states just below 10 THz may
be attributed to harmonic oscillations of the Li+ sublattice that generally do
not contribute to ion conduction. This peak just below 10 THz increases in
intensity at lower temperature (300 to 400 K), indicating that Li+ ions’
vibrations become more harmonic. The zero-frequency density of states
corresponds with anharmonicity and translational motion, which decrease in
intensity at lower temperatures. The supplemental information also provides
the full vibrational density of states spectra for the other ionic species
(see Figure 7 in A.4), which shows that other framework ions such as Ge, P
also exhibit more sharply peaked (more harmonic) vibrations at lower
temperatures. Notably, there are peaks in the vibrational density of states
spectra for P that coincide with the shoulder around 15 THz in the Li+
vibrational density of states. For this reason, this drop in the shoulder
around 15 THz with decreasing temperature could be attributed in part to the
characteristics of P vibrations, which inhibit Li+ diffusion.

A larger overlap across a range of frequencies is indicative of stronger
coupling and momentum transfer to Li+. Therefore, a broader spectrum for the
vibrational density of states for Li+, such as that of LGPS over all the
observed temperatures, will generally translate to faster ion conduction.
Regardless of the slowdown in the ion conduction mechanism along the ab-plane,
the vibrational density of states still supports the fact that LGPS is an
exceptional ion conductor (even at room temperature). Since integrating over
the vibrational density of states is related to vibrational contributions to
entropy [37], a change in the spectrum for Li+ with respect to temperature may
also support a sublattice phase transition, which is generally manifested in a
change in the entropy of the system and may translate to a change in the
entropy of ion migration pathways.

Figure 4: Radial distribution function (RDF), as given by Equation 5, for Li-P
and Li-Ge as a function of temperature from 300 to 800 K. The bin size is 0.05
Å, and the calculation was taken over an average of 100 frames. Structural
changes are evident with the smearing of the two peaks into one peak for Li-P
with increasing temperature.

The radial distribution functions (RDFs), as given by Equation 5 in Figure 4
indicate that the structure indeed changes as a function of temperature.

$g(\vec{r})=\frac{1}{\rho}\langle\sum_{i\neq
0}(\delta(\vec{r}-\vec{r_{i}}))\rangle$ (5)

Specifically, with the Li-P RDF, as the temperature increases from room
temperature, two distinct peaks in the RDF at 3.0 Åand 3.75 Å, indicative of
two inequivalent sites for Li+ network become smeared into one peak at
elevated temperature.

###  3.2 Correlation in LGPS across temperature regimes

The self-diffusivity $D^{*}$, also known as the tracer diffusivity (plotted in
Figure 2) can be computed by a linear fit of the mean-squared displacement
$<\Delta_{self}(t)>$ over time $t$, as shown in Equation 6. The mean-squared
displacement $<\Delta_{self}(t)>$, or self-correlation, is defined in Equation
8, where $N_{c}$ is the total number of hopping ions $i$.

$D^{*}=\frac{1}{6}\lim_{t\to\infty}\frac{<\Delta_{self}(t)>}{t}$ (6)

The total ionic conductivity $\sigma_{i}$ for a solid electrolyte at
equilibrium (with no applied potential) depends on the total diffusivity
$D_{\sigma}$, which accounts for both the self-diffusivity $D^{*}$ and the
contributions of cross-correlations.

$D_{\sigma}=\frac{1}{6N_{c}}<\absolutevalue{\sum^{N_{c}}_{i=1}[\vec{r}_{i}(t)-\vec{r}_{i}(0)]}^{2}>$
(7)

Thus, the equation for total diffusivity $D$ includes both the self-
correlation (Equation 8) and cross-correlation (Equation 9) contributions
embedded in it.

$\Delta_{self}(t)=\frac{1}{N_{c}}\sum^{N_{c}}_{i=1}\absolutevalue{\vec{r}_{i}(t)-\vec{r}_{i}(0)}^{2}$
(8)
$\Delta_{cross}(t)=\frac{1}{N_{c}}\sum^{N_{c}}_{i=1}\sum^{N_{c}}_{i^{\prime}\neq
i=1}\left[\vec{r}_{i}(t)-\vec{r}_{i}(0)\right]\cdot\left[\vec{r}_{i^{\prime}}(t)-\vec{r}_{i^{\prime}}(0)\right]$
(9)

In single-ion conductors such as crystalline electrolytes, the cross-
correlation $\Delta_{cross}(t)$ would represent the correlation between
cations (intra-species correlation) since anions are immobile. An intuitive
depiction of the pairwise correlation between different ions is provided in
A.1. Without computing the cross-correlation directly, the degree of
correlation in the system can be assessed by comparing the self-diffusivity
$D^{*}$ with the total diffusivity $D_{\sigma}$ via the Haven ratio, since the
self-correlation and cross-correlation terms sum to the total diffusivity
$D_{\sigma}$.

$H_{r}\equiv\frac{D^{*}}{D_{\sigma}}$ (10)

The Haven ratio represents the degree to which cross-correlation contributes
to the total diffusivity (and thus total ionic conductivity). In a perfectly
uncorrelated material system, the Haven ratio would be 1 since the total
diffusivity $D_{\sigma}$ is equal to the self-diffusivity $D^{*}$, as a result
of the cross-correlation being negligible.

Figure 5: Haven ratio across a range of temperatures, relating the self-
diffusivity to the total-diffusivity. A Haven ratio greater than 1 is
indicative of destructive cross-correlation, less than 1 is indicative of
constructive cross-correlation, and equal to 1 corresponds with perfectly
uncorrelated ion motion.

While the self-correlation contributes to the self-diffusivity $D^{*}$ plotted
in Figure 5, constructive correlated motion (cross-correlation) dominates at
higher temperature, where the electrostatic repulsion of neighboring ions aids
in concerted ion hops and boosts total diffusivity relative to the self-
diffusivity. However, at low temperatures (in the region of the continuous
sublattice phase transition), ion motion is less correlated, with Haven ratios
closer to 1 (see Figure 5), consistent with the intuition that electrostatic
repulsion from neighboring cations may impede ion motion in a highly
concentrated system such as LGPS below a certain temperature.

Figure 6: Distinctive van Hove correlation, as given by Equation 11, for Li+
in Li10Ge(PS6)2. The bin size was 0.1 Å, and the sampling frequency was 25 ps.
High intensity region approaching $r\rightarrow 0$ are indicative of
collective ion motion. The bands appearing above 3 Åat 300 K indicate that the
Li+ sublattice is not liquid-like.

The distinctive van Hove correlation function, as defined in Equation 11, can
be interpreted as a time-correlated radial distribution function with respect
to an initial reference Li+ ion [38]. In the following equation, $\delta()$ is
a one-dimensional Dirac delta function, $r$ is the radial distance from the
initial reference Li+ ion, and the sum is performed over all the mobile ions
$N_{c}$, all while being normalized by a factor containing the average number
density of the mobile Li+ ions $\rho$.

$G_{d}(r,t)=\frac{1}{4\pi r^{2}\rho N_{c}n}\langle\sum_{i\neq
j}^{N_{c}}\delta(r-\absolutevalue{\vec{r}_{j}(t+t_{0})-\vec{r}_{i}(t_{0})})\rangle_{t_{0}}$
(11)

The peaks (high intensity regions) in the distinctive van Hove correlation
near $r=0$ (in the heat maps in Figure 6), which are especially evident in the
higher temperature cases after several picoseconds, are indicative of
collective motion. The high intensity region near $r=0$ describes a
neighboring Li+ ion approaching the location of the initial reference Li+ ion
after time $t$. This correlation within 1 Åis likely the result of the
electrostatic repulsion of ions playing a role in the concerted ion motion.
However, the correlation between cations (cross-correlation) takes 20 times
longer to set in at 300 K in comparison to 400 K. Interestingly, in contrast
to most correlation effects, which rely on order, this dynamic correlation
relies on disorder to aid ion conduction and persists at elevated temperature
but is not as evident near room temperature, where part of the Li+ sublattice
is more frozen.

The distinctive van Hove correlation function (Figure 6) can be read as a RDF
by examining one slice in time. From this perspective, it is evident that the
Li+ sublattice is more liquid-like at high temperature with broad peaks in the
RDF (uniform color in the van Hove correlation plot above 3 Å). In contrast,
at room temperature, the peaks in the RDF are narrower, giving rise to bands
appearing around 3 and 5 Åin the van Hove correlation plot, suggesting the Li+
sublattice cannot be characterized as liquid-like at 300 K.

###  3.3 Discussion of discrepancy in sublattice phase transition temperature

We observed a discrepancy between the self-diffusivity $D^{*}$ computed by MD
simulations and that measured by NMR [5] in the temperature range between 350
and 450 K, as seen in Figure 2. Liang et al. reported a phase transition in
the Li4GeS4-Li3PS4 pseudo-binary in the range of 300 K to 400 K, heavily
dependent on composition. Using NMR to analyze the local environment around P
in the Li-Ge-P-S quaternary system, the phase transition was attributed to
changes in the Li sublattice and local environment around P from NMR [39]. The
temperature of the phase transition is a very sensitive function of
composition. As a result, the discrepancy between NNMD and experiment could be
due to a small difference in the composition of LGPS studied experimentally by
Kuhn et al. with NMR (in Figure 2) and the nominal LGPS structure simulated in
this study. Undesired orthorhombic LGPS phase was reported with the high
synthesis temperatures by Kuhn et al. [5]. Similarly, there could be impurity
phases containing the Pa or Pc instead of the Z arrangement of PS4 and GeS4
polyhedra arrangement simulated in this study. Such impurity phases such as
orthorhombic LGPS or $\beta$-Li3PS4, which lower the self-diffusivity of LGPS
(relative to the Z tetragonal LGPS simulated in this study) may have resulted
in an altered sublattice phase transition temperature.

In addition to the discrepancy between the self-diffusivity $D^{*}$ calculated
by this study versus experiment, there is also a discrepancy between the self-
diffusivity $D^{*}$ of LGPS calculated by this MD study and other MD studies
driven by neural network potentials [16, 9]. At room temperature, both the
results of Huang et al. and Miwa et al. overestimate the self-diffusivity
$D^{*}$ as compared to the experimental NMR results (see A.8). The discrepancy
may arise from either the neural network architectures - the deep potential
generation scheme [40] employed by Huang et al. or the architecture employed
by Miwa et al. [9]. Another source of possible disagreement are issues related
to convergence of diffusivity statistics with respect to simulation time or
finite size effects (discussed in A.7). The studies of LGPS by Huang et al.
[16] and Miwa et al. [9] use 900-atom and 800-atom supercells, respectively,
whereas this study uses a 3200-atom supercell. The present study is the only
MD simulation that matched the experimental self-diffusivity $D^{*}$ (tracer
diffusivity) at room temperature, as seen in Figure 2. There is a finite size
effect reported in LLZO, another highly concentrated solid electrolyte, where
smaller supercell sizes overestimate the diffusivity [41], likely because
small supercells artificially induce correlation (order). At high temperature
this effect is not as evident because of the liquid-like Li+ sublattice,
whereas at low temperature it may enforce artificial correlation in the ion
hops, resulting in some spurious superionic behavior. In addition, the
vibrational modes in a smaller supercell may be different from those in the
perfect bulk crystal. Similarly for the LLZO prototype system, there appears
to be similar non-Arrhenius behavior approaching room temperature from
elevated temperature [41, 42, 43]. If MD simulations overestimate the
diffusivity due to a finite size effect, such important non-Arrhenius behavior
may be missed. This non-Arrhenius behavior may be indicative of a sublattice
phase transition between two different ion conduction regimes [42], further
demonstrating the importance of simulation techniques that are capable of
bridging across the varied ion conduction regimes that occur at different
temperatures.

##  4 Outlook and Conclusion

Using a fast and accurate NN potential, MD simulations of ionic transport in
LGPS have been performed for nearly 40 ns. Two distinct ion conduction regimes
have been identified at high temperature and at low temperature, resulting in
agreement with historical experimental NMR data and AIMD simulations. A
sublattice phase transition is hypothesized to be responsible, and was
previously hypothesized but had not been confirmed due to the gap between
simulations performed above 500 K and experimental studies performed below 500
K. This sublattice phase transition temperature can be attributed to a drastic
drop in the self-diffusivity along the ab-plane $D_{ab}^{*}$. Below this
temperature, ion conduction in LGPS can be characterized by a Haven ratio
closer to 1, where Li+ ion motion is less correlated and vibrations are more
harmonic. Both LLZO and LGPS are unique, in that both materials systems are
highly concentrated in Li+, enabling a concerted hopping mechanism due to
electrostatic repulsion at elevated temperature. Therefore, a slowdown in ion
conduction may be a shared characteristic among other superionic conductors
that exhibit a concerted ion hopping mechanism. Below a threshold temperature,
the electrostatic repulsion between neighboring mobile ions cannot be overcome
by the vibrational modes (distribution of Li+ velocities) available, and
consequently ion motion is less concerted at lower temperatures. Despite the
change in ion conduction mechanism approaching room temperature from elevated
temperature, LGPS remains an exceptional ion conductor at room temperature.
More broadly, this study demonstrates the need for fast, accurate methods for
driving MD to study ion conduction in crystalline solids at room temperature,
without relying on extrapolation.

## Acknowledgments

The authors thank Simon Axelrod for implementing PaiNN and Dr. Johannes
Dietschreit for many fruitful discussions regarding MD simulations and
analysis. The authors would also like to acknowledge Dr. Andrey Poletayev for
identifying an error in our original analysis of diffusivities. NNMD was run
using the computing resources of MIT Supercloud [44] and Summit at the
Oakridge Leadership Computing Facility (through grant number ALCF-MAT237)
through the Department of Energy (DOE) Advanced Scientific Computing Research
(ASCR) Leadership Computing Challenge (ALCC). DFT work was run with VASP,
using the computing resources of Perlmutter and Cori KNL at the National
Energy Research Scientific Computing Center (through grant number BES-
ERCAP0021172) through the DOE Office of Science’s Energy Research Computing
Allocations Process. The authors also thank the DOE Advanced Manufacturing
Office (grant number DE-EE0009096) and the Toyota Research Institute (TRI)
through their Accelerated Materials Design and Discovery (AMDD) program for
providing financial support for this work.

## Appendix

###  A.1 Intuitive description of correlation

For the sake of computing the diffusivity by multiplication of matrices or
tensors, as can be done when using Equations 8 and 9, the following matrix
interpretation of ion displacement and diffusion may be helpful. The
displacement can be written as the following vector:

$\vec{R}_{i}(t)=\vec{r}_{i}(t)-\vec{r}_{i}(0)$ (12)

The self-correlation, which is used to calculate the self-diffusivity (when
divided by time), can be intuitively thought of as the on-diagonal elements in
the displacement correlation matrix.

$\blockarray{ccccc}&Li^{+}_{1}Li^{+}_{2}Li^{+}_{i}\\\
\block{c[cccc]}Li^{+}_{1}\vec{R}_{1}\cdot\vec{R}_{1}{\color[rgb]{.75,.75,.75}\vec{R}_{1}\cdot\vec{R}_{2}}{\color[rgb]{.75,.75,.75}\ldots}{\color[rgb]{.75,.75,.75}\vec{R}_{1}\cdot\vec{R}_{i}}\\\
Li^{+}_{2}{\color[rgb]{.75,.75,.75}\vec{R}_{2}\cdot\vec{R}_{1}}\vec{R}_{2}\cdot\vec{R}_{2}{\color[rgb]{.75,.75,.75}\ldots}{\color[rgb]{.75,.75,.75}\vec{R}_{2}\cdot\vec{R}_{i}}\\\
{\color[rgb]{.75,.75,.75}\vdots}{\color[rgb]{.75,.75,.75}\vdots}\ddots{\color[rgb]{.75,.75,.75}\vdots}\\\
Li^{+}_{i}{\color[rgb]{.75,.75,.75}\vec{R}_{i}\cdot\vec{R}_{1}}{\color[rgb]{.75,.75,.75}\vec{R}_{2}\cdot\vec{R}_{2}}{\color[rgb]{.75,.75,.75}\ldots}\vec{R}_{i}\cdot\vec{R}_{i}\\\
$ (13)

The cross-correlation, which is added to the self-correlation to calculate the
total diffusivity, can be intuitively thought of as the off-diagonal elements
in the displacement correlation matrix.

$\blockarray{ccccc}&Li^{+}_{1}Li^{+}_{2}Li^{+}_{i}\\\
\block{c[cccc]}Li^{+}_{1}{\color[rgb]{.75,.75,.75}\vec{R}_{1}\cdot\vec{R}_{1}}\vec{R}_{1}\cdot\vec{R}_{2}\ldots\vec{R}_{1}\cdot\vec{R}_{i}\\\
Li^{+}_{2}\vec{R}_{2}\cdot\vec{R}_{1}{\color[rgb]{.75,.75,.75}\vec{R}_{2}\cdot\vec{R}_{2}}\ldots\vec{R}_{2}\cdot\vec{R}_{i}\\\
\vdots\vdots{\color[rgb]{.75,.75,.75}\ddots}\vdots\\\
Li^{+}_{i}\vec{R}_{i}\cdot\vec{R}_{1}\vec{R}_{2}\cdot\vec{R}_{2}\ldots{\color[rgb]{.75,.75,.75}\vec{R}_{i}\cdot\vec{R}_{i}}\\\
$ (14)

###  A.2 Additional calculation details for VASP DFT

The ab initio calculations were carried out as follows: an energy cutoff of
520 eV was used for the plane-wave basis set. The Perdew-Burke-Ernzerhof (PBE)
generalized gradient approximation (GGA) functional was used [21], but a
functional with a higher level of theory could have been used for training
data calculations, embedding into the (NN) potential a higher accuracy, if
desired. A 2 × 2 × 2 $\Gamma$-centered k-point mesh was used to sample the
Brillouin zone. The energy criterion for self-consistency was when the energy
difference between subsequent electronic steps was less than 1 $\mu$eV. For
NEB and structure optimizations (relaxations), the convergence criterion was
set to an atomic force tolerance of 0.01 eV/Å.

###  A.3 Data and model availability

The dataset used for training models and the trained PyTorch models are
available on figshare here. As noted in the main body of the paper, each datum
of this dataset contains the DFT-computed forces and energy corresponding with
each structure (coordinates and lattice). Currently, the PaiNN model can be
used to run MD (with an ASE interface) through the NeuralForceField GitHub
repository, which is a general framework that can be used for running MD with
NN potentials, while the NequIP model can be run using either an ASE interface
or LAMMPS interface with the patch provided through the NequIP GitHub
repository [22]. The raw MD trajectories are available upon request.

###  A.4 Full vibrational density of states for all ions across temperatures

In Figure 7, the peaks corresponding to the vibrational modes of the framework
ions Ge, P, and S in general shift to lower frequencies with increasing
temperature. Notably, the peaks just above 15 THz for P and S broaden at
elevated temperature and smear out the shoulder that is more distinctly
present at room temperature. Such peak broadening can be attributed to less
harmonic vibrations at elevated temperatures since harmonic vibrations would,
in comparison, appear as narrower delta functions.

Figure 7: Vibrational density of states spectra decomposed by element, as
given by Equation 4, at each of the temperatures across the sublattice phase
transition, from 20 ps of simulation with a 5 fs sampling frequency. Select
peaks at specific frequencies for Ge, P, S become more sharply peaked (higher
intensity) at lower temperatures as compared to higher temperatures,
indicative of greater harmonicity at low temperatures

###  A.5 Thermostat parameters and their effect on diffusivity

The Nosé-Hoover ttime parameter as implemented in ASE for the Arrhenius plot
in Figure 2 was set to 20, where the timestep was set to 1 fs. With the NVE
ensemble (no thermostat), the supercell thermally equilibrated well, and the
sublattice phase transition is still evident in Figure 8, demonstrating that
the thermostat does not interfere with the ion hopping dynamics or the
phenomena investigated here. Each of these NNMD runs lasted 10 ns, which is
why convergence may have not been reached for some of the lower temperatures
in Figure 8. The spacing between temperatures was also not consistent since
the only parameter that can strictly be controlled with NVE is the
initialization temperature, whereas the average temperature is subject to
fluctuations.

Figure 8: Arrhenius plot of self-diffusivity $D^{*}$ for Li10Ge(PS6)2,
including AIMD diffusivity measured by Mo et al. [33] and NMR diffusivity
measured experimentally by Kuhn et al. [5] This is MD run without a thermostat
(i.e. NVE ensemble), where lattice parameters used for this simulation are the
0 K lattice parameters.

###  A.6 Accounting for thermal expansion and its effect on diffusivity

Since 0 K lattice parameters are non-physical, the thermal expansion of the
lattice should be accounted for, using the experimental lattice parameters
from elevated temperature to demonstrate that the sublattice phase transition
still occurred. The results from this are shown in Figure 9. Each of these
NNMD runs lasted 10 ns, which is why convergence may have not been reached for
some of the lower temperatures. Indeed, at low temperature, the diffusivities
in Figure 9 are generally about the same as in Figure 2 since the lattice
parameters at low temperature are close to the lattice parameters at 0 K. In
general there is poorer agreement between NNMD and AIMD at higher temperature,
but this is to be expected, given that the reference AIMD by Mo et al. [33]
was run using 0 K lattice parameters.

Figure 9: Arrhenius plot of self-diffusivity $D^{*}$ for Li10Ge(PS6)2,
including AIMD diffusivity measured by Mo et al. [33] and NMR diffusivity
measured experimentally by Kuhn et al. [5] This is MD run with the Nosé-Hoover
thermostat (i.e. NVT ensemble). Lattice parameters used for this are the
experimental lattice parameters that correspond with each temperature on the
Arrhenius plot, accounting for thermal expansion.

###  A.7 Effect of finite simulation size on diffusivity

As seen in Figure 10, there is good agreement in diffusivity values at
temperatures above 500 K across all supercell sizes. However, below 500 K,
which is the region of interest for the sublattice phase transition, the
smaller 2x2x2 supercell (400 atoms), overestimates the diffusivity. Similarly,
as seen in the next section A.8, comparing our work and Huang et al.’s work on
MD simulations with a 900-atom supercell using a NN potential trained on the
same underlying DFT data (PBE GGA) [16], the diffusivity is overestimated at
lower temperatures when a small supercell is used. In both of these instances,
this may be due to a finite size effect, which Huang et al. also observed
[16]. Such a finite size effect could be attributed to smaller supercell
sizes, artificially inducing correlation (order) in comparison to a true
infinite bulk crystal. For this reason, we opted for a larger 4x4x4 supercell
(3200 atoms) for production MD simulations, which would be able to capture the
phenomena of interest —a change in the ion conduction mechanism.

Figure 10: Arrhenius plot of self-diffusivity $D^{*}$ for Li10Ge(PS6)2,
including AIMD diffusivity measured by Mo et al. [33] and NMR diffusivity
measured experimentally by Kuhn et al. [5] This is MD run with the Nosé-Hoover
thermostat (i.e. NVT ensemble) and varied supercell sizes (2x2x2 = 400 atoms,
4x4x4 = 3200 atoms). The lattice parameters used for this simulation are the 0
K lattice parameters.

###  A.8 Comparison with other machine learning potentials for LGPS

As seen in Figure 11, there is a discrepancy between the diffusivity
calculated by NNMD simulations previously reported by Huang et al. [16] and
experiment [5] approaching room temperature. This may be attributed to the
finite size effect discussed in the previous section A.7.

Figure 11: Arrhenius plot of self-diffusivity $D^{*}$ for Li10Ge(PS6)2,
including NNMD literature values measured by Huang et al. [16] and NMR
literature values measured experimentally by Kuhn et al. [5] This is MD run
with the Nosé-Hoover thermostat (i.e. NVT ensemble). The supercell size used
in our NNMD simulation was 4x4x4 (3200 atoms), while the supercell size used
by NNMD from literature (Huang et al.) was 900 atoms. The lattice parameters
used for this simulation are the 0 K lattice parameters.

## References

## References

  * [1] Kamaya N, Homma K, Yamakawa Y, Hirayama M, Kanno R, Yonemura M, Kamiyama T, Kato Y, Hama S, Kawamoto K and Mitsui A 2011 Nature Materials 10 682–686
  * [2] Kato Y, Hori S and Kanno R 2020 Advanced Energy Materials 10 2002153
  * [3] He X, Zhu Y and Mo Y 2017 Nature Communications 8 15893
  * [4] Xu M, Ding J and Ma E 2012 Applied Physics Letters 101 031901
  * [5] Alexander Kuhn, Viola Duppel and V Lotsch B 2013 Energy and Environmental Science 6 3548–3552
  * [6] Weber D A, Senyshyn A, Weldert K S, Wenzel S, Zhang W, Kaiser R, Berendts S, Janek J and Zeier W G 2016 Chemistry of Materials 28 5905–5915
  * [7] Kwon O, Hirayama M, Suzuki K, Kato Y, Saito T, Yonemura M, Kamiyama T and Kanno R 2014 Journal of Materials Chemistry A 3 438–446
  * [8] Kuhn A, Köhler J and Lotsch B V 2013 Physical Chemistry Chemical Physics 15 11620–11622
  * [9] Miwa K and Asahi R 2021 Solid State Ionics 361 115567
  * [10] Boyce J B and Huberman B A 1979 Physics Reports 51 189–265
  * [11] Schweiger L, Hogrefe K, Gadermaier B, Rupp J L M, Martin H and Wilkening R 2022 Journal of the American Chemical Society 144 9597–9609
  * [12] Iwasaki R, Hori S, Kanno R, Yajima T, Hirai D, Kato Y and Hiroi Z 2019 Chemistry of Materials 31 3694–3699
  * [13] Liang X, Wang L, Jiang Y, Wang J, Luo H, Liu C and Feng J 2015 Chemistry of Materials 27 5503–5510
  * [14] Axelrod S, Schwalbe-Koda D, Mohapatra S, Damewood J, Greenman K P and Gómez-Bombarelli R 2022 Accounts of Materials Research 3 343–357
  * [15] Li W, Ando Y, Minamitani E and Watanabe S 2017 Journal of Chemical Physics 147 214106
  * [16] Huang J, Zhang L, Wang H, Zhao J and Cheng J 2021 Journal of Chemical Physics 154 94703
  * [17] Schütt K T, Unke O T and Gastegger M 2021 arXiv URL http://arxiv.org/abs/2102.03150
  * [18] Batatia I, Batzner S, Kovács D P, Musaelian A, Simm G N C, Drautz R, Ortner C, Kozinsky B and Csányi G 2022 arXiv URL https://arxiv.org/abs/2205.06643v1
  * [19] Kresse G and Hafner J 1993 Physical Review B 48 13115
  * [20] Kresse G and Furthmüller J 1996 Physical Review B 54 11169
  * [21] Perdew J P, Burke K and Ernzerhof M 1996 Physical Review Letters 77 3865
  * [22] Batzner S, Musaelian A, Sun L, Geiger M, Mailoa J P, Kornbluth M, Molinari N, Smidt T E and Kozinsky B 2022 Nature Communications 13 2453
  * [23] Jain A, Ong S P, Hautier G, Chen W, Richards W D, Dacek S, Cholia S, Gunter D, Skinner D, Ceder G and Persson K A 2013 APL Materials 1 011002
  * [24] Rycroft C H 2009 Chaos: An Interdisciplinary Journal of Nonlinear Science 19 041111
  * [25] He X, Bai Q, Liu Y, Nolan A M, Ling C and Mo Y 2019 Advanced Energy Materials 9 1902078
  * [26] He B, Mi P, Ye A, Chi S, Jiao Y, Zhang L, Pu B, Zou Z, Zhang W, Avdeev M, Adams S, Zhao J and Shi S 2021 Acta Materialia 203 116490
  * [27] Chen H, Wong L L and Adams S 2019 Acta crystallographica Section B, Structural science, crystal engineering and materials 75 18–33
  * [28] Wong L L, Phuah K C, Dai R, Chen H, Chew W S and Adams S 2021 Chemistry of Materials 33 625–641
  * [29] Wang W, Yang T, Harris W H and Gómez-Bombarelli R 2020 Chemical Communications 56 8920–8923
  * [30] Ang S J, Wang W, Schwalbe-Koda D, Axelrod S and Gómez-Bombarelli R 2021 Chem 7 738–751
  * [31] Schwalbe-Koda D, Tan A R and Gómez-Bombarelli R 2021 Nature Communications 12 5104
  * [32] Hjorth Larsen A, Jørgen Mortensen J, Blomqvist J, Castelli I E, Christensen R, Dułak M, Friis J, Groves M N, Hammer B, Hargus C, Hermes E D, Jennings P C, Bjerre Jensen P, Kermode J, Kitchin J R, Leonhard Kolsbjerg E, Kubal J, Kaasbjerg K, Lysgaard S, Bergmann Maronsson J, Maxson T, Olsen T, Pastewka L, Peterson A, Rostgaard C, SchiØtz J, Schütt O, Strange M, Thygesen K S, Vegge T, Vilhelmsen L, Walter M, Zeng Z and Jacobsen K W 2017 Journal of Physics: Condensed Matter 29 273002
  * [33] Mo Y, Ong S P and Ceder G 2012 Chemistry of Materials 24 15–17
  * [34] Dietschreit J C, Peters L D, Kussmann J and Ochsenfeld C 2019 Journal of Physical Chemistry A 123 2163–2170
  * [35] Peters L D, Dietschreit J C, Kussmann J and Ochsenfeld C 2019 The Journal of Chemical Physics 150 194111
  * [36] Jacox M E 2009 Journal of Physical and Chemical Reference Data 27 115
  * [37] McQuarrie D 2000 Statistical mechanics (University Science Books) ISBN 1891389157 URL https://books.google.com/books/about/Statistical_Mechanics.html?id=itcpPnDnJM0C
  * [38] Zhu Z, Deng Z, Chu I H, Radhakrishnan B and Ong S P 2017 Ab initio molecular dynamics studies of fast ion conductors Computational Materials System Design (Springer International Publishing) pp 147–168 ISBN 9783319682808 URL https://doi.org/10.1007/978-3-319-68280-8_7
  * [39] Liang X, Jiang Y, Cai W, Wu S, Wu S, Wang L, Lei Z, Lei Z, Chen J, Chen J, Lei Y, Lei Y, Yang L, Yang L and Feng J 2020 ACS Applied Materials and Interfaces 12 27029–27036
  * [40] Zhang L, Lin D Y, Wang H, Car R and Weinan E 2019 Physical Review Materials 3 023804
  * [41] Klenk M J and Lai W 2016 Solid State Ionics 289 143–149
  * [42] Dai J, Chen Q, Glossmann T and Lai W 2019 Computational Materials Science 162 333–339
  * [43] Qi J, Banerjee S, Zuo Y, Chen C, Zhu Z, Likhit H C M and Ong S P 2021 Materials Today Physics 21 100463
  * [44] Reuther A, Kepner J, Byun C, Samsi S, Arcand W, Bestor D, Bergeron B, Gadepally V, Houle M, Hubbell M, Jones M, Klein A, Milechin L, Mullen J, Prout A, Rosa A, Yee C and Michaleas P 2018 2018 IEEE High Performance Extreme Computing Conference, HPEC 2018 URL https://doi.org/10.1109/HPEC.2018.8547629