# Topic Modeling in Density Functional Theory on Citations of Condensed Matter
Electronic Structure Packages

Marie Dumaz Corresponding Author<EMAIL_ADDRESS>Department of Physics
and Astronomy, West Virginia University, Morgantown, WV 26506-6315, USA  Lane
Department of Computer Science and Electrical Engineering, West Virginia
University, Morgantown, WV 26506-6315, USA  Camila Romero-Bohórquez
Neuroscience Department, West Virginia University, Morgantown, WV 26506-6315,
USA  Donald Adjeroh  Lane Department of Computer Science and Electrical
Engineering, West Virginia University, Morgantown, WV 26506-6315, USA  Aldo H.
Romero  Department of Physics and Astronomy, West Virginia University,
Morgantown, WV 26506-6315, USA

###### Abstract

With an increasing number of new scientific papers being released, it becomes
harder for researchers to be aware of recent articles in their field of study.
Accurately classifying papers is a first step in the direction of personalized
catering and easy access to research of interest. The field of Density
Functional Theory (DFT) in particular is a good example of a methodology used
in very different studies, and interconnected disciplines, which has a very
strong community publishing many research articles. We devise a new
unsupervised method for classifying publications, based on topic modeling, and
use a DFT-related selection of documents as a use case. We first create topics
from word analysis and clustering of the abstracts from the publications, then
attribute each publication/paper to a topic based on word similarity. We then
make interesting observations by analyzing connections between the topics and
publishers, journals, country or year of publication. The proposed approach is
general, and can be applied to analyze publication and citation trends in
other areas of study, beyond the field of Density Function Theory.

##  1 Introduction

For many years, scientific documents in basic sciences were classified using
the Physics and Astronomy Classification Scheme (PACS) method, which maps a
series of codes with keywords to cover the main topics of the publication.
This is usually done by the author, and provides relative freedom for the
actual classification. Though this is still one of the methods used by several
journals, this step is transitioning into a more focused type of topic
classification[13]. Basically, this change is driven by the fact that the
institute that introduced this approach (the American Institute of Physics,
AIP) has decided to discontinue its support, primarily due to the huge cost of
administration and paper handling during manuscript submission. Therefore, a
more dynamical and flexible method is necessary. A method were different
topics capture the continuous change in the scientific development and
scientific interest and are able to map reviewers and even journal in this
arena. In this respect, the use of topic modeling from machine learning is
gaining a lot of attention [14, 2]. In this paper, we propose an unsupervised
learning technique that focuses on finding text correlations and thus classify
documents according to a small set of topics. These topics are found through
an unsupervised machine learning technique that creates clusters of meaningful
words in the document selection, by finding similarity between them in a
latent space. Though the idea of topic modeling has been used in the past to
classify research topics [10, 15], in this work, we start from a well-defined
database, and apply this methodology on publications related to the use of
electronic structure packages.

While a general method to categorize topics in basic sciences would be ideal,
in this paper, we have decided to focus in a particular research field, where
we can have access to a very well controlled database [6]. Density Functional
Theory (DFT) is the most used methodology to characterize materials at the
atomic scale [9, 8]. It falls into the so called electronic structure methods
and it has the benefit that a large number of computational codes and software
employ this method, which have created communities around them, supporting the
research efforts. We base our document selection on research papers citing
those computational code libraries, in order to reduce our scope to only DFT-
related papers.

Some efforts to analyze DFT-related papers from different subfields have been
made, but they classify the articles based on the PACS numbers [1, 5].
However, our dataset downloaded from the Web of Science database, does not
keep records of the PACS numbers and only retains keywords. Moreover, as DFT
is a methodology with many different applications, in fields other than
physics, not all papers follow the PACS structure. This shows again, the
necessity for a new classification system, based purely on the content of the
documents. Our approach is not only able to classify documents on content but
also can be used to find novel topics and relationships between topics which
can be used to identify new and emerging research fields.

##  2 Methodology

The research papers and their abstracts were obtained using the same
methodology as in previous work [6], through the Web of Science database
(WoS), but expanded to include 2020 articles.

Once we extracted the abstracts from the WoS database, we performed several
pre-processing steps to capture only the most meaningful words out of the
text. First, we deleted certain groups of words, such as “left arrow” or
“vertical bar,” that are only used to create mathematical representations. For
the same reason, we removed all characters between two dollar signs, as this
is a standard command to make equations, which would only confuse the
algorithm. There is one exception to this rule if an underscore is found
between the two dollar signs, as that pattern is used for chemical formulas
that are essential and meaningful notions in sub-fields of DFT. We also remove
publishing and copyright tags, DOI numbers, and all non-alphanumeric
characters.

Then, the leftover sentences get tokenized, which means the text is broken
down into units, creating a list of tokens from the text. From those lists, we
take out all stopwords, or words that are so common they do not help
differentiate different sub-fields. We used the standard NLTK [3] list of
stopwords and added others based on our own experience with the data after
rounds of trial and feedback. For example, “theory” and “method” are not
specific to some areas of physics like magnetic properties or gas absorption.
We also remove all single letters and digits, often used in mathematical
expressions or to quantify a value and don’t hold meaning intrinsically. We
then create bigrams and trigrams, groups of two or three tokens that often
appear one after the other and will now be linked as one word by an
underscore. Examples would be “ab initio” and “density functional theory.”
Finally, the tokens are lemmatized to reduce words to their root, and the
resulting list of tokens is filtered to only keep the ones that are in at
least two documents and in less than 20% of the papers. From this final list
are created a dictionary and a Bag of Words (BoW) which are the input to our
algorithm.

Latent Dirichlet Allocation (LDA) [4] is a topic modeling algorithm which uses
Bayesian statistical methods to generate text, from documents passed in as
input in the BoW. This means that LDA does not take into account the order of
words. In fact, an LDA model creates a set of topics where each topic is a
probability distribution over a corpus of words drawn from the distribution
$\beta_{k}\sim\mbox{Dirichlet}(\eta)$, and each document is a probability
distribution over topics drawn from $\theta_{k}\sim\mbox{Dirichlet}(\alpha)$.

Here, we define $M$ as the total number of documents, $N$ as the total number
of words in a document $m$, and $k$ as the number of topics. As shown in
Figure 1, LDA spans over three layers. The parameters $\alpha$ and $\eta$ are
corpus-level parameters, which means that they are sampled only once per
corpus. $\theta$ denote document-level variables, sampled once per document.
Finally, variables $z_{n}$ and $w_{n}$ are at word-level: they are sampled
once for each word in each document. $w_{n}$ is also the only observable
variable, picked from the bag of words given as input.

Figure 1: Graphical representation of a LDA model [4]

For topic modeling, we implemented a Latent Dirichlet Allocation (LDA)
algorithm with the Python library Gensim [11] that already has a handy set of
functions to create, evaluate and visualize LDA models. We used two bag of
words (BOWs): one with all documents in our dataset from 1990 to 2019, and one
with only documents published in 2020. The 2020 corpus was used as a test set.

We created 16 models with a different number of topics ($k$), from 5 to 80,
while keeping all other parameters null. For this analysis, we chose to tune
the parameters left for 15, 25, 35 and 45 topics. All models were evaluated
with the $C_{v}$ coherence score [12], which is the coherence score with the
strongest correlation to human rankings. We then ran 480 models to tune
hyperparameters including $\alpha$ and $\eta$, for each of the 5 different
values of $k$. For each value of $k$, the model that gives the highest $C_{v}$
score is used through the rest of this article.

To speed up training, we run the parameter tuning on several cores. We should
note that the stochastic optimization process used in online LDA includes the
randomization of data partition, which makes it impossible to ensure exact
reproducibility if coupled with parallelization.

To discover trends and interests, we assign topics to our publications
bibliometric attributes. We assign topics by considering the probabilities
returned by the LDA model as weight. We create a matrix $A$ where each row
represents a document, and each column represents a topic. The elements of
this matrix are the probabilities for a given topic to be assigned to a given
document.

