# Interpreting Sequence-Levenshtein distance for determining error type and
frequency between two embedded sequences of equal length

Robert Logan1†$\ast$, Amy Wangsness Wehe2†, Dori C Woods3, Jon Tilly3,
Konstantin Khrapko3

1Science and Technology Division, Biology and Bioinformatics Department,
Eastern Nazarene College, Quincy, MA 02170

2Health and Natural Sciences Division, Mathematics Department, Fitchburg State
University, Fitchburg, MA 01420-2697

3College of Science, Department of Biology, Northeastern University, 330
Huntington Ave, Boston, MA 02115

†Co-first authors

$\ast$Corresponding author

Robert Logan<EMAIL_ADDRESS>
Abstract

Levenshtein distance is a commonly used edit distance metric, typically
applied in language processing, and to a lesser extent, in molecular biology
analysis. Biological nucleic acid sequences are often embedded in longer
sequences and are subject to insertion and deletion errors that introduce
frameshift during sequencing. These frameshift errors are due to string
context and should not be counted as true biological errors. Sequence-
Levenshtein distance is a modification to Levenshtein distance that is
permissive of frameshift error without additional penalty. However, in a
biological context Levenshtein distance needs to accommodate both frameshift
and weighted errors, which Sequence-Levenshtein distance cannot do. Errors are
weighted when they are associated with a numerical cost that corresponds to
their frequency of appearance. Here, we describe a modification that allows
the use of Levenshtein distance and Sequence-Levenshtein distance to
appropriately accommodate penalty-free frameshift between embedded sequences
and correctly weight specific error types.

Keywords: Levenshtein distance, weighted error, molecular biology, edit
distance, matrix

##  1 Introduction

Levenshtein distance (LD) is a widely used edit distance metric [1]. The LD
algorithm identifies the number of insertions, deletions and substitutions
needed to convert one sequence to another. LD can also assign weights to each
error type. A weighted error has a numerical cost that is inversely associated
with its appearance frequency. Common applications of LD include natural
language processing such as speech recognition, dialect detection, plagiarism
exposure and spell checking  [2] [3] [4] [5]. LD operates under the assumption
of fixed sequence length, with the analytical window frame for computing
distance including the full length of the two sequences. However,
modifications to the LD algorithm allow comparison of embedded sequences,
which experience error-induced frameshift, without any additional distance
cost as described in the Sequence-Levenshtein distance (SLD) modification [6].
However, SLD cannot accommodate weighted errors. The ability to allow for
frameshift without penalty as well as weighted errors has direct relevance in
many molecular biology applications. For example, DNA sequencing platforms can
introduce a characteristic error profile into the nucleic acid sequence such
that certain error types will occur more frequently than others. In these
cases, the more frequently occurring error types should have a smaller error
weight associated with their cost of distance as to reasonably accommodate
them in sequence comparison [7].

Combining the benefits of weighted errors and frameshift accommodated can be
accomplished by interpreting the location and value of the lowest value along
the last column and last row in the completed unweighted LD table. In other
words, by interpreting the SLD position and value on the completed LD matrix,
the error types and frequency of appearance between the sequences can be
determined. This strategy allows for error-specific weights to be added while
also accommodating error-induced frameshift without additional penalty. The
following mathematical conjecture describes the relationship between error
type and frequency and the location of the lowest value along the last column
and last row of the LD matrix. Examples are given for each case, illustrating
how this information can be used for interpreting the error profile between
strings, which can then be used to incorporate weighted LD with frameshift
correction allowance into the sequence analysis.

Let $(a_{i},b_{j})$ be the entry in the $i^{th}$ row and the $j^{th}$ column
of a LD matrix table created by comparing two strings of the same length, $n$.
Then $(a_{n},b_{n})$ will be the lowest right-hand entry in the table,
positioned at the last row and last column. When an entry in an embedded
string of interest is deleted or inserted, there is a frameshift to the left
or to the right, while the window of analysis (the frame) remains the same
since the lengths of the two strings are uniform. Any elements beyond the
string of interest that constitutes its context will then fill in the empty
space in the case of a deletion. In the case of an insertion, elements may be
pushed outside the analytical window upstream or downstream. The position of
each base in the sequence will be denoted as $[S_{1},S_{2},S_{3}\ldots
S_{n-1},S_{n}]$. A LD matrix table is provided in Figure 1 as a reference for
analyzing the changes in the following cases. The green highlighted cells that
appear in the case matrices represent the SLD positions, which are the cells
to be interpreted.

##  2 Results

##  3 Case 1

Given k insertions, l deletions, any number of substitutions, no error at
entry Sn regardless of Sn frameshift, and no insertion(s) between Sn-1 and Sn,
the entry with the lowest value in column n and row n will be the following:

$\left\\{\begin{array}[]{c}\left(a_{n}, b_{n-\left(l-k\right)}\right) if l\geq
k\\\ (a_{n-(k-l)}, b_{n}) if k>l\\\ \end{array}\right.$

Case 1 example

To change “TAGCTAGC" to “TAGTAGCT", the operations include a deletion of “C"
and the addition of the “T" on the 3’ end due to frameshift. The computed LD
matrix between these words and the green highlighted SLD placement and value
is interpreted to reveal a single deletion as shown below in Figure 2.

##  4 Case 2

Given d consecutive deletion-induced, insertion-induced or bona fide
substitution(s) that start at Sn and accumulate upstream and no downstream
frameshift of the substitutions, regardless of error(s) elsewhere, provided an
equal number of insertions or deletions can account for the changes between
the sequences, the entries that share the lowest values in column n and row n
will be in the positions $\left(a_{n}, b_{n-y}\right)$ and $\left(a_{n-y},
b_{n}\right)$ for all $y\in\left\\{0, 1, \ldots, d\right\\}$

Case 2 example

To change “TAGCTAGC" to “TATAGCTA", the operations can include either a
deletion of “GC" and frameshift towards the 5’ end that allows “TA" to enter
the frame effectively making Sn and Sn-1 substituted, or an insertion of “TA"
that pushes “GC" out of the frame on the 3’ end, also making the Sn and Sn-1
substituted. The computed LD matrix between these words and the green
highlighted SLD placement and value is interpreted to reveal either two
deletions or two insertions as shown in Figure 3.

##  5 Case 3

Given d > 1 consecutive deletion-induced, insertion-induced or bona fide
substitution(s) that start at Sn and accumulate upstream and $p $insertion-
induced downstream frameshifts of the substitutions such that at least one
substitution remains at Sn, regardless of error(s) elsewhere, the entries that
share the lowest values in column n and row n will be in positions
$\left(a_{n-(p)}, b_{n}\right)$ and $\left(a_{n-(p+k)}, b_{n}\right) $for
$k\in\left\\{1,2\ldots y\right\\}$. $y $$=$ the number of downstream
substitutions left in the analytical window. If frameshift-inducing deletions
or insertions neighbor the substitutions, behavior will match case 2.

##  6 Case 3 example

To change “TAGCTAGC" to “ATAAGCTG", the operations can include either two
insertions of “A" and one deletion of “A", two insertions of “A" and a
substitution of “A" to “G" or three insertions of “A". The computed LD matrix
between these words and the green highlighted SLD placement and value is
interpreted to reveal either two insertions and one deletion, two insertions
and one substitution, or three insertions as shown in Figure 4.

##  7 Case 4

Given a k $=$ 1 insertion between Sn-1 and Sn and no error at Sn or Sn-1, any
offsetting error(s) elsewhere such that there is no downstream or upstream
frameshift of Sn, the entries that share the lowest values in column n and row
n will be $(a_{n}, b_{n})$, $\left(a_{n}, b_{n-1}\right)$, and $\left(a_{n},
b_{n-2}\right)$

##  8 Case 4 example

To change “TAGCTAGC" to “TAGCAGTC", the operations can include either a
deletion of “T" and an insertion of “T", a deletion of “T" and substitution of
“C" to “T", or a deletion of “T" and “C". The computed LD matrix between these
words and the green highlighted SLD placement and value is interpreted to
reveal either an insertion and deletion, a deletion and a substitution, or two
deletions as shown in Figure 5.

##  9 Case 5

Given l $=$ 1 deletions of Sn-1, no error at Sn or Sn-2 and no insertion
between Sn-2 and Sn-3, and any error(s) elsewhere such that there is no
frameshift of Sn, the entries that share the lowest values in column n and row
n will be $(a_{n}, b_{n})$, $(a_{n-1}, b_{n})$, and $(a_{n-2}, b_{n})$

Case 5 example

To change “TAGCTAGC" to “TAGTCTAC", the operations can include either an
insertion of “T" and a deletion of “G" or an insertion of “T" and a
substitution of “G" to “C" or an insertion of a “T" and a “C". The computed LD
matrix between these words and the green highlighted SLD placement and value
is interpreted to reveal either an insertion and deletion, an insertion and a
substitution, or two insertions as shown in Figure 6.

##  10 Case 6

Given l $=$ 1 deletions of Sn-1 and a substitution error at Sn such that the
substitution doesn’t equal the original Sn-1, and any error(s) elsewhere such
that there is no frameshift of Sn and the sequences don’t match, the entries
that share the lowest values in column n and row n will be $(a_{n-1}, b_{n})$,
and $(a_{n-2}, b_{n})$. If Sn experiences upstream frameshift under these
conditions, the lowest values will follow the rules in conjecture 2. If Sn
experiences downstream frameshift under these conditions, the lowest values
will follow the rules in conjecture 1.

##  11 Case 6 example

To change “TAGCTAGC" to “TAGTCTAT", the operations can include either an
insertion of “T" and a substitution of “G" to “T" or an insertion of a “T" and
another “T". The computed LD matrix between these words and the green
highlighted SLD placement and value is interpreted to reveal either an
insertion and a substitution or two insertions as shown in Figure 7.

##  12 Case 7

Given l > 1 consecutive deletion errors starting at Sn-1, any error(s)
elsewhere except consecutive off-setting insertions between $S_{n-(l+1)}$ and
$S_{n-(l+2)}$ such that there is no frameshift of Sn, the entries that share
the lowest values in column n and row n will be $\left(a_{n-l},b_{n}\right)$
and $\left(a_{n-\left(l+1\right)},b_{n}\right)$.

##  13 Case 7 example

To change “TAGCTAGC" to “ATAAGACC", the minimum number of operations can
include insertions of three “A" letters and a “C" or insertions of three “A"
letters and a substitution of “T" to “C". The computed LD matrix between these
words and the green highlighted SLD placement and value is interpreted to
reveal either an insertion and a substitution or two insertions as shown in
Figure 8.

##  14 Discussion

The position and value of SLD along the last column and last row of the LD
matrix can reveal the error types and appearance frequencies between two
sequences of the same length. Insertions move the SLD from the corner up the
$i$ border (the last column) whereas deletions move the SLD from the corner to
the left along the $j$ border (the last row). An insertion matched with a
deletion does not move the SLD position along the matrix and neither do
substitutions.

Interestingly, an insertion and a deletion pair can produce the same result as
a substitution in a sequence if they occur in the same place. However, an
insertion and a deletion pair are two errors as opposed to a single
substitution error, so distinguishing between these two options matters when
counting errors. If an SLD occurs in the corner of the $i$ and $j$ borders and
the value is 2 or more, it could represent one or more insertion-deletion
pairs, all substitutions or a combination of both between sequences. If an SLD
occurs in the corner of the $i$ and $j$ borders and the value is 1, it
represents a single substitution error. Furthermore, some values of two or
more in the corner of the $i$ and $j$ borders cannot be substitutions, but
rather represent only insertion and deletion pairs, as demonstrated in the
examples for conjectures 4 and 5. Therefore, if there is ambiguity when
interpreting SLD for error type and frequency at the corner of the $i$ and $j$
borders, it can be resolved by applying any known probability of whether a
substitution is expected to occur more, less or the same as an insertion-
deletion pair. Similarly, conjecture 2 describes a scenario where either all
insertions or an equal number of deletions can make one sequence match the
other. In order to make a determination of which of the two error types is
likely responsible for the sequence change for any given analysis, the
expected appearance relationship between them for that specific analysis needs
to be known. An example of using known probabilities to guide decision making
is relying on a specific DNA sequencer’s error hallmark.

Acknowledgements

We wish to thank Boris Yazlovitsky, Greg Shomo, Mariana Levi and the rest of
the Research Computing team at Northeastern University for their support.

Funding

This research was supported by the Eunice Kennedy Shriver National Institute
of Child Health and Human Development (RO1-HD091439 to JLT, DCW and KK).

Declaration of interest

The authors declare that they have no competing interests.

## References

  * [1] Levenshtein V  “Binary codes capable of correcting deletions, insertions, and reversals”  In _Soviet Physics-Doklady_ 10, 1966, pp. 707–710
  * [2] Su Z et al.  “Plagiarism detection using the Levenshtein distance and Smith-Waterman algorithm”  In _3rd International Conference on Innovative Computing Information and Control, ICICIC’08_ , 2008, pp. 0-3  DOI: 10.1109/ICICIC.2008.422
  * [3] Wieling M et al.  “Measuring Foreign Accent Strength in English: Validating Levenshtein Distance as a Measure”  In _Language Dynamics and Change_ 4, 2014, pp. 253–269
  * [4] Yulianto M, Arifudin R and Alamsyah A  “Autocomplete and Spell Checking Levenshtein Distance Algorithm To Getting Text Suggest Error Data Searching In Library”  In _Scientific Journal of Informatics_ 5, 2018, pp. 75  DOI: 10.15294/sji.v5i1.14148
  * [5] Zgank A and Kacic Z  “Predicting the acoustic confusability between words for a speech recognition system using levenshtein distance”  In _Elektronika Ir Elektrotechnika_ 18, 2012, pp. 81–84  DOI: 10.5755/j01.eee.18.8.2628
  * [6] Buschmann T and Bystrykh L  “Levenshtein error-correcting barcodes for multiplexed DNA sequencing”  In _BMC Bioinformatics_ 14, 2013, pp. 1–10  DOI: 10.1186/1471-2105-14-272
  * [7] Logan R et al.  “3GOLD: optimized Levenshtein distance for clustering third-generation sequencing data”  In _BMC Bioinformatics_ 23, 2022, pp. 1–18  DOI: 10.1186/s12859-022-04637-7


# CU-UD: text-mining drug and chemical-protein interactions with ensembles of
BERT-based models

Mehmet Efruz Karabulut1 K. Vijay-Shanker1,∗ Yifan Peng2,∗ ∗co-corresponding
1Computer & Information Sciences, University of Delaware, Newark, DE, USA
2Population Health Sciences, Weill Cornell Medicine, New York, NY, USA

###### Abstract

Identifying the relations between chemicals and proteins is an important text
mining task. BioCreative VII track 1 DrugProt task aims to promote the
development and evaluation of systems that can automatically detect relations
between chemical compounds/drugs and genes/proteins in PubMed abstracts. In
this paper, we describe our submission, which is an ensemble system, including
multiple BERT-based language models. We combine the outputs of individual
models using majority voting and multilayer perceptron. Our system obtained
0.7708 in precision and 0.7770 in recall, for an F1 score of 0.7739,
demonstrating the effectiveness of using ensembles of BERT-based language
models for automatically detecting relations between chemicals and proteins.
Our code is available at https://github.com/bionlplab/drugprot_bcvii.

###### Index Terms:

drug-protein relation, relation extraction; deep learning; ensemble learning

##  I Introduction

Relation extraction is an important task in the biomedical domain that aims to
identify the semantic relationship between the biomedical entities mentioned
in the text. Information about these interactions can be helpful for other
biomedical research such as drug discovery, database curation and question
answering systems. In recent years, there have been lots of efforts to
automatically extract the relations from the biomedical articles (e.g., (1,
2)). BioCreative events have focused on different relation extraction tasks
over the years. One of the tasks in BioCreative VI was ChemProt (3), which
aims to promote the development and evaluation of systems that can
automatically detect relations between chemical compounds and genes/proteins
in PubMed abstracts. This task has been continued as a DrugProt task in
BioCreative VII (4). In this paper, we discuss our participation in this task,
describing ensemble approaches and incorporating multiple BERT-based language
representation models.

##  II Methods

###  II-A Dataset

Organizers of the DrugProt track created a manually annotated corpus of 5,000
abstracts from PubMed (4). The dataset was split into a training set (3,500
abstracts), a development set (750), and a test set (750) (Table I). Following
the settings of previous BioCreative tracks, the test set consists of a large
collection of records containing a subset of a total of 750 Gold Standard
records that was used for evaluation purposes (3).

TABLE I: Statistics of the DrugProt dataset | Training | Development
---|---|---
Documents | 3,500 | 750
Tokens | 1,001,168 | 199,620
Entities | 89,529 | 18,858
Chemical | 46,274 | 9,853
Gene | 43,255 | 9,005
Relations | 17,288 | 3,765
ACTIVATOR | 1,394 | 246
AGONIST | 652 | 131
AGONIST-ACTIVATOR | 27 | 10
AGONIST-INHIBITOR | 10 | 2
ANTAGONIST | 931 | 218
DIRECT-REGULATOR | 2,061 | 458
INDIRECT-DOWNREGULATOR | 1,297 | 332
INDIRECT-UPREGULATOR | 1,351 | 302
INHIBITOR | 5,277 | 1,152
PART-OF | 861 | 258
PRODUCT-OF | 904 | 158
SUBSTRATE | 1,921 | 495
SUBSTRATE_PRODUCT-OF | 25 | 3

###  II-B Preprocessing

Since all relations in the DrugProt corpus are enclosed in one sentence, our
system considers interactions only occur within a sentence. We split the raw
text into individual sentences by Stanza sentence splitter (5). These
sentences are passed onto the BERT framework which tokenizes them using
Wordpiece tokenizer (6).

Given a sentence with multiple entity mentions marked, we treated the relation
extraction problem as a multi-class classification problem. Curators of the
DrugProt corpus annotated thirteen key relations of biomedical importance. A
sentence instance with a marked protein and a chemical mention is considered
one of these thirteen DrugProt relation classes if there is a relationship
between a protein and chemical occurrences. Otherwise, this instance is
negative. There may be multiple interacting entities and non-interacting
entities marked in the same sentence. We consider only a single chemical-
protein pair in that sentence per sentence instance. Therefore, we formulate
the task as predicting one of the fourteen relation types for a sentence
instance.

We consider two entity tagging schemes (Table II). In the first scheme, we
anonymized target named entities in a sentence using pre-defined tags to
prevent overfitting to entities (7). Specifically, given that each instance
was comprised of a sentence and two entity mentions, we replaced each entity
of interest with the tokens DRUG or PROTEIN, depending on what the specific
mention represented. Other drugs and proteins are replaced with DRUG_O and
PROTEIN_O, respectively.

TABLE II: Examples of two entity tagging schemes Original sentence |  human type 12 RDH reduces dihydrotestosterone to androstanediol
---|---
Schema 1 |  DRUG reduces dihydrotestosterone to PROTEIN
Schema 2 |  $<$DRUG-B$>$ human type 12 RDH $<$DRUG-E$>$ reduces dihydrotestosterone to $<$PROTEIN-B$>$ androstanediol $<$PROTEIN-E$>$

In the second scheme, we directly inserted the start and end entity markers
into the original sentence. For example, we used the format $<$PROTEIN-B$>$
human type 12 RDH $<$PROTEIN-E$>$ to enclose the protein mention “human type
12 RDH” in the sentence.

###  II-C Model Development

We hypothesize that an ensemble system that combines results of different
models could lead to better predictive performance than using a single model
(1). The superior performance of ensemble models is achieved when there is
enough diversity among the individual models (8). In this study, we propose
two different ensemble approaches which incorporate multiple BERT-based
language representation models. Here, we chose BERT as the backbone because
recent studies demonstrate it achieved state-of-the-art performance on a
number of natural language understanding tasks, including relation extraction
task (9, 10, 7).

We explore various architectures and pre-trained models to achieve the goal of
diversifying BERT-based models. We will discuss these techniques and the
ensemble algorithms in the following subsections. Due to the space
limitations, we explain only the most relevant parts of the BERT model to the
DrugProt task in the following. We refer the reader to the original paper for
the details.

####  II-C1 Single models

We built three different architectures. The first one uses the last layer of
the [CLS] token for the relation extraction task because it aggregates the
information of the whole sentence. Fig. 1a demonstrates the fine-tuning
process of BERT on the relation extraction task. A sentence instance is
tokenized and embedded as in the original BERT paper. Inspired by the work of
Su et al. (11), we also explore two other methods of summarizing the
information in the last layer: Long Short Term Memory (LSTM) on the last layer
(Fig. 1b) and attention mechanism on all units of the last layer (Fig. 1c).

(a)

(b)

(c)

Figure 1: Architectures of individual models on relation extraction.

To obtain diversity among the models used in our ensembles, we chose three
BERT-based models in our approach, BioBERT (12), PubMedBERT (13), and BioM-
ELECTRA (14, 15), since they differ in some pre-training aspects as well as
masking, vocabulary, and the learning objectives techniques. The pre-training
of BioBERT used PubMed and PMC after the initial pre-training of BERT, whereas
PubMedBERT pre-training using PubMed alone started from scratch. These two
models are BERT-Base with 12 bidirectional self-attention heads. To study the
effectiveness of large models for this task, we also used BioM-ELECTRAL. BioM-
ELECTRAL is pre-trained on PubMed abstracts with no previous pre-training and
is based on the ELECTRA architecture (14).

####  II-C2 Class-weighted loss function

Since we have created a negative instance for each non-interacted chemical and
protein in the same sentence, we have substantially more negative examples in
the training dataset. To mitigate this unbalance, we assigned class weights
inversely proportional to their respective frequencies.

####  II-C3 Majority voting

Majority voting is the simplest yet effective weighting method for output
fusion. We select the relation type that gets the most votes out of five
individual models. If the models cannot agree on the majority vote, we assign
this instance as negative.

####  II-C4 Meta-learning method

We also employ another ensemble method using multilayer perceptron (MLP) as a
meta-learner to learn from individual classifiers. Specifically, we use MLP to
combine the outputs of each individual model in an ensemble (Fig. 2). After
each individual model is fine-tuned, we extract the last output layer of the
[CLS] token from each model. We then use the development set that we reserve
for the ensemble learning to obtain the activation vector for the [CLS] token.
Afterward, we concatenate them into one vector and pass it to another MLP. We
insert a softmax layer on top of the MLP to obtain final predictions.

Figure 2: The architecture of the MLP-based ensemble model. $h^{i}_{[CLS]}$
indicates the [CLS] token output of the final layer for the $i$’th model in
the ensemble.

##  III Results and Discussion

###  III-A Experimental settings

We submitted five runs. For Runs 1 and 2, we used the five individual models
listed in Group A (PubMedBERT, PubMedBERT-CW, PubMedBERT-LSTM, BioBERT-
Attention, Biobert-Tag2) and Group B (PubMedBERT, PubMedBERT-CW, PubMedBERT-
LSTM, BioM-ELECTRAL, BioBERT-Tag2). “Tag2” denotes the second entity tagging
scheme. “CW” indicates that this model was fine-tuned using class-weighted
loss function.

Individual model predictions were combined using a stacking approach involving
a MLP with one hidden layer having 512 units as described in Section II-C4.
Run 3 involves the same Group B and uses both stacking and the majority
voting. Runs 4 and 5 apply only one type of model, BioM-ELECTRAL. Five models
are obtained for Run 4 by training on five different partitions. However, Run
5 involves only one model, allowing us to see the extent of the gain by using
an ensemble method.

While participant teams were provided with the training data set and the
development data set, we combined them together to develop the models. For
Runs 1 and 2, we first used 70% of the total data for training, 20% for
development (tuning the hyper-parameters for the individual models) and 10% to
train the ensemble models. For the final models, we trained the individual
models on the original training data and ensemble models on the development
data with the chosen hyper-parameters.

For Run 3, we used the same approach of 70-20-10 split, but repeated this
process five times (Fig 3). Thus, we obtained a group of five models
corresponding to each 70-20-10 partition. Next, we trained a separate MLP for
each group of models using the 10% as described above. As a result, we have
five predictions for each instance corresponding to each MLP. Finally, we
combined them using majority voting to obtain the final predictions.

Figure 3: 70-20-10 split is performed 5 times.

For Run 4, we fine-tuned five BioM-ELECTRAL models using 80% of total data in
each fold and then combined their predictions using majority voting.

Run 5 involved only one model, and hence 100% of the data was used for its
fine-tuning.

Our models are implemented in Tensorflow 1.14 and all the experiments are
carried out on NVIDIA 2080Ti GPUs. We employ AdamW optimizer with an initial
learning rate of 2e-5. The maximum epoch is set as 10.

###  III-B Results and discussions

Table III and IV shows the results we obtained for 5 Runs on the development
set in iteration 1 (Fig 3) and the official test set, respectively. Runs 3 and
4 have the best F1 performance, which is substantially better than the average
F1 across participant teams. Interestingly, F1 score of the single large model
is almost the same as those of the ensembles of Runs 1 and 2. The best F1
scores of Runs 3 and 4 are due to the increase of Precision even if the Recall
dropped slightly.

TABLE III: Results of our systems on the development set in iteration 1 (Fig 3) Run | System | P | R | F1
---|---|---|---|---
1 | Stacking (MLP) | 0.7630 | 0.7720 | 0.7674
2 | Stacking (MLP) | 0.7700 | 0.7838 | 0.7764
3 | Stacking (MLP)+Majority Voting | 0.7770 | 0.7780 | 0.7770
4 | Majority Voting | 0.7622 | 0.7832 | 0.7726
5 | BioM-ELECTRAL | - | - | -
TABLE IV: Results of our systems on DrugProt test set Run | System | P | R | F1
---|---|---|---|---
1 | Stacking (MLP) | 0.7421 | 0.7902 | 0.7654
2 | Stacking (MLP) | 0.7360 | 0.7925 | 0.7632
3 | Stacking (MLP)+Majority Voting | 0.7708 | 0.7770 | 0.7739
4 | Majority Voting | 0.7721 | 0.7750 | 0.7736
5 | BioM-ELECTRAL | 0.7548 | 0.7747 | 0.7647
| Mean of participant teams | 0.6430 | 0.6291 | 0.6196

Table V shows the results of individual relation types using our best system
(Run 3). Our system obtained the best performance of 0.8902 in F1 for the
ANTAGONIST relation, and 0.8570 for the INHIBITOR relation. As shown in Table
I, INHIBITOR is the most frequent relation in the dataset. Thus its
classification is probably well trained. On the other hand, ANTAGONIST is one
of the least frequent relations. Further error analysis is needed to
understand why our model can effectively extract it with such few cases.
Compared to results of individual relation types by other Runs, we observed
that Runs 3 and 4 have a uniform increase in F1 performance for almost all
individual relation types. On the other hand, Run 1 has the best F1
performance on the AGONIST type with a 2.5% increase over Run 3.

TABLE V: Results of individual relation types using our best system Relation | P | R | F1
---|---|---|---
ACTIVATOR | 0.8083 | 0.8083 | 0.8083
AGONIST | 0.7857 | 0.7623 | 0.7738
ANTAGONIST | 0.8342 | 0.9542 | 0.8902
DIRECT-REGULATOR | 0.7579 | 0.6643 | 0.7080
INDIRECT-DOWNREGULATOR | 0.6976 | 0.7894 | 0.7407
INDIRECT-UPREGULATOR | 0.7260 | 0.7942 | 0.7586
INHIBITOR | 0.8586 | 0.8553 | 0.8570
PART-OF | 0.6781 | 0.7763 | 0.7239
PRODUCT-OF | 0.6285 | 0.7292 | 0.6751
SUBSTRATE | 0.7189 | 0.6348 | 0.6742

##  IV Conclusion

In this manuscript, we describe our submission in the BioCreative VII DrugProt
task. The results demonstrate that our ensemble system can effectively detect
the chemical-protein relations from biomedical literature.

## Acknowledgment

This work is supported by the National Library of Medicine under Award No.
4R00LM013001.

## References

  * (1) Peng, Y., Rios, A., Kavuluru, R. and Lu, Z., 2018. Extracting chemical–protein relations with ensembles of SVM and deep learning models. Database, 2018.
  * (2) Zhang, Y., Lin, H., Yang, Z., Wang, J., Zhang, S., Sun, Y. and Yang, L., 2018. A hybrid model based on neural networks for biomedical relation extraction. Journal of biomedical informatics, 81, pp.83-92.
  * (3) Krallinger, M., Rabal, O., Akhondi, S.A., Pérez, M.P., Santamaría, J., Rodríguez, G.P., Tsatsaronis, G. and Intxaurrondo, A., 2017, October. Overview of the BioCreative VI chemical-protein interaction Track. In Proceedings of the sixth BioCreative challenge evaluation workshop (Vol. 1, pp. 141-146).
  * (4) Antonio Miranda, Farrokh Mehryary, Jouni Luoma, Sampo Pyysalo, Alfonso Valencia and Martin Krallinger. Overview of DrugProt BioCreative VII track: quality evaluation and large scale text mining of drug-gene/protein relations. Proceedings of the seventh BioCreative challenge evaluation workshop. 2021.
  * (5) Qi, P., Zhang, Y., Zhang, Y., Bolton, J. and Manning, C.D., 2020, July. Stanza: A Python Natural Language Processing Toolkit for Many Human Languages. In ACL: System Demonstrations (pp. 101-108).
  * (6) Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, et al. 2016. Google’s neural machine translation system: Bridging the gap between human and machine translation. arXiv preprint arXiv:1609.08144.
  * (7) Peng, Y., Yan, S. and Lu, Z., 2019, August. Transfer Learning in Biomedical Natural Language Processing: An Evaluation of BERT and ELMo on Ten Benchmarking Datasets. In Proceedings of the 18th BioNLP Workshop and Shared Task (pp. 58-65).
  * (8) Sagi, O. and Rokach, L., 2018. Ensemble learning: A survey. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 8(4), p.e1249.
  * (9) Devlin, J., Chang, M.W., Lee, K. and Toutanova, K., 2018. Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805.
  * (10) Su, P., Peng, Y. and Vijay-Shanker, K., 2021, June. Improving BERT Model Using Contrastive Learning for Biomedical Relation Extraction. In Proceedings of the 20th Workshop on Biomedical Language Processing (pp. 1-10).
  * (11) Su, P. and Vijay-Shanker, K., 2020, December. Investigation of bert model on biomedical relation extraction based on revised fine-tuning mechanism. In 2020 IEEE International Conference on Bioinformatics and Biomedicine (BIBM) (pp. 2522-2529). IEEE.
  * (12) Lee, J., Yoon, W., Kim, S., Kim, D., Kim, S., So, C.H. and Kang, J., 2020. BioBERT: a pre-trained biomedical language representation model for biomedical text mining. Bioinformatics, 36(4), pp.1234-1240.
  * (13) Gu, Y., Tinn, R., Cheng, H., Lucas, M., Usuyama, N., Liu, X., Naumann, T., Gao, J. and Poon, H., 2020. Domain-specific language model pretraining for biomedical natural language processing. arXiv preprint arXiv:2007.15779.
  * (14) Alrowili, S. and Vijay-Shanker, K., 2021, June. BioM-Transformers: Building Large Biomedical Language Models with BERT, ALBERT and ELECTRA. In Proceedings of the 20th Workshop on Biomedical Language Processing (pp. 221-227).
  * (15) Clark, K., Luong, M. T., Le, Q. V., and Manning, C. D. 2020. Electra: Pre-training text encoders as discriminators rather than generators. arXiv preprint arXiv:2003.10555.


|  Multi-reference coupled cluster theory using the normal ordered exponential
ansatz
---|---
|  Alexander Gunasekera,1∗ Nicholas Lee,1∗ and David P. Tew1∗
|  Properly spin-adapted coupled-cluster theory for general open-shell
configurations remains an elusive goal in electronic structure theory. In this
contribution we examine Lindgren’s normal-ordered exponential ansatz using
spin-free excitation operators, with the aid of automatic equation generation
software. We present a size-extensive reformulation of the unlinked working
equations, and analyse the performance of the method with single and double
excitations for simple molecular systems in terms of accuracy and size
consistency.

##

††footnotetext: 1 Physical and Theoretical Chemistry Laboratory, University of
Oxford, South Parks Road, Oxford, OX1 3QZ, UK††footnotetext: ∗ Corresponding
author email<EMAIL_ADDRESS>∗
Corresponding author email<EMAIL_ADDRESS>∗
Corresponding author email<EMAIL_ADDRESS>
##  1 Introduction

Single-reference coupled-cluster theory1 is firmly established and is widely
used in high-accuracy electronic structure calculations. Multi-reference
coupled-cluster (MRCC) theories, in contrast, remain an area of active
research due to challenges arising from the complexity of defining wave
operators and working equations for correlating open-shell electrons. Multi-
determinant reference states are needed to describe open-shell configurations
which are important in many chemical systems including in diradicals,
transition metal, and molecules at non-equilibrium geometries. Several MRCC
ansätze and approaches have been forwarded and explored for correlating open-
shell configurations2, 3, 4, but no consensus method has yet emerged. The
ansätze can be broadly classified as either a Jeziorski-Monkhorst (JM) ansatz5
or an internally-contracted ansatz6, 7, 8. In the former, a set of excitation
operators are used for each determinant of the reference wavefunction. In the
latter, a single exponential operator is applied to the whole reference state.

In this work we focus on the internally-contracted ansatz. The usual choice,
taken in the the Unitary Group Approach to Coupled Cluster (UGA-CC)9, 10, 11,
12, 13 and in the internally contracted MRCC (ic-MRCC)6, 7, 8 methods, is to
employ a Complete Active Space (CAS) reference wavefunction. Active-to-active
excitations are then omitted from the cluster operator, which ensures that the
Baker-Campbell-Hausdorff (BCH) expansion terminates in the energy and residual
equations.

Many open-shell systems, such as organic radicals and transition metal spin
states, can be well approximated through an appropriately chosen single
configuration state function (CSF), which is an eigenfunction of the spin
operators. In these cases, it is not necessary to perform a full CAS
calculation, whose cost formally scales exponentially with the number of open-
shell electrons. In this work, we explore fully spin-adapted internally-
contracted coupled cluster theory for correlating single configurational state
functions (CSF). These are multi-determinant states with coefficients
determined by spin and spatial symmetry constraints. They are invariant with
respect to rotations among the doubly-occupied orbitals, but not in general
among the open-shell (active) orbitals.

To fully account for correlation-induced orbital relaxation of the CSF
reference, the cluster-operator must contain active-to-active excitations. We
therefore adopt Lindgren’s normal-ordered exponential (NOE) form of the wave
operator14. This choice ensures that there are no contractions among cluster
operators and thus the working equations terminate at finite order in the
cluster amplitudes. The NOE ansatz has previously been considered in the
context of Fock-Space MRCC15, 16, 17, 18, 19, 20 and the Similarity
Transformed Equation-of-Motion Coupled Cluster21, 22, 23, 24, 25, 26 (ST-EOM-
CC) approaches. In these two methods, normal-ordering of the exponential
allows for a partial decoupling of the different valence sectors, simplifying
the working equations. A variation on the NOE ansatz, where only contractions
through spectator orbitals are allowed, have also been considered by Mukherjee
and co-workers in a Combinatoric open-shell coupled cluster theory (COS-CC27,
28, 29).

In this work we propose and examine through numerical testing a rigorously
spin-adapted NOE coupled-cluster approach for correlating arbitrary spin-
states. Spin-adaptation is achieved through the use of spin-free excitation
operators in our cluster operator, in a similar manner to that of UGA-CC
methods. Due to the complexity of obtaining spin-free working equations with
the NOE ansatz, we have implemented our own domain-specific Wicks contraction
engine and automated equation generator, using similar ideas to those of
Hermann and Hanrath30, 31, 32.

##  2 Theory

###  2.1 Preliminaries

We employ the following convention for orbital labelling: $i,j,k,\dots$ for
core (doubly-occupied) orbitals, $a,b,c,\dots$ for virtual (unoccupied)
orbitals, $t,u,v,\dots$ for active orbitals, and $p,q,r,\dots$ for general
orbitals. The core, active, and virtual spaces will be denoted $\mathbb{C}$,
$\mathbb{A}$, and $\mathbb{V}$ respectively.

The CSF, $\Ket{\Phi_{0}}$, is expressed through a linear combination of
creation operator strings acting upon a closed shell reference,
$\ket{\text{cs}}$.33 That is,
$\Ket{\Phi_{0}}=\hat{O}^{S,M}_{N}(\textbf{t})\ket{\text{cs}}$ where
$\hat{O}^{S,M}_{N}(\textbf{t})$ is specified by a genealogical coupling vector
t to give spin $S$ and spin projection $M$.

$\hat{O}^{S,M}_{N}(\textbf{t})=C^{S,M}_{t_{N},\frac{1}{2}}\hat{O}^{S-t_{N},M-\frac{1}{2}}_{N-1}(\textbf{t})a^{\dagger}_{N\alpha}+C^{S,M}_{t_{N},-\frac{1}{2}}\hat{O}^{S-t_{N},M+\frac{1}{2}}_{N-1}(\textbf{t})a^{\dagger}_{N\beta}$
(1)

$C^{S,M}_{t_{N},\pm\frac{1}{2}}$ are the Clebsch-Gordon coefficients. For
example, the operator used to form an open-shell singlet from
$\ket{\text{cs}}$ is given by

$\textsection\textsection\textsection\textsection\textsection\textsection\textsection\textsection\textsection\textsection\textsection\textsection\hat{O}^{0,0}_{2}=\frac{1}{\sqrt{2}}\Big{(}a^{\dagger}_{p\alpha}a^{\dagger}_{q\beta}-a^{\dagger}_{p\beta}a^{\dagger}_{q\alpha}\Big{)}$
(2)

and the $M=0$ component of a triplet by

$\hat{O}^{1,0}_{2}=\frac{1}{\sqrt{2}}\Big{(}a^{\dagger}_{p\alpha}a^{\dagger}_{q\beta}+a^{\dagger}_{p\beta}a^{\dagger}_{q\alpha}\Big{)}$
(3)

We work with pseudo-creation ($\hat{a}^{p\sigma}$) and pseudo-annihilation
($\hat{a}_{p\sigma}$) operators defined with respect to the _closed-shell
vacuum_ , $\ket{\text{cs}}$. They are related to creation and annihilation
operators in the following way34:

$\begin{split}\hat{a}^{i\sigma}=\hat{i\sigma}\quad\quad\hat{a}_{i\sigma}=\hat{i\sigma}^{\dagger}\quad\quad\hat{a}^{a\sigma}=\hat{a\sigma}^{\dagger}\\\
\quad\quad\hat{a}_{a\sigma}=\hat{a\sigma}\quad\quad\hat{a}^{u\sigma}=\hat{u\sigma}^{\dagger}\quad\quad\hat{a}_{u\sigma}=\hat{u\sigma}\end{split}$

such that
$\hat{a}_{i\sigma}\Ket{\Phi_{0}}=\hat{a}_{u\sigma}\Ket{\Phi_{0}}=\hat{a}_{a\sigma}\Ket{\Phi_{0}}=0$,
and $\hat{i\sigma}^{\dagger}$, $\hat{u\sigma}^{\dagger}$ and
$\hat{a\sigma}^{\dagger}$ are creation operators acting upon the core, active
and virtual spin-orbitals, respectively, and $\sigma$ is the spin index. We
also employ the generalised tensor notation of Head-Gordon and White35,
expressing pseudo-creation and pseudo-annihilation operators as
$\hat{a}^{p\sigma}=\hat{a}_{p\sigma}^{\dagger}$ and
$\hat{a}_{p\sigma}=\hat{a}_{p\sigma}$, respectively. These operators are also
used by Nooijen and Bartlett21, 22 in the context of ST-EOM-CC. We work with
n-body spin-free excitation operators $\hat{T}_{n}$ each of which commute with
the spin operator $\hat{S}^{2}$.

$\hat{T}=\hat{T}_{1}+\hat{T}_{2}+\dots=\sum_{\begin{subarray}{c}p\in\mathbb{A}\cup\mathbb{V}\\\
q\in\mathbb{C}\cup\mathbb{A}\end{subarray}}t_{p}^{q}\,\hat{E}_{q}^{p}+\frac{1}{2}\sum_{\begin{subarray}{c}p,q\in\mathbb{A}\cup\mathbb{V}\\\
r,s\in\mathbb{C}\cup\mathbb{A}\end{subarray}}t_{pq}^{rs}\,\hat{E}_{rs}^{pq}+\dots$
(4)

$\displaystyle\hat{E}_{q}^{p}$
$\displaystyle=\sum_{\sigma}\hat{a}^{p\sigma}\hat{a}_{q\sigma}$ (5)
$\displaystyle\hat{E}_{rs}^{pq}$
$\displaystyle=\sum_{\sigma\tau}\hat{a}^{p\sigma}\hat{a}^{q\tau}\hat{a}_{s\tau}\hat{a}_{r\sigma}=\hat{E}_{r}^{p}\hat{E}_{s}^{q}$
(6)

The normal-ordered36 spin-free molecular electronic Hamiltonian $\hat{H}$ is
given by:

$\hat{H}=\sum_{pq}h_{p}^{q}\hat{E}_{q}^{p}+\frac{1}{2}\sum_{pqrs}g_{pq}^{rs}(\hat{E}^{p}_{r}\hat{E}^{q}_{s}-\delta_{r}^{q}\hat{E}_{s}^{p})$
(7)

where $h_{p}^{q}=\braket{p}{\hat{h}}{q}$ and
$g_{pq}^{rs}=\braket{pq}{r_{12}^{-1}}{rs}=(pr|r_{12}^{-1}|qs)$.

###  2.2 Normal Ordered Exponential Coupled Cluster

Lindgren’s general normal-ordered coupled-cluster wavefunction is defined as

$\Ket{\Psi}=\left\\{\mathrm{e}^{\hat{T}}\right\\}\Ket{\Phi_{0}}$ (8)

where the reference state $\Ket{\Phi_{0}}$ is an arbitrary CSF and braces
indicate normal-ordering with respect to the closed-shell vacuum. The normal-
ordered exponential form properly parameterises independent correlation
processes and excludes contractions between cluster operators. Without normal-
ordering, contractions among active-to-active excitations in the cluster
operator, such as $\hat{T}_{t}^{u}$ and $\hat{T}_{tx}^{uy}$, would give rise
to non-terminating working equations. Normal-ordering of the exponential is
therefore necessary to ensure that the energy and amplitude equations truncate
at finite order in the cluster operators. The order at which the equations
truncate is $4+\min(n,2N_{r})$, where $n$ is the number of active orbitals and
$N_{r}$ is the maximum rank of the excitation operators. When applied to a
closed-shell reference,
$\left\\{\mathrm{e}^{\hat{T}}\right\\}=\mathrm{e}^{\hat{T}}$ and NOECC is
equivalent to the standard coupled-cluster ansatz.

In SRCC theory, the linked formulation is used. The BCH expansion terminates
at finite order and the energy and amplitude contributions are term-wise size
extensive, allowing for truncated CC theories that retain size-extensivity.
ic-MRCC is also formulated with linked equations, which is possible because
active-to-active exctiations are excluded. However, using a linked formalism
for NOECC theory requires knowledge of the wave operator’s inverse
$\\{e^{\hat{T}}\\}^{-1}$. The construction of this inverse is non-trivial22.
Our approach is to instead use the unlinked form, but to reformulate the
equations to recover size extensivity at every order of expansion in the
cluster operators.

The unlinked energy and amplitude equations are given by

$\displaystyle E$
$\displaystyle=\braket{\Phi_{0}}{\hat{H}\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}$
(9) $\displaystyle R_{I}$
$\displaystyle=\braket{\Phi_{I}}{(\hat{H}-E)\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}$
(10)

where the first-order interacting space is
$\langle\Phi_{I}|=\langle\Phi_{0}|\hat{E}_{I}^{\dagger}$, with $i$ standing
for one of the cluster excitations. Naive application of the truncated
unlinked formalism results in non-size extensive equations. The unlinked
residual equations can be rewritten as

$\displaystyle R_{I}$
$\displaystyle=\braket{\Phi_{I}}{\hat{H}\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}-\braket{\Phi_{I}}{\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}E$
(11)
$\displaystyle=\braket{\Phi_{I}}{\hat{H}\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}-\braket{\Phi_{I}}{\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{H}\left\\{\mathrm{e}^{\hat{T}}\right\\}}{\Phi_{0}}$
(12)

Eq. 12 is equivalent to a connected form of the residual equations. The proof
is provided in Appendix B. Retaining all terms in Eq. 12 up to a given power
of $\hat{T}$ therefore results in size-extensive truncated NOECC working
equations. Specifically, the linearised equations (l-NOECC) are:

$\displaystyle E$
$\displaystyle=\braket{\Phi_{0}}{\hat{H}\left\\{1+\hat{T}\right\\}}{\Phi_{0}}$
(13) $\displaystyle R_{I}$
$\displaystyle=\braket{\Phi_{I}}{\hat{H}\left\\{1+\hat{T}\right\\}}{\Phi_{0}}-\braket{\Phi_{I}}{\left\\{\hat{T}\right\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{H}}{\Phi_{0}}$
(14)

and the quadratic equations (q-NOECC) are:

$\displaystyle E=$
$\displaystyle\braket{\Phi_{0}}{\hat{H}\left\\{1+\hat{T}+\frac{1}{2}\hat{T}^{2}\right\\}}{\Phi_{0}}$
(15) $\displaystyle R_{I}=$
$\displaystyle\braket{\Phi_{I}}{\hat{H}\left\\{1+\hat{T}+\frac{1}{2}\hat{T}^{2}\right\\}}{\Phi_{0}}-\braket{\Phi_{I}}{\left\\{\hat{T}\right\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{H}\left\\{1+\hat{T}\right\\}}{\Phi_{0}}$
$\displaystyle-\braket{\Phi_{I}}{\left\\{\frac{1}{2}\hat{T}^{2}\right\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{H}}{\Phi_{0}}$
(16)

Since standard CCD only includes terms up to quadratic in the amplitudes, our
truncated q-NOECCD method exactly replicates closed-shell CCD, and q-NOECCSD
differs from CCSD only by neglecting third- and fourth-order terms involving
singles amplitudes. A similar NOECC wavefunction is used in the UGA-SSMRCC
ansatz, but with multiple CSF references rather than one.37, 38 However, the
UGA-SSMRCC method uses somewhat different working equations.

###  2.3 Choice of Cluster Operator

The choice of which excitations to include in the cluster operator is
motivated by the desire that it should generate the entire first-order
interacting space out of the reference state. For excitations from core
orbitals to virtual, this leads to the usual inclusion of singles and doubles
as in standard CCSD approaches. We also include core-to active, active-to-
virtual, and active-to-active excitations in both the singles and doubles. The
active-to-active excitations are important to capture correlation-induced
orbital relaxation of the reference state. The Goldstone diagrams for the
general excitation operators for NOECCSD are depicted in Figure 1. For
reference states with only one active orbital, we exclude purely active-to-
active excitations in the singles and doubles.

$t_{a}^{i}\hat{E}_{i}^{a}$$t_{ab}^{ij}\hat{E}_{ij}^{ab}$$t_{u}^{i}\hat{E}_{i}^{u}$$t_{a}^{u}\hat{E}_{u}^{a}$$t_{v}^{u}\hat{E}_{u}^{v}$$t_{ua}^{ij}\hat{E}_{ij}^{ua}$$t_{uv}^{ij}\hat{E}_{ij}^{uv}$$t_{ab}^{iu}\hat{E}_{iu}^{ab}$$t_{av}^{iu}\hat{E}_{iu}^{av}$$t_{va}^{iu}\hat{E}_{iu}^{va}$$t_{vw}^{iu}\hat{E}_{iu}^{vw}$$t_{ab}^{uv}\hat{E}_{uv}^{ab}$$t_{wa}^{uv}\hat{E}_{uv}^{wa}$$t_{wx}^{uv}\hat{E}_{uv}^{wx}$

Fig. 1: Goldstone diagrams for the singles and doubles amplitudes used in our
NOECCSD. Core and virtual electrons are represented by solid arrows; active
electrons are represented by skeleton arrows.

The active-to-active amplitudes introduce spectator excitations, where a
double excitation involves an electron being annihilated and created in the
same active orbital and acts as a single excitation when applied to the
reference. The same residual equation is therefore obtained by projection with
either the single excitation $\hat{E}_{p}^{q}$ or the double excitation with a
spectator $\hat{E}_{pu}^{qu}$, and there are an insufficient number of
equations to uniquely determine the amplitudes39, 5, 15, 40, 41. In the
current formulation, the issue of redundancy is unresolved. However, in the
linearised approximation l-NOECCSD, the redundant excitations are
interchangeable and the infinite set of amplitudes that solve the residual
equations all result in exactly the same energy. We have also observed that
different amplitude solutions also yield exactly the same energies for the
q-NOECCSD method. We find that is always possible to converge the residual
equations, but this typically requires use of an appropriate level shift and
the Direct Inversion in the Iterative Subspace (DIIS) technique.

##  3 Method

###  3.1 Automated equation generator

The large numbers of terms appearing in the spin-free open-shell coupled
cluster equations make it necessary to automate the generation of working
equations. This is done using an object-oriented Python code that implements a
second quantised operator algebra along with a tensor algebra to represent the
associated coefficients. In this framework, a particular coupled cluster
ansatz may be represented as a sum of tensor products, with associated second-
quantised operator products. The open-shell reference is included through the
application of second quantised operators with specific orbital indices that
create the appropriate configuration state function out of the closed shell
Fermi vacuum. Rigorous spin adaptation of the theory is performed at the level
of the coupled-cluster ansatz by including all terms in the spin summation.
This brute force approach leads to very many terms in the equations, but
produces a fully spin-adapted coupled cluster theory. The equations are
projected onto the excited state manifold generated by the cluster operators
and Wick’s theorem is applied to obtain the terms that contribute to the
energy and amplitude equations.

$a^{t\alpha}$$t_{ef}^{lv}\left\\{a^{e\alpha}a^{f\alpha}a_{v\alpha}a_{l\alpha}\right\\}$$g_{jc}^{dk}\left\\{a^{j\alpha}a^{c\alpha}a_{k\alpha}a_{d\alpha}\right\\}$$\Phi_{iu}^{ab}\left\\{a^{i\alpha}a^{u\alpha}a_{b\alpha}a_{a\alpha}\right\\}$$a_{t\alpha}$
$\Rightarrow$ $\Rightarrow$

Fig. 2: Goldstone diagram construction representing the automated generation
of the residual term
$r^{iu}_{ab}=g_{jc}^{dk}\,t_{ef}^{lv}\,\delta^{u}_{t}\delta^{i}_{k}\delta^{f}_{b}\delta^{c}_{a}\delta^{j}_{l}\delta^{e}_{d}\delta^{t}_{v}$

As an example, consider a doublet with one open-shell electron, with spin
$\alpha$, in an orbital labelled $t$. Each contribution $r_{I}$ to a residual
$R_{I}$ is found by evaluating its contribution to the expectation value when
projected onto the excitation manifold $\hat{\Phi}_{I}={\phi}_{I}\hat{E}_{I}$,

$\Braket{a_{t\alpha}}{\hat{\Phi}^{I}r_{I}}{a^{t\alpha}}$ (17)

For example, the linear contribution $r^{iu}_{ab}$ of amplitude $t^{iu}_{ab}$
with $g^{bj}_{ia}$ to the residual $R^{iu}_{ab}$ is found from

$\Braket{a_{t\alpha}}{\phi_{iu}^{ab}\,\left\\{\hat{E}^{iu}_{ba}\right\\}g_{jc}^{dk}\,\left\\{\hat{E}^{jc}_{kd}\right\\}t_{ef}^{lv}\,\left\\{\hat{E}^{ef}_{vl}\right\\}}{a^{t\alpha}}$
(18)

One of the spin cases is

$\Braket{a_{t\alpha}}{\phi_{iu}^{ab}\,\left\\{a^{i\alpha}a^{u\alpha}a_{b\alpha}a_{a\alpha}\right\\}g_{jc}^{dk}\,\left\\{a^{j\alpha}a^{c\alpha}a_{k\alpha}a_{d\alpha}\right\\}t_{ef}^{lv}\,\left\\{a^{e\alpha}a^{f\alpha}a_{v\alpha}a_{l\alpha}\right\\}}{a^{t\alpha}}$
(19)

where possible full contraction is

$\phi_{iu}^{ab}\,r^{iu}_{ab}=\phi_{iu}^{ab}\,g_{jc}^{dk}\,t_{ef}^{lv}\,\delta^{u}_{t}\delta^{i}_{k}\delta^{f}_{b}\delta^{c}_{a}\delta^{j}_{l}\delta^{e}_{d}\delta^{t}_{v}$
(20)

for a residual contribution of

$r^{iu}_{ab}=g_{jc}^{dk}\,t_{ef}^{lv}\,\delta^{u}_{t}\delta^{i}_{k}\delta^{f}_{b}\delta^{c}_{a}\delta^{j}_{l}\delta^{e}_{d}\delta^{t}_{v}$
(21)

The automated procedure to generate this expression is equivalent to the
Goldstone diagram construction in figure 2.

This and other possible spin components are accounted for automatically,
together with the factor of 4 associated with this diagram. In our pilot
implementation, we do not make any attempt to identify equivalent terms in the
equations by the topologies of their diagrammatic construction. Furthermore,
the equations are not factorized in this naive implementation, giving a
suboptimal scaling of computational cost that remains to be addressed in
future versions of the code. It is expected that the efficiency will be
improved by a substantial factor when these steps are introduced through
future code developments.

###  3.2 Solution of CC Equations

The solution scheme used for these NOECC equations is based on the iterative
quasi-Newton method commonly employed in coupled cluster methods. With no true
spin-free Fock matrix available, the preconditioner was the generalized Fock
matrix, given by

$\displaystyle F^{m}_{n}$
$\displaystyle=\sum_{q}\braket{\Phi_{0}}{\hat{E}^{m}_{q}}{\Phi_{0}}h_{n}^{q}+\sum_{qrs}\braket{\Phi_{0}}{\hat{E}^{mr}_{qs}}{\Phi_{0}}g_{nr}^{qs}$
(22)

Iterative solution of the generated equations is aided by use of a level
shift, chosen to be larger than the difference between the Fock matrix
eigenvalues corresponding to the highest core orbital and lowest virtual
orbital. Direct Inversion in the Iterative Subspace (DIIS) was used to
accelerate convergence42, 43.

##  4 Results

###  4.1 Doublets

We start with doublet electronic states because the single determinant ROHF
reference can be correlated using standard unrestricted and spin-restricted
spin-orbital coupled cluster implementations, which enables us to compare
directly to our spin-free approach. In Table 1 we present results for
unrestricted CCSD (uCCSD), Szalay and Gauss’s spin-restricted CCSD (rCCSD)44
and our normal-ordered exponential NOECCSD, truncated to linear and quadratic
order. Szalay and Gauss’s spin-restricted approach constrains the cluster
amplitudes by imposing the exact $S^{2}$ expectation value, rather than fully
spin-adapting the amplitudes. We also present results for the quadratic spin-
free CCSD without normal-ordering. The geometries were obtained from Szalay
and Gauss44 and we use the same cc-pVTZ basis as that work. In all cases, the
ROHF state is used as the reference determinant. Energies were converged to
$10^{-10}$ $\mathrm{E_{H}}$.

Molecule | uCCSD | rCCSD | l-CCSD | q-NOECCSD | q-CCSD
---|---|---|---|---|---
$\mathrm{{}^{2}CH}$ | $-141.105$ | $-140.973$ | $-137.013$ | $-140.851$ | $-141.049$
$\mathrm{{}^{2}OH}$ | $-230.302$ | $-230.156$ | $-223.133$ | $-230.306$ | $-230.331$
$\mathrm{{}^{2}CN}$ | $-354.724$ | $-353.995$ | $-325.369$ | $-353.749$ | $-354.454$
$\mathrm{{}^{2}NO}$ | $-433.834$ | $-433.354$ | $-402.668$ | $-433.327$ | $-434.265$
$\mathrm{{}^{2}CH_{3}}$ | $-199.002$ | $-198.832$ | $-191.233$ | $-198.833$ | $-198.974$
$\mathrm{{}^{2}NH_{2}}$ | $-220.008$ | $-219.831$ | $-211.911$ | $-219.806$ | $-219.991$

Table 1: Correlation energies in $\mathrm{mE_{H}}$

The rCCSD correlation energies are smaller than uCCSD due to the spin-
restriction imposed on the amplitudes. The spin-free NOECCSD energies are are
in close agreement with rCCSD, differing by only 0.3 m$\mathrm{E_{H}}$. Exact
agreement between spin-restricted and spin-adapted theories is not expected,
even if we had not truncated the NOECCSD equations to quadratic order. The
difference between the two approaches small in comparison to the expected
magnitude of the contribution from the three-body excitations that are absent
in both methods. The non normal-ordered spin-free CCSD correlation energies
are universally lower than for NOECCSD, in some cases by as much as 0.8
m$\mathrm{E_{H}}$. This is because, without normal-ordering, quadratic terms
involving active orbitals result in suprious additional energy contributions.
For example $\hat{t}_{a}^{u}t_{u}^{i}E_{u}^{a}E_{i}^{u}$ gives rise to an
additional core to virtual excitation $t_{a}^{u}t_{u}^{i}E_{i}^{a}$ that
lowers the energy. Normal-ordering has no impact on the linear terms, but
linearised NOECCSD only recovers 90% of the correlation energy.

###  4.2 Open Shell Singlets and Triplets

####  4.2.1 Beryllium Atom

Having demonstrated our spin-free NOECC framework for one-electron doublets,
we now turn to the two-electron cases: the triplet and the open-shell singlet.
The ${}^{3}P$ and ${}^{1}P$ 2s2p excited states of the beryllium atom
exemplify these electronic configurations, in a system simple enough to allow
comparison of the singlet-triplet splitting with full configuration
interaction (FCI).

Table 2 lists computed CSF reference energies, and all-electron FCI and
NOECCSD energies of the singlet ${}^{1}P$ and triplet ${}^{3}P$ states. uCCSD
energies are also included for the triplet state, which has a single reference
determinant. We use the cc-pCVTZ basis and the orbitals for both CSF
references were generated from an ROHF calculation for the $M=1$ triplet
state. The CSF for the open-shell singlet is that generated by Eq. 2. Energies
are converged to $10^{-10}$ $\mathrm{E_{H}}$.

The uCCSD, q-NOECCSD and FCI energies for the single-reference triplet state
are all in excellent agreement, within 60 $\mu\mathrm{E_{H}}$. The agreement
between CCSD and FCI is a consequence of the higher-order excitations all
involving core-valence correlaiton, which is small for Be. We can also
conclude that, one hand, spin-contamination in the uCCSD wavefunction for this
state is very small, and on the other hand, the neglected cubic and higher
NOECCSD contributions must also either be small in magnitude, or cancel.

The singlet CSF reference is almost 0.1$\mathrm{E_{H}}$ above FCI, whereas the
triplet reference was 0.05$\mathrm{E_{H}}$ above the corresponding FCI energy.
This is in part a consequence of using the the ROHF orbitals optimised for the
triplet, but also because the correlation among singlet spin-coupled electrons
is larger than triplet spin-coupled due to the Fermi heap in their joint
probability distribution.45 The effect of both orbital relaxation and
correlation are well captured in the spin-free NOECCSD approach. At linear
order, l-NOECCSD is within 10 m$\mathrm{E_{H}}$ of FCI, which reduces to 0.3
m$\mathrm{E_{H}}$ for q-NOECCSD. The singlet-triplet gap is also accurate to
0.3 m$\mathrm{E_{H}}$ compared to FCI.

$2S+1$ | CSF | uCCSD | l-CCSD | q-NOECCSD | FCI
---|---|---|---|---|---
$3$ | $-14.513129$ | $-14.561352$ | $-14.561504$ | $-14.561358$ | $-14.561412$
$1$ | $-14.356425$ | $-$ | $-14.473941$ | $-14.462727$ | $-14.463057$
$\Delta$ | $0.156703$ | $-$ | $0.087563$ | $0.098630$ | $0.098355$

Table 2: $2s2p$ excited state energies and singlet-triplet gap $\Delta$ of the
beryllium atom in $\mathrm{E_{H}}$

A fully spin-adapted theory should recover the same energy for each of the
$(2S+1)$ spin-projections of a state with total spin $S$. Our approach enables
us to correlate arbitrary CSFs, and as a test of our method, we computed all-
electron l-NOECCD energies of the $M=1$ and $M=0$ components of the triplet
state. The $M=1$ CSF is a single ROHF determinant, whereas the $M=0$ CSF is a
linear combination of two determinants (Eq. 3). The two calculations converged
with energies at every iteration that were identical, confirming that our
method is fully spin-adapted.

####  4.2.2 Oxygen Molecule

Another archetypical singlet-triplet system is provided by the oxygen
molecule. The lowest energy singlet state ${}^{1}\Delta_{g}$ lies 0.0359
$\mathrm{E_{H}}$ above the triplet ${}^{3}\Sigma_{g}^{-}$ ground state. Both
states correspond to the open-shell configuration
$\pi_{g,x}^{1}\pi_{g,y}^{1}$. The $O_{2}^{1,1}$ reference CSF for the
${}^{3}\Sigma_{g}^{-}$ is a single ROHF determinant, whereas the $O_{2}^{0,0}$
reference CSF for the ${}^{1}\Delta_{g}$ state is a linear combination of two
determinants. Table 3 report all-electron NEOCCSD energies for the singlet and
triplet states of $\mathrm{O_{2}}$ at a bond length of 1.2075Å using a cc-
pCVTZ basis. For both CSFs, the orbitals were obtained from a ROHF calculation
on the triplet state. We also report energies for the ${}^{3}\Sigma_{g}^{-}$
ground state computed by ROHF-uCCSD for comparison. All energies are converged
to $10^{-10}$ $\mathrm{E_{H}}$.

${}^{2S+1}\Lambda$ | Reference Energy | uCCSD | l-CCSD | q-NOECCSD
---|---|---|---|---
${}^{3}\Sigma_{g}^{-}$ | $-149.653208$ | $-150.223429$ | $-150.238029$ | $-150.222443$
${}^{1}\Delta_{g}$ | $-149.605647$ | $-$ | $-150.204043$ | $-150.183136$
${}^{1}\Delta_{g}-{}^{3}\Sigma_{g}^{-}$ | $0.047561$ | $-$ | $0.033986$ | $0.039307$

Table 3: Energies (in $\mathrm{E_{H}}$) of the triplet ground state and
lowest-lying open and closed shell singlet excited states of $\mathrm{O_{2}}$
in a cc-pCVTZ basis set

As before, the q-NOECCSD energy is within 1 $\mathrm{E_{H}}$ of the uCCSD
energy, where they can be compared. The computed singlet-triplet gap of 0.0393
$\mathrm{E_{H}}$ is in good agreement with the experimental value of 0.0359
$\mathrm{E_{H}}$, considering that the experimental value is the adiabatic 0-0
energy, whereas we have computed the vertical electronic energy, and that
triple excitations are not fully accounted for in our method. Although not
done here, our approach can straightforwardly be extended to treat the
${}^{1}\Sigma_{g}^{+}$. This is a $\pi_{g,x}^{2}+\pi_{g,y}^{2}$ configuration,
which also cannot be treated through uCCSD approaches.

##  5 Discussion

We discuss the characteristics of the NOECC method using criteria set out by
Köhn and coworkers for MRCC theory: size extensivity; size consistency;
orbital invariance; compatibility with SRCC; satisfying the proper residual
condition.

###  5.1 Size extensivity

The NOECC equations are formulated such that, at each order in the amplitudes,
the unlinked terms in the residual equations exactly cancel. As shown in
Appendix B, solving for the amplitude equations is equivalent to solving for a
corresponding set of equations comprising only connected terms. Since the
energy of the CSF reference is also size extensive, our NOECC method is fully
size extensive.

###  5.2 Size consistency

Size consistency requires that the energy of a molecule composed of two
infinitely separated spin-coupled fragments is equal to the sum of the
energies of the isolated fragments. For size consistency to be satisfied, the
reference CSF must dissociate into the corresponding CSFs for the constituent
fragments, and non-vanishing higher-order excitations present in the NOECC
wave operator for the molecule must be representable through products of
lower-order excitations of the fragments.

The appropriate CSF for spin-coupled, but non-interacting fragments is the one
constructed from Clebsch–Gordon coupling of the spins of the fragments. For
example, a two-electron open-shell triplet and singlet are formed by high-spin
or low-spin coupling two one-electron doublet states, respectively. To assess
the size consistency of the NOECC ansatz, we examine the energies for the
homolytic bond fission of Li2 as the simplest non-trivial test case.

At dissociation, the ${}^{1}\Sigma_{g}^{+}$ and ${}^{3}\Sigma_{u}^{+}$ states
of lithium dimer should have the same energy, and this should be equal to
twice that of an isolated lithium atom in the ${}^{2}S$ state. The results of
all-electron spin-free NOECCD and NOECCSD calculations assessing the extent to
which this is true are presented in Tables 4 and 5, respectively. The orbitals
of the ${}^{3}\Sigma_{u}^{+}$ state of Li2 and the ${}^{2}S$ state of the Li
atom were obtained from ROHF calculations. The ${}^{1}\Sigma_{g}^{+}$ state
was constructed using the localised $1s$ and $2s$ orbitals of the
${}^{3}\Sigma_{u}^{+}$ state, ensuring that the energy expectation value of
the open-shell singlet CSF exactly matches the ROHF energy of the triplet, and
that both are exactly twice that of the doublet. The cc-pCVTZ basis was used
and an interatomic distance of $10^{9}$Å was chosen to represent the system of
the two dissociated, but still spin-coupled, lithium atoms. All energies are
converged to $10^{-10}$ $\mathrm{E_{H}}$. uCCD and uCCSD energies are also
provided for comparison.

| uCCD | l-CCD | q-NOECCD
---|---|---|---
$\mathrm{Li}\,({}^{2}S)$ | $-41.505947$ | $-41.680942$ | $-41.532939$
$\mathrm{Li_{2}}\,({}^{3}\Sigma_{g}^{-}),r=10^{9}$Å | $-83.011893$ | $-83.389282$ | $-83.091473$
$\mathrm{Li_{2}}\,({}^{1}\Sigma_{g}^{+}),r=10^{9}$Å | $-$ | $-83.389289$ | $-83.091788$
$2E\left({}^{2}S\right)-E\left({}^{3}\Sigma_{g}^{-}\right)$ | $0.000000$ | $0.027398$ | $0.025595$
$2E\left({}^{2}S\right)-E\left({}^{1}\Sigma_{g}^{+}\right)$ | $-$ | $0.027404$ | $0.025910$

Table 4: CCD correlation energies in $\mathrm{mE_{H}}$ for $\mathrm{Li}$ and
$\mathrm{Li_{2}}$ at a separation of $10^{9}$Å in cc-pCVTZ basis.

| uCCSD | l-CCSD | q-NOECCSD
---|---|---|---
$\mathrm{Li}\,({}^{2}S)$ | $-41.546366$ | $-41.694641$ | $-41.545784$
$\mathrm{Li_{2}}\,({}^{3}\Sigma_{g}^{-}),r=10^{9}$Å | $-83.092731$ | $-83.389282$ | $-83.091498$
$\mathrm{Li_{2}}\,({}^{1}\Sigma_{g}^{+}),r=10^{9}$Å | $-$ | $-83.389289$ | $-83.091947$
$2E\left({}^{2}S\right)-E\left({}^{3}\Sigma_{g}^{-}\right)$ | $0.000000$ | $0.000001$ | $-0.000070$
$2E\left({}^{2}S\right)-E\left({}^{1}\Sigma_{g}^{+}\right)$ | $-$ | $0.000007$ | $0.000379$

Table 5: CCSD correlation energies in $\mathrm{mE_{H}}$ for $\mathrm{Li}$ and
$\mathrm{Li_{2}}$ at a separation of $10^{9}$Å in cc-pCVTZ basis.

Both l-CCD and q-NOECCD show significant size inconsistency errors of $0.03$
$\mathrm{mE_{H}}$. Size consistency is violated because the cluster operator
includes terms that excite fragment A with spectator orbitals on fragment B,
and vice versa

$\hat{E}_{pu}^{qu}\rightarrow\hat{E}_{p_{A}u_{A}}^{q_{A}u_{A}}+\hat{E}_{p_{A}u_{B}}^{q_{A}u_{B}}+\hat{E}_{p_{B}u_{A}}^{q_{B}u_{A}}+\hat{E}_{p_{B}u_{B}}^{q_{B}u_{B}}$
(23)

The terms
$\hat{E}_{p_{A}u_{B}}^{q_{A}u_{B}}+\hat{E}_{p_{B}u_{A}}^{q_{B}u_{A}}$ are non-
vanishing in the molecule, but are absent in the cluster operators of
fragments and the cluster operator is therefore not additively separable in
the dissociative limit, which for CCD would require

$\hat{E}_{pu}^{qu}\rightarrow\hat{E}_{p_{A}u_{A}}^{q_{A}u_{A}}+\hat{E}_{p_{B}u_{B}}^{q_{B}u_{B}}$
(24)

The doubles spectator excitations correspond to single excitations when acting
on the reference. Since the l-CCSD and q-NOECCSD methods explicitly include
the single excitations $\hat{E}_{p_{A}}^{q_{A}}$, $\hat{E}_{p_{B}}^{q_{B}}$
that were absent in NOECCD, size consistency can be restored numerically
through the combined action of the singles and doubles. The size inconsistency
errors are reduced to 7 $\mathrm{nE_{H}}$ for l-CCSD and 4 $\mu\mathrm{E_{H}}$
for 1-NOECCSD. While the NOECCSD method is not rigorously size consistent, it
is numerically very close to being size consistent in practice. Similar
observations have been made by Hanauer and Köhn8 for the ic-MRCC method, where
they found that the magnitude of size inconsistency errors depended critically
on the way they removed the redundancies from their working equations.

###  5.3 Orbital invariance

The NOECC method proposed is invariant to core-core and virtual-virtual
rotations in the same way as standard coupled-cluster theory. The reference
CSF is not in general invariant to active-active rotations. This is in
contrast to a CASSCF reference function, used in ic-MRCC methods, which is
invariant to active-active rotations46. For any given excitation order, the
NOE cluster operator $\hat{T}$ consists of all possible excitations from the
set of core and active orbitals to the set of active and virtual orbitals. The
NOE wave operator $\\{e^{\hat{T}}\\}$ is therefore invariant to orbital
rotations within each of these sets. If a certain subset of these excitations
were excluded, the wave operator ceases to be orbital invariant. Our proposed
NOECC method includes all excitations of a given type and therefore retains
all orbital invariance properties of the reference CSF.

###  5.4 Compatibility with SRCC

The state-specific NOECC method we have proposed is a single-reference method
in the sense that we correlate a single reference CSF representative of the
electronic eigenstate under consideration. In general, the CSF is a specific
linear combination of many determinants. In the case that the CSF is a single
closed-shell determinant, NOECC reduces identically to conventional spin-free
coupled-cluster theory, truncated at the same order in $\hat{T}$ making it
seamlessly compatible with SRCC methods. Many of the problems associated with
combining results from MRCC methods with SRCC methods stem from the use of a
CASSCF reference, where it is often challenging to keep the reference space
consistent across different regions of the potential energy surface. In the
context of state-specific NOECC, this problem becomes one of choosing an
appropriate CSF reference. Some regions of the potential energy surface will
not be well represented by a single CSF47 and there the state-specific
approach is expected to be less accurate.

###  5.5 Projected Schrödinger equation

A proper residual equation48, 3, 2 is defined as one that is equivalent to
solving a projected Schrödinger equation. In general, the amplitude equations
for methods using the JM ansätze, where each reference determinant has its own
wave operator, do not correspond to a projected Schrödinger equation (a
notable exception being Hanrath’s MRexpT49, 50 method). The NOECC approach
uses an internally-contracted ansatz with a complete first-order interacting
space and therefore does satisfy the condition of having a proper residual
equation.

##  6 Conclusion

We have developed a novel coupled-cluster method to correlate multi-
determinant states using a normal-ordered exponential ansatz, NOECC. The
ansatz is rigorously spin-adapted, opening up a new possibility of recovering
the dynamic correlation of an arbitrary configuration state function, without
spin contamination. In contrast to many existing methods, both high- and low-
spin states of an atom or molecule can be correlated.

Our working equations are derived from a reformulation of the unlinked
coupled-cluster equations, which we prove are order-by-order equivalent to
solving fully connected equations and therefore remain size extensive when
truncated to low orders in the cluster operators. In this way, we circumvent
the requirement for the inverse of the normal ordered exponential, for which
the closed form is not known. The method formally terminates at
$4+\min(n,2N_{r})$ in the cluster amplitudes, where $n$ is the number of open-
shell orbitals and $N_{r}$ is the maximum excitation rank of the cluster
operator, and we have developed code to generate the spin-adapted equations
automatically.

The NOECCSD method truncated at second-order has been examined numerically
using a set of small atoms and molecules, with encouraging results. Our
energies for doublet systems are comparable to those found by Szalay and Gauss
using a spin-restricted formalism and the singlet-triplet energy splitting are
shown to be in excellent agreement with FCI for the 1s2s configuration of
berylium and within 10 kJ/mol of experiment for the oxygen molecule. Numerical
tests of size consistency reveal that, while the method is not rigorously size
consistent, size inconsistency errors are on the order of $\mu\mathrm{E_{H}}$
for the cases tested. In common with many MRCC methods, the NOECC wavefunction
contains spectator excitations that lead to a set of redundant amplitudes.
Although this leads to an infinite family of solutions to our residual
equations, we find that different amplitude solutions yield exactly the same
energies.

Our method is suitable for correlating specific spin states, such as those in
organic radicals and transition metal spin states. Since only one CSF is used
as a reference, it cannot be used to dissociate bonds in the same way as
CASSCF-based MRCC methods. However, our approach can easily be extended to
include linear combinations of CSFs. Since it does not rely on the absence of
active-to-active exctiations in the cluster operator, it can in princple
therefore be used to correlate RASSCF51 or GASSCF52 references to recover the
dynamic correlation of highly multi-reference states, without spin
contamination.

## Notes and references

  * Bartlett and Musial 2007 R. J. Bartlett and M. Musial, _Rev. Mod. Phys._ , 2007, 79, 291\.
  * Köhn _et al._ 2012 A. Köhn, M. Hanauer, L. A. Mück, T.-C. Jagau and J. Gauss, _WIREs Comput Mol Sci_ , 2012, 3, 176–197.
  * Lyakh _et al._ 2011 D. I. Lyakh, M. Musial, V. F. Lotrich and R. J. Bartlett, _Chem. Rev._ , 2011, 112, 182–243.
  * Evangelista 2018 F. A. Evangelista, _J. Chem. Phys._ , 2018, 149, 030901.
  * Jeziorski and Monkhorst 1981 B. Jeziorski and H. J. Monkhorst, _Phys. Rev. A_ , 1981, 24, 1668\.
  * Banerjee and Simons 1981 A. Banerjee and J. Simons, _Int. J. Quantum Chem._ , 1981, 19, 207–216.
  * Evangelista _et al._ 2012 F. A. Evangelista, M. Hanauer, A. Köhn and J. Gauss, _J. Chem. Phys._ , 2012, 136, 204108.
  * Hanauer and Köhn 2011 M. Hanauer and A. Köhn, _J. Chem. Phys._ , 2011, 134, 204111\.
  * Xiangzhu and Paldus 1993 L. Xiangzhu and J. Paldus, _Int. J. Mol. Sci._ , 1993, 48, 269–285.
  * Xiangzhu and Paldus 1994 L. Xiangzhu and J. Paldus, _J. Chem. Phys._ , 1994, 101, 8812–8826.
  * Xiangzhu and Paldus 1997 L. Xiangzhu and J. Paldus, _Recent Advances in Coupled-Cluster Methods_ , World Scientific, 1997, vol. 3, pp. 183–219.
  * Xiangzhu and Paldus 1998 L. Xiangzhu and J. Paldus, _Mol. Phys._ , 1998, 94, 41–54.
  * Xiangzhu and Paldus 1998 L. Xiangzhu and J. Paldus, _Int. J. Quantum Chem._ , 1998, 70, 65–75.
  * Lindgren 1978 I. Lindgren, _Int. J. Quantum Chem._ , 1978, 14, 33–58.
  * Haque and Mukherjee 1984 M. A. Haque and D. Mukherjee, _J. Chem. Phys._ , 1984, 80, 5058–5069.
  * Stolarczyk and Monkhorst 1985 L. Z. Stolarczyk and H. J. Monkhorst, _Phys. Rev. A_ , 1985, 32, 725\.
  * Stolarczyk and Monkhorst 1985 L. Z. Stolarczyk and H. J. Monkhorst, _Phys. Rev. A_ , 1985, 32, 743\.
  * Stolarczyk and Monkhorst 1988 L. Z. Stolarczyk and H. J. Monkhorst, _Phys. Rev. A_ , 1988, 37, 1908\.
  * Stolarczyk and Monkhorst 1988 L. Z. Stolarczyk and H. J. Monkhorst, _Phys. Rev. A_ , 1988, 37, 1926\.
  * Kaldor 1991 U. Kaldor, _Theor. Chim. Acta_ , 1991, 80, 427–439.
  * Nooijen and Bartlett 1996 M. Nooijen and R. J. Bartlett, _J. Chem. Phys._ , 1996, 104, 2652–2688.
  * Nooijen 1996 M. Nooijen, _J. Chem. Phys._ , 1996, 104, 2638–2651.
  * Nooijen and Bartlett 1997 M. Nooijen and R. J. Bartlett, _J. Chem. Phys._ , 1997, 106, 6441–6448.
  * Nooijen and Bartlett 1997 M. Nooijen and R. J. Bartlett, _J. Chem. Phys._ , 1997, 106, 6449–6455.
  * Nooijen and Bartlett 1997 M. Nooijen and R. J. Bartlett, _J. Chem. Phys._ , 1997, 107, 6812–6830.
  * Nooijen _et al._ 2014 M. Nooijen, D. Ondrej, D. Datta, L. Kong, K. R. Shamasundar, V. F. Lotrich, L. M. Huntington and F. Neese, _J. Chem. Phys._ , 2014, 140, 081102\.
  * Jana _et al._ 2002 D. Jana, U. S. Mahapatra and D. Mukherjee, _Chem. Phys. Lett._ , 2002, 353, 100–110.
  * Datta and Mukherjee 2008 D. Datta and D. Mukherjee, _Int. J. Quantum Chem._ , 2008, 108, 2211–2222.
  * Datta and Mukherjee 2009 D. Datta and D. Mukherjee, _J. Chem. Phys._ , 2009, 131, 044124\.
  * Hermann and Hanrath 2020 N. Hermann and M. Hanrath, _J. Chem. Phys._ , 2020, 153, 164114\.
  * Hermann and Hanrath 2021 N. Hermann and M. Hanrath, _Mol. Phys._ , 2021, 120, e2005836.
  * Hermann and Hanrath 2022 N. Hermann and M. Hanrath, _J. Chem. Phys._ , 2022, 156, 054111\.
  * Helgaker _et al._ 2013 T. Helgaker, J. Olsen and P. Jorgensen, _Modern Electronic Structure Theory_ , Wiley-Blackwell, 2013.
  * Isaiah and Bartlett 2009 S. Isaiah and R. J. Bartlett, _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_ , Cambridge University Press, 2009.
  * Head-Gordon _et al._ 1998 M. Head-Gordon, P. E. Maslen and C. A. White, _J. Chem. Phys._ , 1998, 108, 616–625.
  * Kutzelnigg and Mukherjee 1997 W. Kutzelnigg and D. Mukherjee, _J. Chem. Phys._ , 1997, 107, 432–449.
  * Maitra _et al._ 2012 R. Maitra, D. Sinha and D. Mukherjee, _J. Chem. Phys._ , 2012, 137, 024105.
  * Das _et al._ 2012 S. Das, S. Pathak, D. Datta and D. Mukherjee, _J. Chem. Phys._ , 2012, 136, 164104.
  * Mukherjee _et al._ 1975 D. Mukherjee, R. K. Moitra and A. Mukhopadhyay, _Mol. Phys._ , 1975, 30, 1861–1888.
  * Pal _et al._ 1984 S. Pal, M. D. Prasad and D. Mukherjee, _Theor. Chim. Acta_ , 1984, 66, 311–332.
  * Lindgren and Mukherjee 1987 I. Lindgren and D. Mukherjee, _Phys. Rep._ , 1987, 151, 93–127.
  * Pulay 1980 P. Pulay, _Chem. Phys. Lett._ , 1980, 73, 393–398.
  * Scuseria _et al._ 1986 G. E. Scuseria, T. J. Lee and H. F. Schaefer III, _Chem. Phys. Lett._ , 1986, 130, 236–239.
  * 44 P. G. Szalay and J. Gauss, 107, 9028–9038.
  * Tew _et al._ 2007 D. P. Tew, W. Klopper and T. Helgaker, _J. Comp. Chem._ , 2007, 28, 1307–1320.
  * Kong 2009 L. Kong, _Int. J. Quantum Chem._ , 2009, 110, 2603–2613.
  * Marti-Dafcik _et al._ 2024 D. Marti-Dafcik, N. Lee, H. G. A. Burton and D. P. Tew, _Spin-coupled molecular orbitals: chemical intuition meets quantum chemistry_ , 2024.
  * Kong 2008 L. Kong, _Int. J. Quantum Chem._ , 2008, 109, 441–447.
  * Hanrath 2005 M. Hanrath, _J. Chem. Phys._ , 2005, 123, 084102.
  * Hanrath 2006 M. Hanrath, _Chem. Phys. Lett._ , 2006, 420, 426–431.
  * Malmqvist _et al._ 1990 P. A. Malmqvist, A. Rendell and B. O. Roos, _J. Phys. Chem._ , 1990, 94, 5477–5482.
  * Ma _et al._ 2011 D. Ma, G. Li Manni and L. Gagliardi, _J. Chem. Phys._ , 2011, 135, 044128\.

##  Appendix A Hanauer’s lemma

We gratefully acknowledge the communication of this result and the proof,
which were made known to us by Matthias Hanauer and Andreas Köhn through
private communications.

Hanauer’s lemma states that the product of Hamiltonian and normal ordered wave
operator $\hat{H}\\{e^{\hat{T}}\\}$ can be expressed as the product of
$\\{e^{\hat{T}}\\}$ with a sum of fully connected terms $\tilde{H}_{p}$

$\hat{H}\\{e^{\hat{T}}\\}=\\{e^{\hat{T}}\\}\sum_{p=0}^{\infty}(-1)^{p}\tilde{H}_{p}$
(25)

The proof uses a derivation found in Lindgren’s paper14 that extracts the
connected terms from the product of the Hamiltonian and the normal ordered
exponential form of the wave operator:

$\displaystyle\hat{H}\\{e^{\hat{T}}\\}$

The contracted terms can be separated into products of the connected parts and
the remaining cluster operators:

$\displaystyle\left\\{\wick{\c{\hat{H}}\c{e}^{\hat{T}}}\right\\}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n!}\\{\wick{\c{\hat{H}}\,\\{\c{\hat{T}}^{n}\\}}\\}$
(26) $\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n!}\sum_{k=1}^{n}{n\choose
k}\\{(\hat{H}\,\hat{T}^{k})_{C}\,\\}\\{\hat{T}^{n-k}\\}\\}\\}$ (27)
$\displaystyle=\sum_{(n-k)=0}^{\infty}\sum_{k=1}^{\infty}\frac{1}{k!}\frac{1}{(n-k)!}\\{(\hat{H}\,\\{\hat{T}^{k}\\})_{C}\,\\{\hat{T}^{n-k}\\}\\}$
(28)
$\displaystyle=\sum_{k=1}^{\infty}\frac{1}{k!}\\{(\hat{H}\,\\{\hat{T}^{k}\\})_{C}\,\\{e^{\hat{T}\\}}\\}$
(29)