This matrix is particularly useful for the analyses made of topic
distributions over countries, journals and years. This way, when computing the
number of publications per topic per country, for example, we can use those
weights to perform a more precise analysis. For these analyses, we create a
matrix where each row is a topic and each column is a country. If a document
was created in the USA for example, then each topic would add their respective
weight for this document in the “USA” column. These analyses are also all
normalized, which means each value in the matrix is divided by the sum of all
values in the $A$ matrix, which is equal to the total number of documents.

##  3 Results

Here, we describe our dataset, experiments and results.

To create a dataset of research papers focused on DFT, we used the graphical
interface “Web of Science” (WOS) and downloaded all articles citing some of
the most commonly used computational packages to study crystalline systems. We
acquired all papers from the first citation up to 2019 for our training set;
and all citations in the 2020 year for our test set. [6]

###  3.1 Distribution of documents and topics

Figure 2: Document distribution over topics (for publications from 1990-2019)

Figure 2 shows the distribution of documents in our dataset over topics. This
means that it displays how many documents contain each topic. Note, we use the
probabilities of a topic appearing in a document as weight, as described in
the Methodology section. We add up the weights for all documents in each
topic, and divide the results by the total number of documents. The figure
gives us an idea of how well distributed the topics identified from each model
is, and what topic is more common within the documents.

We can see that the model with 25 topics is extremely unbalanced, with only 4
topics accounting for up to 80% of the documents. Next, the LDA models for $k$
= 15 and $k$ = 45 seem to be somewhat well distributed, as no topics is
contained in more than 18% and 11% of the documents. This can easily be
deduced via the y-axis of the histograms, that are also a quick way to get a
general idea of a distribution: if the y-axis range is larger, it means that
at least one topic takes a bigger part of the documents, which means the
distribution is less uniform. Similarly, we can study the range of the
attributed normalized number of documents. The model with 25 topics has the
bigger range at 0.255, followed by 15 topics at 0.175 and 45 topics at 0.105.
The smallest range, and hence the model that has a more uniform distribution
is for $k$ = 35 with 0.086. Indeed, to get 90% of the documents, it takes the
first 8 topics (or 53% of the topics) with the highest proportion for $k$ =
15, 24% of the topics for $k$ = 25, 57% for $k$ = 35 and 37% for $k$ = 45.
Moreover, the model with 35 topics is also the one that has the highest
minimum percentage of documents.

The diversity in the models for $k$ = 35 and $k$ = 45 also seems to imply that
15 topics is not enough for our corpus. In fact, we find evidence for topics
specialization that further confirm the benefits of a higher number of topics.
For example, in the $k$ = 15 model, Topic 3 contains the following words:
“formation”, “activity”, “hydrogen”, “metal”, “reaction”, “mechanism”,
“process”, “catalyst”, “oxidation”, and “li”. We can map this topic to two
different topics in the $k$ = 35 model: it shares the words “activity”,
“metal”, “reaction”, “catalyst”, “oxidation” with Topic 6 and the words
“reaction”, “formation”, “mechanism”, “process” with Topic 27. We can
distinctly see that it splits into two categories, one is more focused on
specific processes and applications (Topic 6) and one is more general (Topic
27).

We can conclude from these analyses that the most balanced model is for $k$ =
35. For these two reasons, we will mostly focus on this particular LDA model
and its results in the rest of this section.

Figure 3: Top 10 words making up each topic in the LDA model for $k$ = 35

###  3.2 Closer look at the Topics

Figure 3 shows the word clouds for the 35 topics in our LDA model. Each word
cloud represents the 10 words with the highest probability in the topic. The
size of each word is proportional to their probability within their topics.
Bigger words are terms with a higher probability and a bigger impact in the
topic. We also added topic names to give a quick idea of what each word cloud
represents. Even though some topics overlap, most topics are rather
independent and clearly defined. For example, in Topic 13, we can find words
relating to magnetic problems, while Topic 7 was more about hydrogen and
hydride. Hydrogen is not a magnetic atom so those two subjects are not
expected to be clustered together, and indeed, they were not: they are quite
different. Topic 1 on alloy, however, can be related to magnetism as some
alloy metals are magnetic and yet, it is once again put into a separate group
(Topic 13), which shows how good the LDA model is at detecting different word
relationships. Topic 14 is another good example of terms that could apply to
many different fields with words such as “data”, “simulation”, “model” or
“approach”, but are associated to each other in a unique topic.

Figure 4: T-SNE reduction for the LDA model for $k$ = 35

We can also analyze the topics of our model by visualizing them with the t-SNE
reduction method [7]. Figure 4 displays the data reduction to 2-dimensions,
where each point is a document, positioned depending on their similarity with
other documents, and each color is their dominant topic. Note that some topics
are not represented as they are not a dominant topic in any document, this is
the case for topics 0 and 29 for example. Since the topics are well
distributed over the documents, we can see a lot of topics being represented,
instead of a few overshadowing the entire map. We can also notice that some
topics are more dispersed, such as the subject of new methods (e.g., Topic
31). That is expected as researchers are developing new processes and data
simulations in many different fields. However, other topics are clustered, and
clearly defined like the first topic on alloying and crystal formation. This
figure also helps visualizing the similarity between topics. The closer two
topics are, the more interconnected they are. For example, Topic 27 (energy
barriers/reaction paths/diffusion/activation processes) is close to Topic 6
(catalysis/oxidation) and Topic 3 (polymers/molecular properties) but is on
the opposite side of Topics 1 (alloying/crystal formation) and 13 (magnetic
properties). Indeed, molecular properties and activation processes like
catalysis are chemistry related while magnetism and alloying are more often
studied in physics.

Figure 5: Document distribution over topics, for $k$=35 topics. Left:
1990-2019 corpus. Right: 2020 corpus

##  4 Temporal evolution of topics

Figure 5 compares the distribution of documents over topics, for the 1990-2019
corpus on the left (that was already shown in Figure 2), and for the 2020
corpus on the right. The two plots are similar and both display a good
distribution of the documents. While Topic 10 (temperature/pressure crystal
phase transitions) and Topic 14 (atomistic modelling/classical
potential/classical molecular dynamics) were the most common subjects up to
2019, a number of other topics caught up to them in 2020, including Topic 31
(development of new methods) which now has the highest number of assigned
documents. This could also mean that while Topic 8 (materials
growth/interfaces/electron transport) and Toipc 10 were very important up to
2019, and still are, Topic 4 (monolayers/topological properties), Topic 24
(electrochemical processes/lithium batteries) and Topic 31 (development of new
methods) are gaining more momentum recently. The similarity between the two
plots also indicate that our model is generalizing correctly and can
efficiently attribute topics to new documents.

Figure 6: Normalized and weighted number of documents over the years, per
topic for $k$ = 35. Number inside each plot (top left) denotes the best
fitting power factor.

We computed, for each topic independently, the normalized and weighted number
of documents per year. Every year, there is an exponential amount of DFT-
related articles being published, and consequently, almost all topics also
grew exponentially. To compare them, we calculated the power law factor of
each curve since 2000. Even though some topics’ evolution are better fitted
with an exponential curve, the power law function yields the highest mean
$r^{2}$ over all topics. The results are shown in Figure 6, where the number
in each top left corner of the subplots are the power factor best fitting for
each topic. We chose to only consider publications after 2000 as this is the
period in which DFT software really started to to get globally used.

A total of nine topics have a power factor lower than 1 and while they still
grew in importance since 2000, a few of them are starting to stabilize and
even loose some interest from the DFT community. For example, Topic 18 on DFT
schemes and Topic 19 on nanotubes, which have the lowest factors recorded,
have not gained any new momentum in the recent years and Topic 18 also
observed a decrease of publications in 2017 for the first time ever. On the
contrary, Topics 24 related to electrochemical processes and lithium
batteries, as well as Topic 4 (monolayers/topological properties), Topic 6
(catalysis/oxidation) and Topic 31 (development of new methods) are growing at
a significantly high rate and seem to be among the most important subjects in
the field of DFT today. Other subjects are all developing at intermediate
rates, representing the diversity of DFT and how important the methodology is
for different topics.

Figure 7: Comparison of the two topic with the lowest and highest growth, fit
to the best exponential curve and powerlaw curve. Left: Topic 18 with lowest
power law factor. Right: Topic 24 with highest power law factor.

We report in Figure 7 the two topics with the smallest and highest growth out
of the 35 in our model. We fit a power law and an exponential curve to both of
the curves and reported the factors in the top left corner. Topic 24
(electrochemical processes/lithium batteries) is an example of a growth better
fitted with an exponential function, getting a $r^{2}$ of 0.998, even though a
power law function also get a high $r^{2}$ score of 0.988. On top of the
important difference in growth, we can also notice how “bumpy” the evolution
of Topic 18 (DFT schemes) is. The many irregularities, compared to the smooth
curve of Topic 24 might also be a sign of waning (or at least, inconsistent)
interest.

Figure 8: Normalized and weighted number of documents over the years, per
topic, for $k$ = 35

Figure 8 also shows the growth of topics over the years, in the form of a heat
map. We notice that no important increase happened before 2000, which confirms
our choice to study the exponential factors only after 2000. Moreover, most
topics only start gaining momentum after at least 2005. It is interesting to
note that the three topics with the highest exponential factors (Topic 4, 6
and 24) grew extremely fast in a few years, which is the reason for the high
factor. Topic 24 especially only started to differentiate itself from the
others just in 2018. It is also interesting how out of those three topics, two
of them (Topics 6 and 24) are not part of the topics with the highest
proportion of documents. However, given their high power factor, we can only
predict that they will become more important in the next few years. Similarly,
by looking at the trajectory of some subjects such as Topics 1 or 17, we can
confidently forecast that the interest in those subfields will grow in the
coming years. The 2020 column in this figure seems to agree with those
predictions as Topics 1, 4, 6, 17 and 24 all gained in proportion during the
last year. On the other hand, we can also easily determine which topics are
stabilizing and are loosing momentum. For example, Topic 3 (polymers/molecular
properties) and Topic 18 (DFT schemes/approximations in DFT) have had interest
from the scientific community for at least a decade and yet, are not growing
as fast as others.

It is interesting to compare the evolution of Topics 14 and Topic 31, as they
are related to each other. Topic 14 contains words such as “data”, “model” and
“simulation” which represent the community’s effort to incorporate code,
machine learning and new computing methods into their work. This topic has
been growing steadily over the years, as it becomes a bigger part of today’s
experiments. Secondly, Topic 31 relating to the improvement and optimization
of code and processes, only started gaining momentum in 2014, about 10 years
after Topic 14. However, its growth was faster, which places both of the
topics at the top of the most published subjects in 2020 in DFT. We can also
mention the evolution of Topic 4 related to graphene and monolayers. It
started to grow a few years after 2010, which is the year that the Nobel Prize
in physics was awarded to Andre Geim and Kostya Novoselov for their work
regarding the two-dimensional material graphene.

Figure 9: Normalized and weighted number of citations over the years, per
topic, for $k$ = 35

##  5 Trends in citations

Looking at the number of citations per year, per topic in Figure 9, we notice
the same trends as noted in the previous figure, only pushed a few years back.
Indeed, papers published today are citing older research articles so a raise
in publications in 2019 means a raise in citations of older papers, typically
at least 3 years prior. This gap is surely due to the research process of
reading the literature, putting together a new idea and hypothesis, getting
data and experiences and finally drawing conclusions and then publishing.

Hence, looking at the trends of citations over the years might give us more
insight and more predictive value than only the number of published papers.
For topics that are only starting to attract more interest, like Topics 6 and
24, we can trace back the beginning of the trend to the late 2000s or early
2010s whereas topics that have been in place for longer such as Topics 8, 10
or 14 can be attributed to discoveries in the late 1990s.

Figure 10: Normalized and weighted number of documents per country (only top
15 shown), per topic, in different time periods, for $k$ = 35

###  5.1 Spatial evolution of topics

We also studied the spatial evolution of topics in DFT over the years,
especially, by countries. To do so, we identify countries using the
affiliations of all authors. However, we only count a country once for each
paper.

Figure 10 divides the heat maps of topics per country over four different time
periods. Before 2005, the USA were dominating the field of DFT and published
the highest number of publications in a larger number of topics. Then, the
proportions slowly shifted to China, which now generates more articles, and
diversified its interests. China started to get involved in DFT after 2005
through the subject of crystal phase transitions (Topic 10). Pressure crystal
phase transitions is the main strength of DFT and one of the most active
subfields, which explains why it would be the first topic to attract interest.
Topic 10 stayed an important part of work in China until now, as it was still
one of the main published subjects by authors from China between 2015 and
2019, along with Topics 4 and 12. Nowadays, the USA still contributes to many
different topics in DFT, and grew interest in Topic 14, but decreased its
number of publications in Topic 5 (adsorption/diffusion).

We can accurately predict that China will soon dominate the field of DFT as
its number of publications will only keep growing. Figure 11, representing the
normalized distribution of documents per topics over countries, for the year
2020 confirms the influence of China in DFT. Moreover, even though most
software counted in this work are developed in Europe, they are the most used
in China and the USA. Germany, the third country with the highest number of
publications, and country of origin to Turbomole, the sixth most used code in
our dataset is barely even competing with China and the USA after 2005.
Similarly, the impact of the United Kingdom and France kept decreasing over
the years, even though they were very active in several different topics
before 2005.

Figure 11: Normalized and weighted number of documents per country (only top
15 shown), per topic, for articles published in 2020, for $k$ = 35

##  6 Topics distribution by journals and publishers

Figure 12: Normalized and weighted number of documents per journal (only top
15 shown), per topic, in different time periods, for $k$ = 35

Similarly to figure 10, figure 12, shows the normalized number of publications
(computed with the probabilities as weight) per topic over journals, dividing
the papers over four 5-year periods based on their publication year. In the
same manner, these plots express the diversification of DFT publications over
different journals, while still highlighting the importance of a very small
number of prominent members. However, while figure 10 displays signs of
diversification amongst countries as soon as 2005, most papers were still
published in Physical Review B until a lot more recently, and no other journal
seems to be taking over its number of publications.

The 15 journals displayed have the highest number of publications and are
arranged in order from left to right. Even though some journals like Surface
Science are specializing in one topic, and sticking with it over time, we can
notice how most journals are present in several different subjects. Two key
journals, Physical Review B and Journal of Physical Chemistry C are both
influential in DFT and were expected to be the largest contributors to the
field. Moreover, as DFT is mostly used in physics and chemistry, it is not
surprising that those two journals both focused on condensed matter analysis
each represent one of the two domains.

In fact, the three journals in positions 2, 3 and 4 are related to the the
chemistry side of DFT. Even though, they have similar targets, they evenly
distribute the publications between them. The Journal of Chemical Physics
focuses on Topic 14 and 18, Journal of Physical Chemistry C contributes more
to Topics 5, 6, 8 and a few others while Physical Chemistry Chemical Physics
has a more diverse range of subjects. Hence, it seems that the Journal of
Chemical Physics contributes more to methodologies and not material
characteristics, judging by the words present in Topics 14 and 18
(modeling/approximations). Even though physics researchers mostly publish
their articles in one particular journal, chemists have more options available
depending on their interests.

Figure 13: Normalized and weighted number of citations per topic and per
journal (only top 15 shown), with only open access journals considered, for
$k$ = 35

Figure 13 seem to confirm our deductions from the previous figure. However,
even though journals publish papers in a lot of different topics, they seem to
be cited for only a few, which might steer their editorial agenda in the
future. In fact, the four topics most cited in Physical Review B, are also the
four topics most published in the same journal. With such hypothesis we could
predict a raise in the number of papers published in Topics 16 and 27, as they
are more cited and published, as of 2019. We should note that even tough the
values are normalized, a few highly cited articles might skew the results.

Figure 14: Normalized and weighted number of documents per topic and per
publisher (only top 10 shown), with only open access journals considered, for
$k$ = 35.