This leads to the classic result that39, 14

$\displaystyle\hat{H}\\{e^{\hat{T}}\\}=\sum_{k=0}^{\infty}\frac{1}{k!}\\{(\hat{H}\,\\{\hat{T}^{k}\\})_{C}\,\\{e^{\hat{T}}\\}\\}=\\{\\{e^{\hat{T}}\\}(\hat{H}\\{e^{\hat{T}}\\})_{C}\\}$
(30)

where the subscript $C$ denotes connected terms. An entirely analogous
argument gives the result

$\\{e^{\hat{T}}\\}\hat{H}=\\{\\{e^{\hat{T}}\\}(\\{e^{\hat{T}}\\}\hat{H})_{C}\\}$
(31)

noting that, within a normal ordered product, operators consisting of even
numbers of creation and annihilation operators commute. The same argument can
be followed when $\hat{H}$ is replaced by the connected composite
$\tilde{H}_{0}=(\hat{H}\\{e^{\hat{T}}\\})_{C}$ to give

$\begin{split}\\{e^{\hat{T}}\\}\tilde{H}_{0}&=\\{\\{e^{\hat{T}}\\}(\\{e^{\hat{T}}\\}\tilde{H}_{0})_{C}\\}\\\
&=\\{\\{e^{\hat{T}}\\}\tilde{H}_{0}\\}+\\{\\{e^{\hat{T}}\\}(\\{e^{\hat{T}}-1\\}\tilde{H}_{0})_{C}\\}\\\
&=\\{\\{e^{\hat{T}}\\}\tilde{H}_{0}\\}+\\{\\{e^{\hat{T}}\\}\tilde{H}_{1}\\}\end{split}$
(32)