Figure 14 is similar to Figure 12, but the x-axis shows publishers (and not
journals). The other difference is that we created this plot with only
articles published in open-access journals. It is interesting to see how the
two figures have different interests. For example, open-access journals seem
to have a strong interest in Topic 31 while subscription journals in Figure 12
only published a small amount of their articles relating to this subject. The
other three most important topics in open-access journals are Topics 4, 8 and
10 which are also well published in paid journals but to different levels.
While Topics 4 and 8 grew in proportion, Topics 10, 13 and 14 all saw their
share of documents decrease in open-access publishers, even though they still
remain relevant. Moreover, it seems that open-access journals are more focused
on a few topics, instead of spreading to many different subjects, like paid
journals. We can also notice the rather good correlation with the evolution of
topics over time as the most published papers in open-access journals are also
the ones that are the most published in 2020 (see Figure 8).

##  7 Discussion

Bibliometrics and topic modeling are powerful tools that can be used to study
and evaluate the evolution and impact of any given field. In this paper, we
developed several analyses related to the evolution of topics with emphasis on
the specific field of DFT. Here, evolution was considered in terms of time
(over the years), and space (across countries), and also with respect to
journals, publishers, institutions, and subjects based on information from
over 120,000 papers in Density Functional Theory. We find that automatically
classifying papers based on their abstract is a promising option for improved,
faster, and more objective categorization. In fact, most of the 35 topics
resulting from our LDA model were consistent with our document corpus and
could be attributed a name.

Even though most packages are developed in Europe, they are mostly used in the
USA and China, that are dominating the amount of publications created in the
DFT field. While the two countries share interests in a lot of different
subjects, American scientists work with many of the most important DFT
packages, and researchers in China use mostly VASP, which is a licensed
package. This difference might come from the intrinsic set of tools each
software has, and the subfields they serve the most. Indeed, China seems to
direct most of its work to materials science and applications as they focus on
2D materials (Topic 4), chemical reactions (Topic 6) or batteries (Topic 24)
while researchers in the USA published most of their work in simulations
(Topic 14) and new methodologies (Topic 31).

Moreover, two of the most recurrent topics in our analysis, and that receive
the most interest from influential journals and countries are related to code
development, simulations and modeling. As our work focuses on computational
packages, it is not surprising that some papers focus on the methods and
processes of software. However, the large amount of interest in the subjects
suggests a recent and strong effort in the physics and chemistry fields to
develop and optimize code packages and libraries.

The topic modeling done in this article was created using the abstracts of
papers and already resulted in coherent and clear topics. Future work could
include replicating our methodology on the full-text of all 120,000 papers in
our dataset. Further, the approach used is quite general, and can easily be
adapted to apply to other fields of study, beyond DFT.

## Additional Information

The authors declare no competing interests.

## Data Availability

The dataset that supports the findings of this study are available in Figshare
with the identifier https://doi.org/10.6084/m9.figshare.12494654.v2.

## References

  * [1] Alberto Aleta, Sandro Meloni, Nicola Perra and Yamir Moreno  “Explore with caution: mapping the evolution of scientific interest in physics”  In _EPJ Data Science_ 8.1 Springer Berlin Heidelberg, 2019, pp. 27
  * [2] Bhagyashree Vyankatrao Barde and Anant Madhavrao Bainwad  “An overview of topic modeling methods and tools”  In _2017 International Conference on Intelligent Computing and Control Systems_ , 2017, pp. 745–750  IEEE
  * [3] Steven Bird, Ewan Klein and Edward Loper  “Natural language processing with Python: analyzing text with the natural language toolkit”  ” O’Reilly Media, Inc.”, 2009
  * [4] David M Blei, Andrew Y Ng and Michael I Jordan  “Latent dirichlet allocation”  In _Journal of machine Learning research_ 3.Jan, 2003, pp. 993–1022
  * [5] Matteo Chinazzi, Bruno Gonçalves, Qian Zhang and Alessandro Vespignani  “Mapping the physics research space: a machine learning approach”  In _EPJ Data Science_ 8.1 Springer Berlin Heidelberg, 2019, pp. 33
  * [6] Marie Dumaz, Reese Boucher, Miguel AL Marques and Aldo H Romero  “Authorship and citation cultural nature in Density Functional Theory from solid state computational packages”  In _Scientometrics_ 126.8 Springer, 2021, pp. 6681–6695
  * [7] Geoffrey E Hinton and Sam Roweis  “Stochastic neighbor embedding”  In _Advances in Neural Information Processing Systems_ 15, 2002
  * [8] Pierre Hohenberg and Walter Kohn  “Inhomogeneous electron gas”  In _Physical review_ 136.3B APS, 1964, pp. B864
  * [9] Walter Kohn and Lu Jeu Sham  “Self-consistent equations including exchange and correlation effects”  In _Physical review_ 140.4A APS, 1965, pp. A1133
  * [10] Michael Paul and Roxana Girju  “Topic modeling of research fields: An interdisciplinary perspective”  In _Proceedings of the International Conference RANLP-2009_ , 2009, pp. 337–342
  * [11] Radim Řehůřek and Petr Sojka  “Software Framework for Topic Modelling with Large Corpora”  In _Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks_ Valletta, Malta: ELRA, 2010, pp. 45–50
  * [12] Michael Röder, Andreas Both and Alexander Hinneburg  “Exploring the space of topic coherence measures”  In _Proceedings of the eighth ACM international conference on Web search and data mining_ , 2015, pp. 399–408
  * [13] Arthur Smith  “From PACS to PhySH”  In _Nature Reviews Physics_ 1.1 Nature Publishing Group, 2019, pp. 8–11
  * [14] Ike Vayansky and Sathish AP Kumar  “A review of topic modeling methods”  In _Information Systems_ 94 Elsevier, 2020, pp. 101582
  * [15] Weizhong Zhao, Wen Zou and James J Chen  “Topic modeling for cluster analysis of large biological and medical datasets”  In _BMC bioinformatics_ 15.11, 2014, pp. 1–11  BioMed Central

## Author Contributions

M.D. and C.R.-B. downloaded the data from the Web of Science database and
preformed data cleaning and data pre-processing. M.D. executed model training
and testing. M.D., D.A. and A.H.R. did the data analysis and interpretation,
created the data figures and wrote the main text. All authors reviewed the
manuscript.

## Figures

Fig. 1 Graphical representation of a LDA model [4]

Fig. 2 Document distribution over topics (for publications from 1990-2019)

Fig. 3 Top 10 words making up each topic in the LDA model for $k$ = 35

Fig. 4 T-SNE reduction for the LDA model for $k$ = 35

Fig. 5 Document distribution over topics, for $k$=35 topics. Left: 1990-2019
corpus. Right: 2020 corpus

Fig. 6 Normalized and weighted number of documents over the years, per topic
for $k$ = 35. Number inside each plot (top left) denotes the best fitting
power factor.

Fig. 7 Comparison of the two topic with the lowest and highest growth, fit to
the best exponential curve and powerlaw curve. Left: Topic 18 with lowest
power law factor. Right: Topic 24 with highest power law factor.

Fig. 8 Normalized and weighted number of documents over the years, per topic,
for $k$ = 35

Fig. 9 Normalized and weighted number of citations over the years, per topic,
for $k$ = 35

Fig. 10 Normalized and weighted number of documents per country (only top 15
shown), per topic, in different time periods, for $k$ = 35

Fig. 11 Normalized and weighted number of documents per country (only top 15
shown), per topic, for articles published in 2020, for $k$ = 35

Fig. 12 Normalized and weighted number of documents per journal (only top 15
shown), per topic, in different time periods, for $k$ = 35

Fig. 13 Normalized and weighted number of citations per topic and per journal
(only top 15 shown), with only open access journals considered, for $k$ = 35

Fig. 14 Normalized and weighted number of documents per topic and per
publisher (only top 10 shown), with only open access journals considered, for
$k$ = 35.


# Discrete Newtonian dynamics with Nosé-Hoover thermostats

Søren Toxvaerd  Department of Science and Environment, Roskilde University,
Postbox 260, DK-4000 Roskilde<EMAIL_ADDRESS>
###### Abstract

Almost all Molecular Dynamics (MD) simulations are discrete dynamics with
Newton’s algorithm first published in 1687, and much later by L. Verlet in
1967. Discrete Newtonian dynamics has the same qualities as Newton’s classical
analytic dynamics. Verlet also published a first-order expression for the
instant temperature which is inaccurate but presumably used in most MD
simulations. One of the motivations for the present article is to correct this
unnecessary inaccuracy in $NVT$ MD dynamics. Another motivation is to derive
simple algorithms for the Nosé-Hoover $NVT$ dynamics (NH) with the correct
temperature constraint. The simulations with NH discrete Newtonian dynamics
show that the NH works excellent for a wide range of the response time $\tau$
of the NH thermostat, but NH simulations favor a choice of a short response
time, even shorter than the discrete time increment $\delta t$ used in MD, to
avoid large oscillations of the temperature.

##  I Introduction

Newton published PHILOSOPHIÆ NATURALIS PRINCIPIA MATHEMATICA ($Principia$) in
1687 Newton1687 ; Toxvaerd2023  and he began $Principia$ by postulating the
discrete Newtonian dynamics in Proposition I, where he used the algorithm for
discrete dynamics to derive his second law for classical mechanics. The
classical analytic dynamics of $N$ interacting objects is the exact solution,
obtained by solving the coupled second-order analytic differential equations
for the accelerations of the objects. The exact analytic solution for the
dynamical evolution of the positions of the objects is characterized by being
time-reversible and symplectic, and by having three invariances for a
conservative system: the conservation of momentum, angular momentum, and the
energy of the system. But Newton’s discrete dynamics is also time reversible
and symplectic and has the same three invariances Toxvaerd2023  , and if the
exactness of classical mechanics is given by these five qualities, then his
discrete dynamics is also exact. Furthermore, there exists a shadow
Hamiltonian where discrete positions are located on the analytical
trajectories for the shadow Hamiltonian Toxvaerd2023 ; Toxvaerd1994 ;
Toxvaerd2012 . This means that there is no qualitative difference between
Newton’s analytic and discrete dynamics.

The discrete Newtonian algorithm was much later in 1967 rediscovered by Loop
Verlet Verlet1967  and used in his pioneering Molecular Dynamics (MD)
simulation of a system of Lennard-Jones particles. The algorithm has been
reformulated several times and appears under a variety of names, the most
well-known is the Leapfrog version. Almost all MD simulations today in Natural
Science are performed by Newton’s discrete algorithm. A review of the discrete
Newtonian dynamics is given in reference No.2 Toxvaerd2023  with the proofs of
the invariances for a conservative system. The proof of the conservation of
the energy is also given here in the appendix.

Loop Verlet also derived a first-order expression $\textbf{v}_{0,i}(t)$ for
the discrete velocity $\textbf{v}_{i}(t)$ for object No. $i$ at time $t$
obtained by a symmetric Taylor expansion Toxvaerd2023appendix

$\textbf{v}_{i}(t)=\textbf{v}_{0,i}(t)-\frac{1}{6}\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}^{\prime}(t)+\mathcal{O}(\delta
t^{4}\textbf{f}_{i}^{\prime\prime\prime}(t))$ (1)

with

$\textbf{v}_{0,i}(t)=\frac{\textbf{r}_{i}(t+\delta t)-\textbf{r}_{i}(t-\delta
t)}{2\delta t}$ (2)

from its positions $\textbf{r}_{i}(t-\delta t)$ and $\textbf{r}_{i}(t+\delta
t)$, and this expression is used in today’s MD Tildesley ; Frenkel . The
velocity $\textbf{v}_{0,i}(t)$ in Eq. (2) is the first term for the velocity
in Verlet’s forward and backward Taylor expansion of the discrete position
$\textbf{r}(t)$, and if the velocity is approximated by this expression, it
entails a systematic and sometimes significant error in the calculated
temperature and in the heat capacity of a system Toxvaerd2023a . It is the one
of the reasons for this article, where the correct expression for the velocity
and the kinetic energy is used to revise the formulation of the $NVT$ Nosé-
Hoover algorithm (NH) for $NVT$ simulations Hoover1985 . Constant temperature
MD has been simulated for many decades and the method is reviewed in
Harish2021 . Many of the algorithms for discrete dynamics at constant
temperatures are, however, with broken time symmetry and irreversible
dynamics. However, the NH dynamics, which is derived from the steady state
Lioville equation for an extended phase space with friction with a friction
coefficient $\eta(t)$, is time reversible Hoover1985 ; Hoover1991 .

The article starts in Section 2 by presenting the expressions for the
velocities, kinetic energy, and temperature by the use of the energy
invariance derived in the appendix, and the NH algorithms for $NVT$ discrete
dynamics are presented in Section 3. The algorithms are used in Section 4 to
obtain $NVT$ state points for a system of Lennard-Jones particles (LJ).
Section 5 summarizes the results.

##  II Discrete Newtonian dynamics

The classical discrete dynamics between $N$ spherically symmetrical objects
with masses $m^{N}=m_{1},m_{2},..m_{i},..,m_{N}$ and positions
r${}^{N}(t)=$r${}_{1}(t)$, r${}_{2}(t),..,$r${}_{i}(t),..$r${}_{N}(t)$ is
obtained by Newton’s discrete algorithm. Let the force, $\textbf{f}_{i}$ on
object No. $i$ be a sum of pairwise forces $\textbf{f}_{ij}$ between pairs of
objects $i$ and $j$

$\textbf{f}_{i}=\sum_{j\neq i}^{N}\textbf{f}_{ij}.$ (3)

Newton’s discrete algorithm is a symmetrical time-centered difference whereby
the dynamics is time reversible and symplectic. The Verlet formulation of the
algorithm is

$\textbf{r}_{i}(t+\delta t)=2\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t)+\frac{\delta t^{2}}{m_{i}}\textbf{f}_{i}(i).$ (4)

The energy in analytic dynamics is the sum of potential energy
$U(\textbf{r}^{N}(t))$ and kinetic energy $K(t)$, and it is an invariance for
a conservative system with analytic dynamics. There is, however, a problem
with the determination of the kinetic energy in the discrete dynamics since
the velocities at time $t$ are not known. Traditionally one uses Eq. (2) and
the expressions

$\displaystyle K_{0}(t)=\sum_{i}^{N}\frac{1}{2}m_{i}\textbf{v}_{0,i}(t)^{2}$
(5) $\displaystyle E_{0}(t)=U(\textbf{r}^{N}(t))+K_{0}(t)$ (6)

for the velocity, kinetic energy $K(t)$, potential energy
$U(\textbf{r}^{N}(t))$ and energy $E(t)$ in MD. But the total energy $E_{0}$
obtained by using Eq. (6) with $K(t)=K_{0}(t)$ for the kinetic energy
fluctuates with time although it remains constant, averaged over long time
intervals. This is due to the fundamental quality of Newton’s discrete
dynamics, where the positions and momenta appear asynchronous and with a
discrete change in momentum at time $t$. The energy invariance in discrete
Newtonian dynamics (D) is derived in the appendix. The kinetic energy
$K_{D}(t)$ in the time interval $[t-\delta t/2,t+\delta t/2]$ is

$\displaystyle
K_{\textrm{D}}=\frac{1}{2}(K_{\textrm{D}_{+}}+K_{\textrm{D}_{-}})=$
$\displaystyle=\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))^{2}}{\delta
t^{2}}+\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))^{2}}{\delta t^{2}}]$
$\displaystyle=\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}(\textbf{v}_{i}(t+\delta
t/2)^{2}+\textbf{v}_{i}(t-\delta t/2)^{2}).$ (7)

###  II.1 The temperature in discrete Newtonian dynamics

The constant kinetic energies in the time intervals $[t-\delta t/2,t]$ and
$[t,t+\delta t/2]$ in between a force impulse at time $t$ are related
Toxvaerd2013 . It is easy to derive the relation

$\textbf{v}_{0,i}(t)^{2}=\frac{1}{2}\textbf{v}_{i}(t+\delta
t/2)^{2}+\frac{1}{2}\textbf{v}_{i}(t-\delta t/2)^{2}-\frac{1}{4}(\delta
t/m_{i}\textbf{f}_{i}(t))^{2}$ (8)