Here we have defined the fully connected term
$\tilde{H}_{1}=(\\{e^{\hat{T}}-1\\}\tilde{H}_{0})_{C}$. Similarly, by defining
the fully connected composites
$\tilde{H}_{n+1}=(\\{e^{\hat{T}}-1\\}\tilde{H}_{n})_{C}$ and rearranging, we
find a recursion relation

$\begin{split}\\{\\{e^{\hat{T}}\\}\tilde{H}_{n}\\}&=\\{e^{\hat{T}}\\}\tilde{H}_{n}-\\{\\{e^{\hat{T}}\\}\tilde{H}_{n+1}\\}\\\
&=\\{e^{\hat{T}}\\}\tilde{H}_{n}-\\{e^{\hat{T}}\\}\tilde{H}_{n+1}+\\{\\{e^{\hat{T}}\\}\tilde{H}_{n+2}\\}\\\
&=\\{e^{\hat{T}}\\}\sum_{p=n}^{N}(-1)^{p-n}\tilde{H}_{p}+(-1)^{N-n+1}\\{\\{e^{\hat{T}}\\}\tilde{H}_{N+1}\\}\end{split}$
(33)

For sufficiently large $N$, $\tilde{H}_{N+1}$ will be of sufficiently high
excitation rank as to completely annihilate the reference function when
applied in a normal ordered product:
$\\{\\{e^{\hat{T}}\\}\tilde{H}_{N+1}\\}\ket{\Phi_{0}}=0$. Therefore, in a
finite basis we can take $N$ to $\infty$ and remove the last term in equation
33, giving

$\begin{split}\\{\\{e^{\hat{T}}\\}\tilde{H}_{n}\\}\ket{\Phi_{0}}&=\\{e^{\hat{T}}\\}\sum_{p=n}^{\infty}(-1)^{p-n}\tilde{H}_{p}\end{split}$
(34)

Substituting the case for $n=0$ back into Eq. 30 yields the desired result

$\hat{H}\\{e^{\hat{T}}\\}=\\{\\{e^{\hat{T}}\\}\tilde{H}_{0}\\}=\\{e^{\hat{T}}\\}\sum_{p=0}^{\infty}(-1)^{p}\tilde{H}_{p}$
(35)

##  Appendix B The Connectedness of the residual equation

Using Hanauer’s lemma, the unlinked energy expression can be written as

$\begin{split}\braket{\Phi_{0}}{\hat{H}\\{e^{\hat{T}}\\}}{\Phi_{0}}&=\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}\sum_{p=0}^{\infty}(-1)^{p}\tilde{H}_{p}}{\Phi_{0}}=\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}\hat{\Lambda}}{\Phi_{0}}\\\
&=\sum_{K}\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\\\
&=\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{\Lambda}}{\Phi_{0}}\\\
&\quad+\sum_{K\neq
0}\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\\\
&=\braket{\Phi_{0}}{\hat{\Lambda}}{\Phi_{0}}+\sum_{K\neq
0}\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\end{split}$
(36)