from Eq. (4) between the first-order expression Eq. (2) for the square of the
velocity at time $t$ and the well-defined expression for the square of the
velocities $\textbf{v}_{i}(t-\delta t/2)$ and $\textbf{v}_{i}(t+\delta t/2)$
in the time intervals $[t-\delta t/2,t]$ and $[t,t+\delta t/2]$. The
corresponding expression for the kinetic energy, $K_{0}(t)$, and the
traditional value for the temperature used in MD simulations

$k_{\textrm{B}}T_{0}(t)=\frac{<2K_{0}(t)>}{N_{f}}$ (9)

used in MD simulations for a system with $N_{f}$ degrees of freedom is less
than the mean kinetic energy, $K_{D}(t)$, and the temperature

$k_{\textrm{B}}T_{\textrm{D}}(t)=\frac{<2K_{D}(t)>}{N_{f}}.$ (10)

In the discrete time interval $[t-\delta t/2,t+\delta t/2]$ the relation is

$K_{0}(t)=K_{D}(t)-\sum_{i}^{N}\frac{1}{8}\delta
t^{2}/m_{i}\textbf{f}_{i}(t)^{2}.$ (11)

The correct temperature $T_{D}$ is bigger than $T_{0}$ and the difference
between $T_{D}$ and $T_{0}$ increases with density, temperature, and the
discrete time step $\delta t$. The difference is relatively small and of the
order of a few percent at state points with low pressure, but it is
significant for state points with high densities, pressure, temperatures, the
strength of the repulsive forces, or for large time increments $\delta t$
Toxvaerd2023a .

##  III  NVT ensemble simulations with the Nosé-Hoover thermostats

Many NVT ensemble simulations in MD are with the Nosé-Hoover thermostat (NH-
MD) where a friction $\eta(t)\textbf{p}_{i}(t)$ acts simultaneously with the
force $\textbf{f}_{i}(t)$ Hoover1985 . The Hamilton formulation of classical
dynamics with the friction is Hoover1985 ; Hoover1991

$\displaystyle\dot{\textbf{r}}_{i}(t)=\textbf{p}_{i}(t)$ (12)
$\displaystyle\dot{\textbf{p}}_{i}(t)=\textbf{f}_{i}(t)/m_{i}-\eta(t)\textbf{p}_{i}(t)$
(13) $\displaystyle\dot{\eta}(t)=\alpha^{-1}(K(t)-K),$ (14)

with a restoring friction $\eta(t)\textbf{p}_{i}(t)$ which constrains the
kinetic energy $K(t)$ to $K=2N_{f}k_{\textrm{B}}T$. The damping factor
$\alpha$ is expressed as

$\alpha=N_{f}\tau/2,$ (15)

where $\tau$ is a characteristic ”response” time of the thermostat
Toxvaerd1991 .

The discrete Newtonian dynamics with the time symmetrical Leapfrog formulation
with NH-MD is

$\displaystyle\textbf{r}_{i}(t+\delta t)=\textbf{r}_{i}(t)+\delta
t\textbf{v}_{i}(t+\delta t/2)$ (16) $\displaystyle\textbf{v}_{i}(t+\delta
t/2)=\textbf{v}_{i}(t-\delta t/2)+\delta t/m_{i}\textbf{f}_{i}(t)$
$\displaystyle-\eta(t)(\textbf{v}_{i}(t+\delta t/2)+\textbf{v}_{i}(t-\delta
t/2))/2,$ (17)

and the discrete friction variable $\eta(t)$ can either be updated from two
sets of discrete values $\eta(t-\delta t),\eta(t)$ by

$\eta(t+\delta t)=\eta(t-\delta t)+\delta t/\tau(T_{D}(t)-T)\textrm{ Method
I},$ (18)

or by

$\eta(t+\delta t)=\eta(t)+\delta t/\tau(T_{D_{+}}(t)-T)\textrm{ Method II},$
(19)

where $k_{\textrm{B}}T_{D_{+}}=\sum_{i}^{N}m_{i}\textbf{v}(t+\delta
t/2)^{2}/N_{f}$. Eq. (18) is the time centred analog to Newton’s discrete
central difference algorithm. The time symmetry and time reversibility is
maintained in the NH-MD discrete dynamics with both methods, and both methods
have been used in NH-MD Toxvaerd1991 ; Hoover1990 .

The NH dynamics can be extended to include more than one thermostat and this
has been used in simulations of systems with organic chain molecules. The
thermostats can equilibrate specific modes in the systems or act on all
$N_{f}$ degrees of freedom Frenkelchap6 . In the present formulation of NH
dynamics and with $n$ thermostats which act with all $N_{f}$ modes

$\displaystyle\eta_{1}(t+\delta t)=\eta_{1}(t-\delta t)+\delta
t/\tau_{1}(T_{D}(t)-T)$ $\displaystyle\eta_{2}(t+\delta t)=\eta_{2}(t-\delta
t)+\delta t/\tau_{2}(T_{D}(t)-T)$ $\displaystyle.........,$ (20)

and with a total friction $\eta=\Sigma_{i}^{n}\eta_{i}/n$.

Eq. (17) can be rearranged to

$\textbf{v}_{i}(t+\delta t/2)=\frac{\textbf{v}_{i}(t-\delta
t/2)(1-\eta(t)/2)+\delta t/m_{i}\textbf{f}_{i}(t)}{1+\eta(t)/2},$ (21)

and the $NVT$ time reversible NH-MD dynamics is obtained by the NH-MD
algorithm: Eqn. (16) and (21) with either Eq. (18), and/or Eq. (19) , and/or
Eq. (20).

Figure 1: NH-MD with $T_{Thermostat}=T_{0}=1.00$ $(NVT_{0}$, old NH-MD), and with $T_{Thermostat}=T_{D}=1.00$ $(NVT_{D}$, present NH-MD). The figure shows five hundred (representative) time steps for the simulations of $10^{6}$ time steps. The temperatures $\bar{T}_{D}$ and the root mean square (rms) fluctuations are given in Table I.  Figure 2: Five evolutions of $T_{D}(t)$ in a LJ system of $N=2000$ particles with $\delta t=0.005$ and with a NH thermostat with $T=1.0$ (Method I, Eq. (18)) together with the corresponding temperature without a thermostat (NVE). (Four NH-evolutions are parallel shifted with succeding $\Delta T=0.1$, and the last evolution with $\tau=1000$ is shifted with $\Delta T=0.2$.) All six simulations are started from the same NVE start configuration with different values of $\tau$. The figure shows the last 50000 time steps for the simulations of $10^{6}$ time steps. The temperatures $\bar{T}_{D}$ and the root mean square (rms) fluctuations are given in Table I.  Table 1:  Temperatures in the LJ system. $\tau$ | $T_{D}\pm rms$ | $T_{0}\pm rms$ | $\eta\pm rms$ | Method
---|---|---|---|---
- | 1.000000000 $\pm 9.5\times 10^{-3}$ | - | NVE |
0.0001 | 0.999999993 $\pm 1.7\times 10^{-2}$ | 0.999014 $\pm 1.7\times 10^{-2}$ | $-2.06\times 10^{-6}\pm 0.18$ | I
0.0001 | 1.000000025 $\pm 1.8\times 10^{-2}$ | 0.999023 $\pm 1.8\times 10^{-2}$ | $1.7\times 10^{-6}\pm 1.4\times 10^{-2}$ | II
0.001 | 1.000000970$\pm 1.7\times 10^{-2}$ | 0.999025 $\pm 1.7\times 10^{-2}$ | $-2.81\times 10^{-6}\pm 2.1\times 10^{-3}$ | I
0.01 | 1.000000052$\pm 1.8\times 10^{-2}$ | 0.999024 $\pm 1.8\times 10^{-2}$ | $2.39\times 10^{-6}\pm 2.1\times 0.1277$ | I
0.1 | 0.999998637$\pm 1.8\times 10^{-2}$ | 0.999023 $\pm 1.8\times 10^{-2}$ | $-1.68\times 10^{-6}\pm 4.0\times 10^{-2}$ | I
1 | 1.000005194 $\pm 1.8\times 10^{-2}$ | 0.999030 $\pm 1.8\times 10^{-2}$ | $3.76\times 10^{-7}\pm 1.3\times 10^{-2}$ | I
1 | 0.999995501 $\pm 1.9\times 10^{-2}$ | 0.999029 $\pm 1.9\times 10^{-2}$ | $-9.16\times 10^{-7}\pm 3.2\times 10^{-3}$ | II
10 | 0.999986336 $\pm 1.9\times 10^{-2}$ | 0.999011 $\pm 1.7\times 10^{-2}$ | $1.29\times 10^{-6}\pm 4.3\times 10^{-3}$ | I
100 | 0.999883749 $\pm 2.9\times 10^{-2}$ | 0.998908 $\pm 2.9\times 10^{-2}$ | $2.3\times 10^{-6}\pm 2.3\times 10^{-3}$ | I
1000 | 0.999458904 $\pm 0.12108$ | 0.998477 $\pm 0.12083$ | $1.73\times 10^{-5}\pm 3.2\times 10^{-3}$ | I
1000 | 0.999516366 $\pm 0.1236$ |  | $1.50\times 10^{-5}\pm 3.2\times 10^{-3}$ | II
Figure 3:  The temperature fluctuations of $T_{D}(t)$ in the temperature
interval [0.95, 1.05] for NH-dynamics with one and with two thermostats. The
fluctuations of the temperature with one thermostat and $\tau=0.0001$ are also
shown in green in Figure 1. The fluctuations in the temperature interval
[0.95, 1.05] with two thermostates are shown above, and for two sets of
frictions. With red is for $\tau_{1}=0.0001$ and $\tau_{2}=0.1$. With blue is
$T_{D}(t)$ with $\tau_{1}=\tau_{2}=0.0001$ with different start values of
$\eta_{1}(0),\eta_{1}(-\delta t)$ and $\eta_{2}(0),\eta_{2}(-\delta t)$,
respectively.  Figure 4: Temperature evolution in the LJ system by changing
the thermostat temperature from $T_{1}=1.00$ to $T_{2}=2.00$. With blue is
$T_{D}(t)$ for a relative long response time $\tau=1.00$ and with red is for a
short response time $\tau=0.0001$. The lower inset shows the temperature
evolution at the start of the calibration and the upper inset is at the end of
the calibration of the temperature.