In the above, we have replaced the sum of connected terms
$\sum_{p=0}^{\infty}(-1)^{p}\tilde{H}_{p}$ with $\hat{\Lambda}$ for the sake
of brevity.

Similarly, the unlinked residual equation can be written as

$\begin{split}\braket{\Phi_{I}}{\hat{H}\\{e^{\hat{T}}\\}}{\Phi_{0}}&=\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{\Phi_{0}}\braket{\Phi_{0}}{\hat{\Lambda}}{\Phi_{0}}\\\
&\quad\quad+\sum_{K\neq
0}\braket{\Phi_{i}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\\\
&=\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{\Phi_{0}}\Big{(}\braket{\Phi_{0}}{\hat{H}\\{e^{\hat{T}}\\}}{\Phi_{0}}\\\
&\hskip 71.13188pt-\sum_{K\neq
0}\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\Big{)}\\\
&\quad\quad+\sum_{K\neq
0}\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\\\
\end{split}$ (37)

where we substituted the result from Eq 36 in the second line. Rearranging and
simplifying gives the result that our residual equation becomes

$\begin{split}R_{I}&=\braket{\Phi_{I}}{\hat{H}\\{e^{\hat{T}}\\}}{\Phi_{0}}-\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{{\Phi_{0}}}\braket{\Phi_{0}}{\hat{H}\\{e^{\hat{T}}\\}}{\Phi_{0}}\\\
&=\sum_{K\neq
0}\Big{(}\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}-\braket{\Phi_{I}}{\\{e^{\hat{T}}\\}}{\Phi_{0}}\braket{\Phi_{0}}{\\{e^{\hat{T}}\\}}{\Phi_{K}}\Big{)}\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}\end{split}$
(38)

Since this is a sum of the connected quantities
$\braket{\Phi_{K}}{\hat{\Lambda}}{\Phi_{0}}$, we conclude that solving for
$R_{I}=0$ is equivalent to solving a totally connected form of the equations.
Our equations are therefore size-extensive order-by-order in $\hat{T}$.