##  IV $NVT$ simulations

The Nosé-Hoover friction $\eta(t)$ thermostats the system’s instantaneous
kinetic energy and temperature $T(t)$ using kinetic energy, Eq. (18) or
Eq.(19). Most $NVT_{0}$ simulations with NH-MD use Verlet’s expression
$T_{0}(t)$, Eq. (9), for the temperature. The exact expression for the kinetic
energy and the instant temperature $T_{D}(t)$ by the discrete dynamics is,
however systematically higher and given by Eq. (10). The difference between
the two canonical simulations, $NVT_{0}$, and $NVT_{D}$, increases with the
discrete time increment $\delta t$ and the force according to Eq. (8). It is
of the order a few percent or less for many MD $NVE$ simulations at moderate
densities and force fields Toxvaerd2023a , but can be much larger for stronger
forces.

The differences between the two ensemble simulations, $NVT_{0}$ and $NVT_{D}$,
where the instant temperature is either given by $T_{0}(t)$ or by $T_{D}(t)$,
are determined for a Lennard-Jones system (LJ). The Lennard-Jones system is a
simple model of a system with the classical gas-liquid-solid thermodynamic
behavior. Systems of $N=2000$ LJ particles were simulated for various state
points $\rho,T$. The differences between $NVT_{0}$ and $NVT_{D}$ at the high
density $[\rho=1.40,T=1.00]$ fluid state point, and with $\delta t=0.010$ and
$\tau=0.0010$ units  are shown in Figure 1. The NH thermostat constrains a
specified kinetic energy either to a given temperature $T$ by Eq. (9):
$NVT_{0}$ with $T_{0}(t)=T$, or by Eq. (10): $NVT_{D}$ with $T_{D}(t)=T$. The
mean canonical temperatures $\bar{T}_{0}$ and $\bar{T}_{D}$, respectively, are
in both cases equal to the input temperature $T=1$, and in both cases the
difference between $NVT_{0}$ and $NVT_{D}$, obtained from $10^{6}$ time steps
are $\approx$ 4- 5 per cent. For $NVT_{0}$: $T_{0}=1.0000000$, green curve in
Figure 1; $T_{D}=1.049$, blue curve. For $NVT_{D}$: $T_{D}=1.0000000$, red
curve in Figure 1; $T_{0}=0.954$, magenta curve. (The figure shows the
difference in a short time interval of five hundred steps).

###  IV.1  Nosé-Hoover simulations of a Lennard-Jones system

The two NH algorithms for $NVT$ dynamics ( I and II) are tested by MD
simulations of a LJ system with the thermostats at various state points in the
fluid state of the LJ system. The first examples are for a liquid state point
$(T,\rho)=(1,0.60)$ near the coexisting liquid with the density
$\rho_{l}=0.589$ at $T=1.00$ Watanabe2012 . The temperatures $T_{D}(t)$ for
various values of the NH-response time $\tau$ are shown in Figure 2, and the
mean temperatures and their root mean fluctuations are given in Table I. The
thermostats work excellent for a wide range of the values of the response time
$\tau$ but the system’s instant temperature $T_{D}(t)$ exhibits oscillations
for values of the response time $\tau>>1$, whereas the thermostat works
excellent for smaller values, and even for response times $\tau<\delta t$.

Table I contains also data for NH dynamics with Method II for $\tau=0.0001,1$
and 1000, respectively, and there are no differences between the mean
temperatures obtained by the two algorithms. The thermostats were tested for
different temperatures and densities, including a fluid system at the high
density $\rho=1.40$ and the NH dynamics work excellently for both methods. The
table for the $NVT_{D}$ with $T_{D}=1.00$ also lists the mean values
$\bar{T}_{0}$. They deviate marginally from $\bar{T}_{D}$ at this state point
with pressure $P\approx 0$ near a liquid in coexisting with gas where the LJ
particles are near the minimum potential energy state, and where the forces
are $\approx 0$.

The instant temperature $T_{D}(t)$ with $NVT_{D}$ dynamics, and for the
smallest value of $\tau=0.0001$ is shown in green in Figure 2. The
temperatures exhibit some ”zones” with different amplitudes of the
fluctuations. The observed ”zones” with a single thermostat friction can be
caused by a nonergodicity of the NH-dynamics, and it can be removed by an
additional friction Martyna1992 . The simulations of the LJ system at
$T,\rho=1,0.6$ with two independent NH thermostats are shown in Figure 3. The
fluctuation zones are removed by the use of two independent thermostats.
However, extended investigations of NH-MD of systems with this nonergoditisity
show, that there are systems for which it is not possible to achieve
ergodisity within computationally achievable time by the use of two NH
thermostates Patra2014 .

A typical situation in NH-MD is a calibration of a system from a temperature
$T_{D}(0)=T_{1}$ to another temperature $T_{2}$. The temperature evolution in
the NH thermostated LJ system by changing the thermostat temperature from
$T_{1}=1.00$ to $T_{2}=2.00$ is shown in Figure 4 for two different values of
the response time $\tau$. With blue is $T_{D}(t)$ for for a relative long
response time $\tau=1.00$ and with red is for a short response time
$\tau=0.0001$. The instant temperature oscillates in both cases around the new
target temperature $T_{2}=2.00$, but the system with the short response time
converges rapidly to instant temperatures near $T_{2}$ whereas the temperature
oscillates with large amplitudes for many thousand time steps in the case with
a long response time. The lower inset shows the evolution shortly after the
start of the calibration and the upper inset are for the temperatures after
$\approx$ hundred thousand time steps.

##  V Summary

Almost all MD simulations are with Newtonian discrete dynamics with Newton’s
Verlet- or the equivalent Leapfrog algorithm, and the discrete Newtonian
dynamics has the same qualities as Newton’s classical analytic dynamics
Toxvaerd2023 . Verlet also proposed an approximative expression, Eq. (2), for
the temperature at $t_{n}$ after n time steps Verlet1967 , which presumably is
used in most MD simulations, but the first-order expression is not useful at
positions where an object is exposed to a strong force e.g. at an ellastic
collision. The expression for the velocity at the discrete dynamics has been
improved by including higher order terms in the Taylor expansion, Eq. (1), by
a spline fit to the discrete positions, or by using a higher-order predictor-
corrector algorithm. A. Rahman used a higher-order predictor-corrector
algorithm in the very first MD simulation Rahman1964 , and also the Runge-
Kutta algorithm is occasionally used in $NVT$ MD dynamics Orel1993 . But
almost all $NVT$ dynamics are with Newton’s discrete algorithm, a thermostat,
and by using Eq. (2) for the velocity at $t_{n}$.

The approximative expressions for the velocity and kinetic energy at the time
step $t_{n}$ can be replaced by the exact expressions for velocities and
kinetic energies in the time intervals at $t_{n}$. One of the motivations for
the present article is to correct this unnecessary error in the determination
of the temperature in Newton’s discrete dynamics. Another motivation is to
derive simple algorithms for the Nosé-Hoover $NVT$ dynamics with the correct
kinetic constraint on the discrete Newtonian dynamics.

The simulations with NH-MD discrete Newtonian dynamics show that NH-MD works
excellent for a wide range of the response time $\tau$, but the NH-MD favor a
choice of a short response time, even shorter than the discrete time increment
$\delta t$ in the MD, in order to avoid large oscillations of the temperature.

##  VI Appendix

The energy invariance

The energy invariance, $E_{\textrm{D}}$ in Newton’s discrete dynamics (D) can
be seen by considering the change in kinetic energy, $\delta K_{\textrm{D}},$
potential energy, $\delta U_{\textrm{D}},$ and the work, $W_{\textrm{D}}$ done
by the force in the time interval $[t-\delta t/2,t+\delta t/2].$ The loss in
potential energy, $-\delta U_{\textrm{D}}$ is defined as the work done by the
forces at a move of the positions Goldstein . An expression for the work,
$W_{\textrm{D}}$ done in the time interval by the discrete dynamics from the
position $(\textbf{r}_{i}(t)+(\textbf{r}_{i}(t-\delta t))/2$ at $t-\delta t/2$
to the position $(\textbf{r}_{i}(t+\delta t)+\textbf{r}_{i}(t))/2$ at
$t+\delta t/2$ with the change in position $\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t-\delta t))/2$ is Toxvaerd2023

$-\delta
U_{\textrm{D}}=W_{\textrm{D}}=\sum_{i}^{N}\textbf{f}_{i}(t)(\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t-\delta t))/2.$ (22)

By rewriting Eq. (4) to

$\textbf{r}_{i}(t+\delta t)-\textbf{r}_{i}(t-\delta
t)=2(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t))+\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t),$ (23)

and inserting in Eq. (22) one obtains an expression for the total work in the
time interval

$-\delta
U_{\textrm{D}}=W_{\textrm{D}}=\sum_{i}^{N}[(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))\textbf{f}_{i}(t)+\frac{\delta t^{2}}{2m_{i}}\textbf{f}_{i}^{2}].$ (24)

An expression for the mean kinetic energy $K_{\textrm{D}}$ of the discrete
dynamics in the time interval $[t-\delta t/2,t+\delta t/2]$ is

$\displaystyle K_{\textrm{D}}=K_{\textrm{D}_{+}}+K_{\textrm{D}_{-}}=$
$\displaystyle\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t/2)-\textbf{r}_{i}(t))^{2}}{\delta(t/2)^{2}}+\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t/2))^{2}}{\delta(t/2)^{2}}]$
$\displaystyle=\frac{1}{2}\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))^{2}}{\delta
t^{2}}+\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))^{2}}{\delta t^{2}}].$ (25)

, and with the change in the time interval

$\displaystyle\delta K_{\textrm{D}}=K_{\textrm{D}_{+}}-K_{\textrm{D}_{-}}$
$\displaystyle=\sum_{i}^{N}\frac{1}{2}m_{i}[\frac{\textbf{(}\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t))^{2}}{\delta
t^{2}}-\frac{\textbf{(}\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))^{2}}{\delta t^{2}}].$ (26)

,

By rewriting Eq. (4) to

$\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t)=\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta t)+\frac{\delta
t^{2}}{m_{i}}\textbf{f}_{i}(t)$ (27)

and inserting the squared expression for $\textbf{r}_{i}(t+\delta
t)-\textbf{r}_{i}(t)$ in Eq. (25), the change in kinetic energy is

$\delta K_{\textrm{D}}=\sum_{i}^{N}[(\textbf{r}_{i}(t)-\textbf{r}_{i}(t-\delta
t))\textbf{f}_{i}(t)+\frac{\delta t^{2}}{2m_{i}}\textbf{f}_{i}(t)^{2}].$ (28)

The energy invariance in Newton’s discrete dynamics is expressed by Eqn. (24),
and Eq. (28) as Toxvaerd2023

$\delta E_{\textrm{D}}=\delta U_{\textrm{D}}+\delta K_{\textrm{D}}=0.$ (29)

## Acknowledgements

The cooperation over many years with my dear friend Luis F. Rull is gratefully
acknowledged.

## References

  * (1) Newton, I.,1687, PHILOSOPHIÆ NATURALIS PRINCIPIA MATHEMATICA. LONDINI, Anno MDCLXXXVII.
  * (2) Toxvaerd, S., 2023, Comprehensive Computational Chemistry 3, 329 (Elsevier, Amsterdam, 2023).
  * (3) Toxvaerd, S. 1994 Phys. Rev. E, 50, 2271 .
  * (4) Toxvaerd, S., Heilmann, O. J., and J. Dyre, J. C., 2012 J. Chem. Phys. 136, 224106.
  * (5) Verlet, L., 1967, Phys. Rev., 159, 98.
  * (6) See Reference No. 2, Appendix.
  * (7) Allen, M. P., Tildesley, D. J., 1987, Computer Simulation of Liquids (Oxford Science Publications, Oxford, 1987).
  * (8) D. Frenkel, D., Smit, B., 2023, Understanding Molecular Simulation (Academic, New York, 2023).
  * (9) Toxvaerd, S., 2024, Phys. Rev. E, in press.
  * (10) Hoover, W. G., 1985, Phys. Rev. A, 31, 1695.
  * (11) Harish M. S., and Patra P. K., 2021, Mol. Phys., 47, 701.
  * (12) Hoover, W. G., 1991, Computational Statistical Mechanics (Elsevier, Amsterdam, 1991).
  * (13) Toxvaerd, S., 2013, J. Chem. Phys., 139, 224106.
  * (14) Toxvaerd, S., 1991, Mol. Phys., 72, 159.
  * (15) Holian, B. L., De Groot, A. J., Hoover, W. G, and Hoover, C. G., 1990,  Phys. Rev. A, 41, 3592.
  * (16) Frenkel, D., Smit, B., 2002, Understanding Molecular Simulation Chapter 6 and Appendix L, (Academic, New York, 2002). .
  * (17) Units in MD of LJ systems: lengts in unit of $\sigma$, energy unit $\epsilon$, time unit $t^{*}=\sigma\sqrt{m/\epsilon}$. Temperature is $k_{\textrm{B}}T/\epsilon$. The LJ forces are cutted and shifted at $r_{cut}/\sigma=2.5$, for cut-and shifted forces in MD: Toxvaerd, S., and Dyre, J. C., 2012, J. Chem. Phys. 134, 081102.
  * (18) Watanabe, H., Ito, N., and Hu, C-K., 2012,  J. Chem. Phys., 136, 204102. Mol. Phys., 87, 1117.
  * (19) Martyna, G. L., Klein, M., and Tuckerman, M., 1992,  J. Chem. Phys., 97, 2635.
  * (20) Patra, P. K., Bhattacharya, B., 2014,  Phys. Rev. E, 90, 043304.
  * (21) Rahman, A., 1964, Phys. Rev., 136, A 405.
  * (22) Janežič, D., and Orel, B., 1993, J. Chem. Inf. Comput. Sci., 33, 252.
  * (23) Goldstein, H., 1980, Classical Mechanics,( Addison-Wesley Press Second Ed. 1980), Chap. 1.