l-
. :l4
v.uz? ‘

$523., ..
..... - n.

. a
ham“

2‘:
9.3 an.

.5

 

iii. 1.
1.3... i .
.3 -v . .

 

i .; V ‘ f .1.
. : 1m; .dwzwrv. .1 Spam.» ‘ .wzmumé
£1.71 \mm..w\..e , . . .. . . . ‘ A. , , .

u... Jfit.

This is to certify that the
dissertation entitled

STUDIES OF SOME PROBLEMS RELATED TO ATOMIC
ORDERING, MOLECULAR MOTION AND PAIR
DISTRIBUTION FUNCTION

presented by

VALENTIN A. LEVASHOV

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Physics and Astronomy

 

 

Mgr/home

Major Professor’ 5 Signatbre

28%?

Date

 

MSU is an Affinnative Action/Equal Opportunity Institution

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINB return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

J WE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/CIRC/DateDue.p65~p.15

__ _.—..__

STUDIES OF SOME PROBLEMS RELATED TO
ATOMIC ORDERING, MOLECULAR MOTION
AND
PAIR DISTRIBUTION FUNCTION

By

Valentin A. Levashov

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Physics & Astronomy

2004

 

 

I“ mamtmww\nn 'IF—wm

‘ll' “"1!- "TH“' "—r w

1-— — v-..-

 

 

 

 

 

 

 

 

ABSTRACT
STUDIES OF SOME PROBLEMS RELATED TO ATOMIC ORDERING,
MOLECULAR MOTION AND PAIR DISTRIBUTION FUNCTION
By

Valentin A. Levasl'rov

In this thesis the results of my work on three out of four projects on which I was
working during my Ph. D. under supervision of Prof. M .F . Thorpe are summarized.

The first project was devoted to the study of properties of a model that. was
developed to reproduce the ordering of ions in layered double hydroxides. In the
model two types of positive ions occupy the sites of triangular lattice. The ordering
of ions is assumed to occur due to the long-range Coulomb interaction. The charge
neutrality is provided by the negative background charge, which is assumed to be
the same at every site of the lattice. General properties of the model in 1d and 2d
were studied and the phase diagrams were obtained. The obtained results predict
multiple phase separations in this system of charges that can, in particularly, affect.
the stability of the layered double hydroxides.

Some properties of the atomic pair distribution function (PDF) were studied dur-
ing my work on the second project. Traditionally PDF was used to study atomic
ordering at small distances. while it was assumed that at large distances PDF is fea-
tureless. Puzzled by the observation that PDF calculated for the crystalline Ni does
not. decay at large distances we studied the behavior. in particularly the origin of
decay, of PDF at large distances. The obtained results potentially could be used to
measure the amount of imperfections in crystalline materials and to test instrumental
resolution in X—ray and neutron diffraction experiments.

During my work on the third project we were developing a technique that would

allow accurate calculation of PDF for the flexible molecules. Since quantum me.-

' "AIL—fl

 

 

 

chanical calculations are complicated and computationally demanding in calculations
of PDF for molecules in liquid or gaseous phases, classical methods. like molecular
dynamics are usually emplcwed. Thus. quanttnn mechanical effects, like zero-point
atomic motion, are usually ignored. However. it is necessary to take into account.
the effect of atomic zero—point motion if there is a desire to extract fine structural
details from the PDF. We developed a method that allows incorporation of the ef-
fect of atomic zero-point motion into the results of classical MD simulations without
performing full quantum mechanical calculations. This technique could be used to
correct classically calculated PDFs and thus to achieve better agreen‘ient between

modeled and experimental PDFs.

 

 

To my grrrmlfat/rer, grandmother and mother"

iv

Acknowledgements

First of all, I would like to thank Prof. Michael F. Thorpe for being my thesis
adviser. It was quite an experience working with him. For me, it was the first
time working with somebody for such a long time. Thus, I was able to observe how
another person thinks and works. Of course, I was able to make such observations
in some other situations, but this time when the main objective of the work was "to
think” about some particular topic. these observations were especially pronounced. I
think that due to him I was able to realize on a deeper level what it means to really
understand a problem and what it means to understand obtained results. There were
situations when I would get some results and I would not. really care about what they
meant. His desire to understand the n'ieaning of the results sometimes would really
open new horizons in the problems. I was very often amazed by his ability to feel
what should be the level of agreement l;)etween results obtained in different ways.
i.e. how different curves should fit. Personally. I would often be satisfied with a fit
that could be significantly worse. But a disagreement occurs when son'iething is done
in a wrong way. Thus his feeling for the level of agreement that we should get was
allowing us to find errors in calculations. I think that. in these situations I got rather
Significant lessons. In most cases he would give a good shot at where the error was.
but sometimes he would just say: ”I do not believe in this result. You have to go
and find an error.” Of course, it is very depressing to search for errors when you do
not know what is wrong or why it is wrong. However, when the error was found (in
most cases there were errors) it would make the results look really beautiful. In this
sense I was able to realize that the beauty of the studied problems lies in the details.
In other words. it is hard to see the beauty if one looks at the problem superficially
or if one would just see the final answer.

In order to see the beauty of the correct answer it. is necessary to see beforehand

 

what an incorrect answer looks like. I would say that. in these details all of the
problerns that I was working on were interesting. They were somewhat non-trivial.
Thus, it was often difficult and depressing, but in the same time enjoyable to work on
them. On the other hand, I often was suffering from lack of understanding of what
the problem that we were working on was really about and of what results we wanted
to obtain.

Often, I also had a hard time trying to understand why it was important to work on
this problem and what its scientific significance was. However, my discussions with
different people revealed that it is in principal the situation that occurs in science
relatively often.

Sometimes we were arguing on this and some other issues. As it seems to me now,

It Was very important for me because I was learning to fight for my own opinion. At
the Same time, I guess that it was rather difficult for him to tolerate such an attitude
and I should acknowledge his tolerance and diplomacy on these issues.

While working on different. projects I was communicating with (,lifferent people. I
“’Ould like to acknowledge them in a “problem by problem” fashion.

In connection with the “Coulomb lattice gas” problem I would like to acknowledge
the following people. I would like to thank Prof. .l.B. Parkinson for his explanations
on the details of the program that was used previously to run Monte Carlo simulations.
I had a very joyful collaboration with Prof. B.W. Southern. His calm and deep
a'ttitude made the study of this problem much more systematic. In particular it was
his idea (and realization) to calculate the properties of the 1d chain with interaction
beYOnd nearest neighbors with transfer matrix method. This allowed us to have a
baSically exact solution for this case. Together, we also gained a deeper understanding
of the devils staircase formalism that I would not have been able obtain by myself.
He also employed Bonsall-Maradudin formalism to calculate exactly the energies of

t . . . . . .
he grven configurations in 2d that was also a Significant development. Especially, I

vi

would like to acknowledge his support in preparing the materials for the publication
[2].

Since this paper is presented in this thesis without any significant changes as a
second chapter I should say that some parts of it that are related to the properties of
the linear chain are written by Prof. B.W'. Southern.

I would like to acknowledge Prof. S.D. Mahanti who carefully read the final version
of the paper and with whom I had detailed discussion of this paper.

I also would like to thank Prof. Thomas J. Pinnavaia for reading the introduction
to the paper. Especially, I also would like to acknowledge Dr. M.V. Chubynsky with
whom I had several very important discussions. In particular the idea about the
Phase separation came to my mind during one of these discussions and it was the first
time when I pronounced the words: “there is phase separation”. This conversation
took place at a very in‘iportant moment and if we would not have had it things could
have developed in a completely different way in general.

While working on problems related to the pair distribution function (PDF) I had
long discussions with Valeri G. Petkov and Emil Bozin. In particular, as a result of

One Of such discussions, I had the idea to write a paper on the role of instrumental
r esOlution function and for the first time I calculated PDF up to 200 A. It was the first
tulle when I was surprised that it. does not decay. The importance of this observation
Was realized a bit later due to M .F . Thorpe and Simon J.L. Billinge and after that we
Calculated PDF at distances up to 10000 A. Even before, as I learned a bit later, this
bellavior was observed in the lab of S..I.L. Billinge (I saw a picture of a PDF on the
Screen of E. Bozins computer with non-decaying PDF at 150 A). However, it seems
that there was a belief that eventually it should decay.

At some moment, Prof. Phillip M. Duxbury advised me to look for a cormection
bet:vveen the non-decaying behavior of the PDF in real space and statistics of energy

1eVels in I; space.

vii

 

 

 

 

 

 

 

 

 

., .1
r' \
-. a it
u
l\
.u.«
I."
' V
.
,' t‘
- r,»
. ‘ n O
n .
I ._
7
l...
‘ .S‘
. ‘ b ‘
.“I—IA-
‘7 1“
~..__
'b -
V“
. ., .
~...__
‘\
_.‘ .
“ n
‘,
|I
..
A“
, Y
‘ O
.\‘
r
.
s
I -
_ ‘ .
“\‘H
. .‘
I ‘\.
I ll
._ ,l
I!
I
l V‘.
n' .
-‘N
4
I
»
l.‘ a
‘.
q
,

Due to this advice I later found a. paper [57] and found out that the behavior of
PDF is related to the Gauss circle problem. There were also few other conversations.

I would like to thank Jean Soo Chung for his explanation on how programs for
calculations of the PDF work and Xiangyun Qiu who was operatively providing the
experimental data for me when I needed them. Of course, all this interest in the PDF
here is in many senses due to S.J.L. Billinge with whom I also had useful discussions.

I also had several useful discussions with Prof. Robert I. Cukier and Prof. James
E. Jackson, on some issues related to molecular dynamics simulations and issues
related to molecular modelling in general.

I would like to acknowledge Dr. M. Lei for being my office mate for the last year
and half. He was the first person with whom I was sharing all the ideas that were
CCurling to my mind. He also was helping me. to find little mistakes in my computer
pr Ograms and manage some other computer issues. Often he was just saying that my
Program does not work because in the description of variable I put. a dot. instead of
a Collima. This helped save me lots of time. \Ve in gmreral discussed many different
iSSUes and it made life in the office much more plausible, at least for me.

There are several people who helped me a lot on computer related issues. Very
often I was asking Dr. Tibor Nagy for advice. During my first two years here h’lykyta
Chllbynsky were also helping me. I had also few discussions with Dr. Andrew J.
Rader. Of course, I should thank the computer staff of the Physics Dermrtment.

especially Dr. George J. Perkins, JeHrey R. Bowes, Nicholas J. Kreucher, Margaret
A' Wilke and Robert T. Neary.

I would like to acknowledge the members of my guidance committee: Prof. J.
Bass, Prof. S.J.L. Billinge, Prof. PM. Duxbury, Prof. J.E. Jackson, and Prof. S.D.
NI a—h anti.

I also would like to thank Prof. D. Tomanek for bringing me to MSU.

I would like to express my sincere gratitude to all of those whose names I would not

viii

list here, but whose presence was extremely important during these six years. They
are those whom I was asking for advice, those who were helping me in different life
situations, those with whom I was talking about different. things, those with whom I
was playing tennis and volleyball, scuba diving, playing underwater hockey, dancing,
partying and (after all) drinking in Peanut Barrel. I would like to acknowledge all
of you because you were turning my existence into living. Thank you very much!
There are also a few memories that would determine my behavior in some complex
situations and I am happy that I have them.

I also would like to thank those who were correcting my English (including the
writing center of MSU) in the present thesis.

Finally I would like to acknowledge my mother and other relatives for their care

and belief in me.

J

 

mu

i Owls

gr,

 

Table of Contents

Abstract 1
Acknowledgements V
List of tables xiv
List of figures (Images in this thesis are presented in color) xv
1 Introduction 1

1.0.1 Chapter ‘2: Charged Lattice Gas with a Neutralizing Background 4

1.0.2 Chapter 3: Absence of Decay in the Amplitude of Pair Distri-

bution Functions at Large Distances .............. 5
1.0.3 Chapter 4: Quantum Correction to the h‘lolecular Pair Distri-

bution Function Calculated Classically ............. 6

2 Charged Lattice Gas with a Neutralizing Background 7
2.1 Introduction ................................ 7

2.2 Model ................................... 11

2.3 Monte Carlo l\Iethod for Simulations of the System .......... 14

2.4 Linear Chain ............................... 20

 

 

 

i Xi'nch'

at 121!

2.4.1 Nearest Neighbors Interaction .................. 20

2.4.2 Finite Range of Interaction .................... 25
2.4.3 Infinite Range Coulomb Interaction ............... 30
2.5 Triangular Lattice ............................. 33
2.6 Conclusion ................................. 40

3 Absence of Decay in the Amplitude of Pair Distribution Emotions

at Large Distances 42

3.1 Introduction ................................ 42

3.2 PDF in d-Dimensional Space ....................... 50
3.2.1 Continuous Definition of PDF .................. 50

3.2.2 Definition of PDF Through Bins ................. 52

3.3 Random Case ............................... 54
3.4 Random Case with Excluded Volume .................. 57
3.4.1 General Discussion ........................ 57

3.4.2 Results of Simulations in 2d ................... 60

3.4.3 Results of Simulations in 3d ................... 63

3.5 Crystals .................................. 66

3.5.1 General Discussion ........................ 66

3.5.2 Exact Solution for 1d Crystal .................. 70

3.5.3 Results of Simulations in 2d ................... 72

3.5.4 Results of Simulations in 3d ................... 76

xi

 

 

3‘" 3m; 1.: "was

 

 

3.5.5 Speculations on the Origin of a-Divergence. .......... 78

3.5.6 Order and Disorder. Similarities and Differei‘ices. ....... 80
3.5.7 The Role of the Spherical Geometry. .............. 81
3.6 Perfect S (q) ................................ 84
3.7 Conclusion ................................. 86

Quantum correction to the Pair Distribution Function calculated

classically. 87
4. 1 Introduction ................................ 87
4.2 Pair Distribution Function for a single molecule ............ 91
4.3 An example of l\/‘Iorse potential. The idea of the method. ....... 92

4.3.1 Molecule N2 ............................ 100

4.3.2 Correction to the heat capacity ................. 106
4.4 Vibrations of large molecules ....................... 112

4.5 TINKER package for molecular modelling and the CfiHH molecule . . 116

4.6 Calculation of PDF for CGH 14 molecule ................. 119
4.6.1 Low temperatures ......................... 119
4.6.2 High temperatures ........................ 128

4.6.3 Pair Distribution Function from Molecular Vibrations and Par-

tial PDFs .............................. 133
4.6.4 Correction to the Heat Capacity for C611“ molecule ...... 139
4.7 Conclusion ................................. 1412

xii

 

'2 Summary

‘ V
A ‘\ o-
I“ H ‘
I .4 .
.- viii \
. .
. \‘ V
h]
o
l.
4

statute:

i am

\ h
a .4
A L \

 

5 Summary
5. 1 Charged lattice gas with a neutralizing background ..........
5.2 Behavior of the PDF at large distances .................

5.3 Quantum corrections for the PDF calculated classically ........

Appendices

A Behavior of pdf at large distances
A.1 Random case in d-dimensions ......................
A. 1.1 Ensemble averaging ........................
A.1.2 Integral approach ..........................
A.2 Amorphous case in d-dimensions: integral approach ..........

A.2.1 Difference with the random case from the integral approach:
detailed derivation ........................

A . 3 1d crystal .................................

B Cut-off for the Morse potential

Bibliography

xiii

1.44

116

147

147

149

153

155

157

159

 

4.1

4.2

List of Tables

Comparison of energies at special concentrations obtained in simula—
tions on lattices with sample sizes 18 x 18, in the first row, ~ 30 x 30
in the second row and exact values calculated with method of Bonsall
and Maradudin in the third row. .

The distances between some atoms in Long and Short conformations.

Boltzman weights for Cb- H 14 molecule in Long and Short conformations
at different temperatures.

xiv

118

128

2.2

2.3

2.4

2.6

List of Figures

The sketch of the structure of aluminum substituted nickel layered
double hydroxides Ni1_IA1I(OH)2(CO3)L/2 - yH-ZO ..........

Dependence of the concentration i' on MC sweeps at a fixed values of
the chemical potential )1. at temperature (T/ J) : 0.08 for the triangular
lattice with 30 x 30 sites in the sample. The values of ii are shown near
their corresponding curves .........................

Histograms that show distributions of concentrations for the curves
from F ig.2.2. Data were collected over 10000 MC sweeps. The sharp
peak corresponding to (,u/ J) = 0.315 has height 4000. The peak cor-
responding to (,u/J) : 0.606 has a height 1700. ............

Dependence of the chemical potential )1 on the average concentration
:7: for the case of a linear chain with interaction limited to nearest
neighbors only. The plateaus seen in sinnilations for (T/ J ) : 0.10
correspond to the regions of phase separation. .............

Energy and Helmholtz free energy as functions of the average concen-
tration :i‘ at temperature (T/ J) = 0.05 for a linear chain with interac-
tion limited to nearest neighbors only ...................

Phase diagram for a linear chain with nearest neighbor interactions.
The solid lines indicate the borders of the phase separation regions.
The dashed curve shown for :1: < 1 / 2 only is the spinodal curve. Both
curves were obtained with the transfer matrix method. The points on
the right show the results of the simulations performed on lattices with
different sample sizes ............................

XV

16

19

24

 

2.7

2.8

2.9

2.10

2.11

2.12

2.13

2.14

The linear chain sample of 120 sites with interactions up to fifth neigh-
bors included at if = 3/8 = 0.375 and temperature ("T/J) = 0.01. The
second row from the top is the continuation of the first row and so on.
The picture illustrates the tendency of the system to phase separation:
in the top row a: z 1 / 2, while in the bottom row .1: a: 1 / 3. Separation
cannot really occur due to the same background charge at every site
of the chain that corresponds to .77 = 3/8 .................

Dependence of the energy per site on concentration, at very low tem-
peratures, for the linear chain with interactimi extending up to first,
second, third, fourth, fifth and sixth neighbors. Points from simula—
tions are plotted on top of the exact solid curves obtained with transfer
matrix technique. Dotted curve shows limiting case of infinite range
coulomb interaction. ...........................

Energy per site as a function of concentration .1" for a linear chain
with Coulomb interactions. The solid curve was obtained by knowing
equilibrium structures predicted with the devil’s staircase formalism
on the lattice sample with 1000 sites. Circles give the results of the
MC simulations at low temperature for samples with 180 sites .....

Dependence of the chemical potential ,u. on concentration it for the
linear chain with Coulomb interaction in the region :7: > 1 / 2. The solid
line shows the limit of zero temperature obtained from the energy curve
on Fig.2.9. The inset shows the results from the lattice sample with
960 sites at very low temperatures in the whole range 0 < :i. < 1. . . .

Phase diagram for the linear chain with Coulomb interactions in the
region 17 > 1 / 2 obtained from MC simulations ..............

Energy per site as a function of concentration (E for the triangular
lattice with Coulomb interactions. The solid curve connects the points
obtained by simulations in previous work[1]. Dashed vertical lines give
the position of concentrations that are stable in MG simulations at
very low tei‘nperatures. .........................

Ground state configurations for the triangular lattice with Coulomb
interaction for .r = 1/2, 3/5, 2/3, and 3/4. ...............

Lattice structures obtained from simulations at three different values
of p at (T / J ) = 0.01. Concentrations of black ions are close to :r =
1/2 = 0.5, :r = 3/5 = 0.6, .2: = 6/7 2 0.857 for the pictures from left
to right. ..................................

xvi

26

[O
\J

30

31

34

36

2.15

3.2

3.3

3.5

3.6

3.7

3.8

Dependence of chemical potential ,u on concentration :1." for the trian-
gular lattice with Coulomb interaction. FOr the cases of temperatures
(T/J) = 0.2 and (T/J) = 0.1 the size of the sample was 9 x 9 and no

points are shown since their density is high. ..............

Phase diagram for the triangular lattice with Coulomb interactimr in
the region .r > 1/2 obtained with MC simulations ............

Experimental F (q) and 0(1) of Ni at 15K measured by neutron scat-
tering. ...................................

Only those sites contribute to the value of 0,-(7') that are in a vicinity
of the distance 'r with respect to site 2‘. .................

PDF of Ni at 15K at different distances. The blue solid lines are results
of calculations while the red dotted lines on the two panels on top show
the result of neutron scattering experiment taken from Fig. 3.1 . . . .

PDF for the amorphous Si: no 2 2.4 A, a = 0.1 A. A: with respect to
1 atom, B: averaged over 10 atoms. C: averaged over 20000 atoms. . .

PDF calculated for a square lattice 0 : 010 according to definitions
through bins (blue lines) and Gaussiaus (red lines). Fig.3.1 ......

Scaled PDF function g;,,(‘r) E Gan/(7H),) vs. reduced distance r/a

Distributions of points for two different excluded distances. The left
figure shows distribution that was obtained by choosing excluded dis-
tance in such a way that excluded area covers 1 / 4 of the total area. In
the right figure, excluded area covers 1 / 2 of the total area. ......

Only those points that are in 0 vicinity of the distance r with respect to
site 2' contribute to the value of (3.,(7‘). Excluded volume around every
point reduces the volume accessible to the other points, that may want
to enter the annulus due to fluctuations, and thus reduces the size size
of fluctuations in the number of points ..................

xvii

37

40

46

49

53

 

     
 

‘ i
re! ..‘
(1“. {A ‘.

 

\‘l‘l.f‘ ,
,‘n‘u‘

\ 'v "i ‘ "'
613' ‘l .

\.

\
l“ MCA‘A‘ ..

v
‘-"\v“ I"
“d". hu'

 

  

 

3.9

3.10

3.11

3-12

3-13

3.14

Dependence of g§_,(E) and 9.3,(0) on 0 for configurations with different
excluded volumes. For the curves (with circles) shown, excluded vol—
ume covers 0.05 , 0.10, 0.20 and 0.525 of the total area of the box. The
curve with squares corresponds to the amorphous distribution of sites.
It basically coincides with the curve that has the fraction of excluded
volume 0.525. The curve with triangles corresponds to completely ran-
dom distribution. .............................

Dependence of the slope (linear part) of the curves 9.3,» versus a (see
F ig.3.9) on excluded volume. The value of excluded volume 0.525 is
very close to the limiting accessible value. The inset shows how the
ratio of the asymptotic value (large a/A) of 9.3, to the value of 93., for
the random case depends on excluded volume. Note the scale on the
y-axis. ...................................

Dependence of g§,~(E) and 934(1) on a for configurations with different
excluded volumes. For the curves (with circles) shown, excluded vol-
ume covers 0.02, 0.05 , 0.10, 0.20 and 0.30 of the total volume of the
box. The curve with squares corresponds to the distribution of sites in
amorphous Si. It basically coincides with the curve that has the frac-
tion of excluded volume 0.35 (not shown). The curve with triangles
corresponds to completely random distribution. ............

Dependence of the slope (linear part) of the curves 9.3,- versus 0' (see
F ig.3.11) on excluded volume. The value of excluded volume 0.525 is
very close to the limiting accessible value. The inset shows how the
ratio of the asymptotic value (large a/A) of 93,- to the value of g3, for
the random case depends on excluded volume. .............

Illustration for the Gauss’s circle problem. Sites that lie inside of the
circle of radius 7‘ shown as black filled circles, while sites that lie outside
of the circle of radius r as open circles. Only those sites contribute to
the value of G,( r) that are in a vicinity of the distance ‘r with respect
to site 2'. ..................................

Figure shows the dependence of the ratio g?L(0)/g'f,,(a) on 0. At
small values of a, when different peaks in C(r) do not overlap, this
ratio linearly decreases as 0 increases. As peaks start to overlap (large

values of a), the rate of decay in the average amplitude of PDF decreases.

xviii

62

6:1

67

71

3.15 Inset shows how F(r) depends on r for r E (0; 1000)a, 7' E (9000; 10, 000)u,
and 7' E (99, 000; 100, 000)a. It shows that on average amplitude of os-
cillations in F (r) does not depend on 7”. At least, it is impossible to see
it by eye. The dashed curve on the main figure shows how < 172(1) >
averaged over few intervals of r of length 1000a depends on the posi-
tion of the averaging interval. At the beginning and at the end of every
averaging interval circles are plotted. .................. 73

3.16 Dependencies of g§(L)/g§,(R) on 0//\ for triangular and square lattices.
Triangles and squares represent the results of numerical calculations for
triangular and square lattice respectively. The solid fitting curves are
given by formulas g§(tr2’)/g§,-(R) = 2.483(0/)\)‘1/3 — 18.264(0/)\) for
triangular lattice and 93(sq1‘)/g§,-(R) = 1.849(0/A)*1/314.478(a/)\) for
square lattice. Dashed lines highlight the effect of excluded volume.
Expressions for them are given by the first divergent term in the for-
mulas for the solid curves. They do not contain the second term (that
might be) caused by excluded volume. From the figure, it follows that
the effect of the excluded volume is bigger for triangular lattice than
for square lattice .............................. 74

3. 1 7 Results of numerical simulations for FCC and orthorhombic lattices.
For rectangular lattice b/a = 2 and c/a. = 3. Points represent results
of simulations. Solid fitting curves are given by formulas described in
the text. Dashed lines highlight the effect of excluded volume. Ex-
pressions for them are given by the first divergent term in the formulas
for the solid curves. Thus they do not contain the second term (that
might be) caused by excluded volume. It can be seen from the figure
that the effect of excluded volume is bigger for FCC lattice than for
orthorhombic lattice. ...........................

“J

\l

3- 1 8 Triangular or square PDFs can be defined through triangular or square
densities, i.e. by counting the number of points inside triangular or
square annuluses and dividing this number by the area of the annulus. 82

3-19 Dependence of pm -— p0 on distance r/a. for triangular PDF on square
lattice. Since the amplitudes of oscillation persist as distance increases
the figure demonstrates that the use of the polygons (any polygon can
be divided into triangles) instead of the circle changes the way in which
PDF should be defined. ......................... 83

3-20 Reduced scattering intensity obtained by Fourier transformation from
PDF calculated for the values of lattice parameters similar to those for ‘I
the Ni crystal at 20K and 300K: (1 2 3.5321, 020 = 0.05Aand 0300 = 0.1.4. 84

xix

4.1

4.2

4.3

4.4

4.5

The blue curve shows the Morse potential for the N2 molecule. The
red curve shows the harmonic approximation to the Morse potential
at small values of 1:. The blue and red horizontal lines are the energy
levels for the Morse and harmonic potentials respectively ........

The squares of the wave functions of the Morse potential for the N2
molecule. The little figures on the left from the top to the bottom
show normalized wave functions for n = 0, 1, 2, 3, 7. The figures on the
right from the top to the bottom show normalized wave functions for
n = 10, 20, 30, 40, 50. There are 61 energy levels for the chosen values
of parameters ................................

The Pair Distribution Functions or the probabilities to find the particle
at a given distance in the Morse potential or its harmonic approxima-
tion at temperatures 100K, 1000K, 2000K and 10,000K. The red and
orange curves represent QM (424,425) and GM (413,414) results for
the Morse potential. The blue and green curves represent QM and CM
results for the harmonic approximation to the Morse potential. The
black curves represent CM results for the Morse potential corrected by
the convolution (4.34) with 0'2 obtained from the harmonic approx-

COT‘T

imation to the Morse potential. .....................

The dependencies of 02 on temperature T that were calculated in four
different ways. The red and orange curves show QM (4.24, 4.25, 4.28)
and CM (4.13, 4.14, 4.17) results for the Morse potential. The blue and
green curves show QM (4.9) and CM (4.12) results for the harmonic
approximation to the Morse potential. The black curve shows the CM
results on the Morse potential with correction (4.31) that comes from
the harmonic approximation to the Morse potential ...........

Average energies for a particle in the Morse potential and its harmonic
approximation as a function of temperature. The red and orange curves
represent the QM and CM solution for the Morse potential. The blue
and green (dashed) curves show the QM and CM solutions for the har-
monic approximation to the Morse potential. The black dashed curve
shows the CM solution for the Morse potential with correction (4.31)
that comes from the harmonic approximation to the Morse potential.
It is assumed that U0 = 0. ........................

96

99

101

103

 

4.6

4.7

4.8

4.9

The dependencies of the heat capacities for the particle in the i\-=Iorse
potential and its harmonic approximation on temperature. The red and
orange curves represent QM and CM solutions for the Morse potential.
The blue and green (dashed) curves show the QM and CM solutions
for the harmonic approximation to the Morse potential. The black
dashed curve shows the CM solution for the Morse potential with the
correction (4.31) that comes from the harmonic approximation to the
Morse potential. All these curves were obtained by the differentiation
(4.39) of the corresponding energy curves on Fig.4.5. .........

The sketch of the CGH 14 molecule. The carbon atoms are shown as
black filled circles. The hydrogen atoms are shown as open circles.
The numeration of the atoms coincides with those used in the Fig.4.9
and Fig.4.10.4.11,4.12 ...........................

The geometry of the (76H 1.; molecule in Long conformation .......

The probability for finding a particular pair of atoms at a given distance
at temperatures 150K on the left and at 500K on the right. The blue
curves were obtained from the classical MD trajectories. The red curves
result from the convolution of the blue curves with the corresponding
correction Gaussian. Every pair of atoms has its own 030,..r(T).

The dependencies of the average distances between $01116. pairs of atoms
on temperature obtained from MD simulations. In the beginning of ev~
ery MD run (at a particular temperature) the molecule was always in
equilibrium Long conformation. At low temperatures these distances
only slightly depend on temperature, as molecule remains in the Long
conformation. As temperature increases instantaneous distance be—
tween a pair of atoms can have significantly different values that occur
when the molecule can be in different conformations (or identical atoms
interchange their positions). Thus, as temperature increases, there oc-
cur sharp changes in the average distances for some pairs of atoms.

ooooooooooooooooooooooooooooooooooooooo

The dependencies of a2 on temperature T for some pairs of atoms.
The orange triangles show the results obtained from MD simulations,
the blue and green curves show the QM and CM results obtained from
the eigenfrequencies and eigenvectors of molecular vibrations. The red
circles show the MD results corrected by adding to the orange triangles
the difference between the blue and green curves. ...........

110

117

118

120

121

4.12

4.13

4.14

4.15

4.16

4.17

The dependencies of 02 on temperature T for some pairs of atoms.
The orange triangles Show the results obtained from MD simulations,
the blue and green curves show the QM and CM results obtained from
the eigenfrequencies and eigenvectors of molecular vibrations. The red
circles show the MD results corrected by adding to the orange triangles
the difference between the blue and green curves. ...........

PDF for the CGH 1.; molecule in Long conformation at low tempera-
tures. The brown curves show the results of MD simulations before
the correction is applied. The blue curves show the results of MD
simulations corrected by convolution (4.33, 4.34). ...........

The blue and green solid curves show total corrected PDFs for the
molecule in Long and Short conformations. The red dashed curve
shows Boltzman combination (4.70) of PDFs from Long and Short con-
formations. At all temperatures when the molecules remains in one or
another conformation the combined PDF basically coincides with PDF
in Long conformation due to the values of the Boltzman weights. . . .

The brown curves on both figures show MD results obtained on the
Long (top) and Short (bottom) conformations. The figure on top also
shows the corrected PDFs obtained by the correction of MD results
from the Long conformation with the set of own from the Long (the
blue solid curve) and Short (the red dashed curve) conformations. The
bottom figure also shows the PDFs obtained by the correction of MD
results from the Short conformation with the sets of own. from Short
(blue solid curve) and Long (red dashed curve) conformations.

Corrected PDFs for the C6H14 molecule at different ten'iperatures. The
inset shows the region between 3 A and 6.5 A on a bigger scale.

The panel on top shows PDFs with respect to carbon atom number
1 created by the other carbon atoms only at 298K. Blue and green
curves were obtained from molecule in Long and Short conformations
correspondingly. The widths of peaks were calculated in QM approach
using eigenfrequencies and eigenvectors of molecular vibrations. As
distance increases there are peaks that correspond to the carbon atoms
number 2, 3, 4, 5 and 6. The red curve is a linear combination of
blue and green curves with corresponding Boltzman weights. The red
curve in the bottom panel is the same as the red curve in the tOp
panel. The grey curve, that also shows PDF with respect to the carbon
atom number 1 due to the other carbon atoms, was obtained from MD
simulations. The black curve was obtained from the grey curve using
our correction scheme. ..........................

125

127

130

 

 

 

 

4.18

4-19

4.21

The top panel shows partial PDFS with respect to the carbon atom
number 1 calculated for the molecule in Long conformation from eigen-
frequencies and eigenvectors of molecular vibrations in the QM ap—
proach. The blue curve shows partial PDF due to the other carbon
atoms only (hydrogen atoms are ignored). The red curve is due to
all the other atoms (carbons and hydrogens). Thus, the difference be-
tween two curves is due to the hydrogen atoms. The bottom panel
shows the same curves as the top panel, but for the molecule in Short
conformation. ...............................

The top panel shows partial PDFS with respect to the carbon atom
number 1 created by all other atoms (carbons and hydrogens). The
blue and green curves were calculated from eigenfrequencies and eigen-
vectors of molecular vibrations in Long and Short conformations cor-
respondingly. The red curve is the linear combination of the blue and
green curves with corresponding Boltzman weights. The red curve in
the bottom panel is the same as the red curve in the top panel. The
grey curve was extracted from MD simulations in which transitions
between conformations occur. The black curve was obtained from the
grey curve using our correction method. ................

The blue and green curves in the top panel were calculated from eigen-
frequencies and eigenvectors of molecular vibrations for the molecule
in Long and Short conformations correspondingly. They show the total
PDFS of the molecule. The red curve is a linear combination of the blue
and green curves with corresponding Boltzman weights (wL(T) = 0.79
and wS(T) = 0.21). The red curve in the bottom panel is the same
as the red curve in the top panel. The grey curve, that also shows
the total PDF of the molecule, was obtained from MD simulations.
The black curve was obtained from the grey curve using our correc-
tion / convolution method ..........................

The blue and green curves in the top panel were calculated from eigen-
frequencies and eigenvectors of molecular vibrations for the molecule
in Long and Short conformations correspondingly. They show the total
PDFS of the molecule. The red curve is a linear combination of the blue
and green curves with corresponding Boltzman weights (wL(T) = 0.62
and w3(T) = 0.38). The red curve in the bottom panel is the same
as the red curve in the top panel. The grey curve, that also shows
the total PDF of the molecule, was obtained from MD simulations.
The black curve was obtained from the grey curve using our correc—
tion/convolution method ..........................

xxiii

136

137

138

 

4.22 The dependencies of the energy and the heat capacity on ten‘iperature

for the C611 14 molecule. Blue and green curves represent the QM and
CM results for the harmonic approximation to the real potential. The
orange curves and circles represent the result of the MD sinmlations.
The red curves and circles show the corrected MD results. See more
detailed description of the figure in the text. ..............

Dependencies of (7510»!le on temperature (see Fig.4.4) for different
values of the upper cut off in the integrals (4.14, 4.17). Thus, if T S
10000 (K) any cut-off in the interval (1523.0) Awould lead to basically
the same result. The inset shows how < UMCM(T) > (in temperature
units) depends on temperature. Thus for the potential energy any
upper cut off (see (4.38)) in the interval (1523.0) Ais also suitable.

xxiv

140

158

 

Chapter 1

Introduction

During my Ph.D. years I was working on several different projects all in the area
of con'lputational / theoretical condensed matter physics. Thus, it. is somewhat difficult.
to write an introduction to the thesis because, as it seems to me, it should unite all
these different projects under the roof of one idea. But it is hard to unite things that
are different. It also means that an introduction, if it has to be written, should be
rather general. So let it be...

Any object in this world interacts in one way or another with other objects. In
fact, We know that an object exists from the way it interacts with us or from the way its
footSteps interact with us. Thus, objects are distinguished by their interactions. Due
to t-he interactions things in the world become organized. Thus, there are three words:
E”listence, interaction and organization that from my point of view, are inseparable
in nature and as follows in the essence of all sciences.

Generally speaking, all different sciences addresses the same questions: which
ObJects exist, how they interact with each other and how objects are organized. We
leaIn about the objects from the way they interact with each other and from the
Way they are organized. Almost always when we want to learn about an object.

“'6 artificially put it in conditions that will reveal its ability to interact in one way

 

or another; i.e., reveal its properties. In other words we do experiments. On the

other hand, if we know about interactions, we use this knowledge to construct more
complicated objects or to predict the behavior of more complicated objects in some
new Situation. For example, if the interaction between two masses is known, the
structures of solar systems and galaxies can be predicted.

N'Iay be the last statements are not true with respect to all sciences—there are so
many that I cannot think about all of them. But it seems to be true in sciences like
physics, chemistry, biology and sciences that are directly related to them, like micro—
biology. Mathematics is somewhat different, because at first it states the existence of
different objects, like points, lines, angles. Then it. establishes interactions between
thern - axioms. From these objects and interactions more con'iplicated objects and
relations between them are coristriicted—theorems. Thus, from this point of view,
these three words are also applicable to mathematics.

All objects are immersed in space or space contains all the objects. It also allows
ObleCt-s to interact with each other. Often. while studying properties of different
Objects, we do not. pay attention to the space in which they are immersed. However
properties of the space and properties of the objects in space are as inseparable as
those three words. Thus, by learning about. the objects we also learn about the space.

One of the properties of space that is always involved in the nature of interactions
and motion of objects is the dimensionality of space. Everyday observations of this
World suggest that it is 3-dimensional. Thus, in order to describe the position of an
Object in space we need 3 real numbers: :13, y, 2.. However, often there is no need to
“39 all three coordinates. Thus. in order to find our location in the city we do not need
all three coordinates— two is enough. Because of it most city maps do not contain
infOrmation about the height of a particular city point with respect to sea level. The

properties of complex objects, even if they are constructed from objects with very

SIInple interactions between them, can be rather complicated. These properties can

 

also drastically depend on the dimensionality of space, as in Ising model in which
features of phase transitions are different in different dimensions.

Evolution of sciences is basically the development of ideas and tools that are used
to study objects, to model their properties and to construct new objects.

From this extremely general point of view my thesis can be separated into two
parts.

The first part is devoted to the modeling of properties of an object. that consists
of many objects with simple known interactions between them (Chapter 2).

The second part is about a particular technique that is being used to study the
properties of complex objects: the ordering and interactions of simple oliijects from
whicll these complex objects are constructed. This second part, in its turn, can be
divided into smaller pieces.

Although this technique was used to study ordering of objects and their proper—
ties for almost a hundred years, it turns out that it has some properties, related to
the general properties of space, that never were carefully discussed. The results of
the Study of these properties are presented in chapter 3.

Footsteps of more complicated objects sometimes are rather complicated. Thus,
it Order to extract information about them from their footsteps it might be necessary
to IIlodel the possible footsteps. To the modeling of these footsteps, in the frame of a
pa'I'ticular technique that is used to study objects properties, chapter 4 of this thesis
is devoted.

Since, as it was written in the very beginning, I was working on several almost
indEzpendent projects, it is natural to write introduction for every project indepen-
dently. Thus, every chapter in this thesis is about a particular project, it has its own
introduction that explains the motivation for this research project, describes direc-
t"10113 of the investigation and methods that were used in the study. In the remaining

part of this introduction there are short descriptions of every project/ chapter with

 

the summaries of the main results.

1.0.1 Chapter 2: Charged Lattice Gas with a Neutralizing
Background

In this part of my thesis we consider a model that was first introduced by
Y.Xia.o, M.F. Thorpe and .18. Parkinson in [1] to describe the ordering of two
different types of positive ions in the metal planes of layered double hydroxides
Ni 1 _xAlx(OH)2(C03)x/2 - yHgO. The ordering is assumed to occur due to long-range
Coulomb interactions, and overall charge neutrality is provided by a negative back-
ground representing the hydroxide planes and C03" anions. The previous study [1]
was restricted to the ground state properties. Here we use a l\-"Ionte Carlo technique
to extend the study to finite temperatures. The model predicts that at some values
of the concentration :1), the system can exhibit an instability and phase separate. In
order to evaluate the precision of these Monte Carlo procedures, we first study a linear
Chain with finite ranged interactions where exact scflutions can be obtained using a
tr a«I’lsfer matrix method. For a linear chain with infinite-ranged interactions, we use a
devil’s staircase formalism to obtain the dependence of the energy of the equilibrium
COlrlf‘lgurations on at. Finally we study the two dimensional triangular lattice using the
same Monte Carlo techniques. In spite of its simplicity, the model predicts multiple
first order phase transitions. The model can be useful in applications such as the
H10Cleling of the ordering of intercalated metal ions in positive electrodes of lithium
baL’Cteries or in graphite. The obtained results were published in PhysRevB [2] and

the text in this chapter is basically identical to the text in this publication.

 

1.0.2 Chapter 3: Absence of Decay in the Amplitude of Pair

Distribution Functions at Large Distances

In this chapter the behavior of pair distribution function (PDF) at large dis-
tances is addressed. Traditionally the PDF, that can be obtained from the Fourier
transform of powder diffraction data, was used to describe short-range correlations
in atomic positions. Amplitudes of peaks in experimental PDF decay as distance in-
creases. Thus, it was always assumed that the PDF at large distances is featureless.
Puzzled by the observation that PDF calculated for crystalline materials does not de-
cay at large distances, if instrumental resolution is ignored, we studied the behavior
of PDF at very large distances. To the best. of our knowledge, the origin of the PDF
decay at large distances has never been carefully discussed. It turns out, surprisingly,
that the increase in the number of neighbors at large distances does not lead to the
decay of the PDF independently from the type of the material. In other words, the
PDF calculated with respect to one atom does not decay at large distances not only
for the crystalline, but also for amorphous materials. We find that this behavior in
aInOl‘phous materials is caused by random fluctuations in the radial number density.
ThUS PDF in amorphous materials decays due to ensemble averaging over different
centI‘al atoms. We achieve an accurate qualitative description of fluctuations for the
(338 when atoms are distributed randomly in space. Differences with the amorphous
Case are discussed. The case of perfect, single component, crystals is the most inter-
esting because in it all atomic positions are equivalent and there is no need to average
Over different atoms. Thus total measurable PDF for perfect crystals does not decay
at, large distances if instrumental resolution is ignored. However, the origin of this
behavior in crystals is significantly more complicated. It turns out. that this behavior
of the PDF is related to the still unsolved problem that C.F. Gauss formulated more
than a hundred fifty years ago and that has give rise to the whole area in mathematics

Qfilled lattice point theory. Further investigation of this case is obviously needed. For

generality we discuss the case of the PDF in d-dimensional space. Our results can
be used to measure the amount of dislocations in crystalline materials and to test

instrumental resolution in scattering experin'ients.

1.0.3 Chapter 4: Quantum Correction to the Molecular Pair

Distribution Function Calculated Classically

In this chapter we present some ideas and the developed technique that may
significantly improve precision of PDF calculations for complex molecules. This tech-
nique allows, in particularly, the incorporation of the purely quantum effect of zero-
point motion into the pair distribution function calculated classically for molecules
using Monte Carlo or Molecular Dynamics simulations. This correction may sig-
nificantly improve agreement. between modeled and experimental data, esl‘)ecially at
small distances. Thus it may allow more definite conclusions about inter- and intra-
molecular motion, including flexilnlity, and also about mutual orientations of different

molecules.

Chapter 2

Charged Lattice Gas with a

Neutralizing Background

2.1 Introduction

This chapter is devoted to the study of a model that was introduced earlier[1]
by Y.Xiao, M.F. Thorpe and J .B. Parkinson, to describe the possible ordering of
metal ions that can occur in aluminum substituted nickel layered double hydroxides
Ni1_IAlx(OH)2(CO3)x/g - yHgO[3, 4, 5, 6, 7].

The Ni ions in Ni(OH)2 occupy the octahedral holes between alternate pairs of
OH planes and thus form a triangular lattice identical to that adopted by the OH
ions, as shown on Fig.2.1. The two planes of OH ions, with the plane of Ni atoms
between them, form a brucite like layer of the host structure[8].

Ni(OH)2 can exist in two polymorphous crystal structures denoted as a and ,3.
Both structures consist of brucite—like layers, that are well ordered in the ,8~phase
and randomly stacked in the (Jr-phase. The interlayer spacing (gallery) in the a-phase
usually is significantly larger than in the fi-phase due to the large number of water

molecules and anionic species that. can penetrate into the gallerieslgl.

 

 

 

 

 

 

 

 

 

Figure'Z .:1 The sket< h (1fthe strurtme of aluminum substituted nickel layered double
hydroxides Ni1__,.Al (CH) (C () )r/v 'Nll2()

Nickel hydroxide Nl(()ll)g in the ,d-phase has li1een extensively used as a. material
for the positive electrode in rechargeable alkaline batteriesllO]. However it has been
shown that electrodes based on the (1- -phase hydroxide have a bigger charge capacity,
lower charge and higher discharge voltages [11, 12]. Unfortunately the (,r-phase reverts
to the fi-phase in the alkaline 111edia (KOH-for example) which is used in batteries.
Thus the stabilization of the c1-phase of N i(OH)2 in an alkaline media is an important
goal for potential applications. To enhance the Ni(()ll)2 stability many studies of
the partial substitution of metal ions (Al-for example) for Ni in the lattice of nickel
hydroxide have been carried out [153, H. 15, .16].

When Nig’ ions are sul’1stituted by Ali” ions in the metal sheets, [CO3]2 a11-
ions accumulate in the galleries in such an amount that total charge neutrality
is preserved. The amount of water in the galleries depends on the preparation

method: the general formula of the aluminum substituted layered nickel. hydroxide is

Nil—xAlx(OH)2(COB)x/2 'yHQO-

It is natural to expect that the stability of these compounds is composition
dependent and also depends on the preparation technique. Different authors were
able to synthesize layered hydroxides with different concentrations of aluminum. In
particular the range 0 S :1: S 0.4 has been reported [8]. Not all ranges of composition
:1: are accessible due to the limited number of [CO3l2‘ ions and water molecules that
can penetrate into the gallery.

Possible orderings of the Ni and Al ions in the metal planes can also affect the
stability of the compound. Several authors have reported observations of in-plane
ordering of metal ions [3, 4, 5, 6, 7]. Ordering of ions was observed near the values
of :1: equal to 1/4 and 1/3 that are in registry with the host geometry of triangular
lattice. However, it seems that there are no detailed experimental studies of the
ordering as a function of composition at.

The ordering of metal ions in alloys is often considered within the framework of
a lattice gas model where only interactions between neighbors that are not separated
by large distances is taken into acc01.111t because of the relatively short screening
length[l7, 18, 19] caused by free electrons.

It is generally accepted that in layered hydroxides the Coulomb interaction be—
tween positively charged metal planes, negatively charged hydroxide planes and neg-
atively charged anions [COgl2‘ in the galleries are important. The screening length
in layered hydroxides should be significantly larger than in metal alloys, because di-
electric screening caused by water and other polar molecules is weaker than screening
caused by free electrons. Thus a model that takes into account long range Coulomb
interaction and interaction between positive metal ions in the plane with negative
hydroxides layers and negative anions in the galleries might be more suitable than a
lattice gas model to describe the ordering of metal ions in layered double hydroxides.

The role of ordering due to Coulomb interactions has been discussed previously by

Thompson [20].

We have previously suggested a simple model to describe the ordering of metal
ions in layered hydroxidesfl]. In this model two kinds of positively charged metal
ions occupy the sites of a triangular lattice. The lattice is immersed in a negatively
charged background which represents the hydroxide layers and negative anions in the
galleries. It was assumed that the background charge is the same at every site of the
triangular lattice. Thus the total charge at every site is formed by the positive charge
due to the metal ion and the negative background charge. The interaction potential
between sites was assumed to be a long ranged 1/7' Coulomb type.

In the previous work[1] the dependence of the ground state energy of this model
system on the concentration of Al was studied assuming a homogeneous concentration
of metal ions in the plane. Equilibrium ordering configurations of ions that can occur
at each concentration in the range 0 S .17 g 1 were calculated and compared with
corresponding X—ray diffraction patterns.

In this study a new interpretation of the previous results is suggested. It will be
shown that at some concentrations .17 the system is unstable with respect to phase sep-
aration into phases with concentrations 1‘1 and 5172 such that. 1'1 < :r < .rg. The phase
diagram of the system is calculated in the (T, 1') plane using the grand canonical
ensemble by introducing a chemical potential n. In the case of the layered hydrox—
ides, the chemical potential under consideration is not related to the voltage on the
electrodes and represents only a useful way to obtain the phase diagram.

The model is quite general and can be employed to describe ordering and first
order phase transitions in ionic systems with long range interactions. It may have
some application to the ordering of intercalated Li ions in rechargeable Li-batteries
[20, 21, 22]. Predicted phase separations can lead to the staging when homogenous
planes with different concentrations of metal ions will form. In plane long range

interaction in this case can be similar to that occuring in staged graphite intercalation

10

 

 

compounds[23, 24].

The chapter is organized as follows. In section 2.2 the model is defined. Then in
section 2.3 the discussion of the details of the Monte Carlo (MC) method that was
used to obtain the phase diagrams is presented. Section 2.4 describes an application
of the method to the linear chain. Subsections 2.4.1 and 2.4.2 are devoted to the case
of finite ranges of interaction where an exact solution can be obtained using t'ra'nsfw
matrix techniques. The case of a linear chain with infinite range Coulomb interaction,
in which the energies of equilibrium configurations can be calculated exactly using
the devil’s staircase method, is discussed in subsection 2.4.3. Finally in section 2.5
the case of the two dimensional triangular lattice is studied. Results are summarized

in the conclusion.

2.2 Model

Consider a system composed of two types of positive ions which occupy the sites
of some lattice. Every site of the lattice is occupied either by a black ion with charge
Q, or by a white ion with charge Qw. The concentration of black ions is .L‘ and the
concentration of white ions is l — .r. In addition to these two types of positive ions
there is also a negative compensating uniform background Charge q, at every site of
the lattice, that ensures charge neutrality in the system. Hence at any site i a total
charge is equal to either Ab = (2;, + q or Aw = Qw + q. The Hamiltonian of the system

of charges can be written in terms of the pairwise interactions V01

H: Ewan). (2.1)

i<j

For the Coulomb interaction V2) = l/Rdj, where Rij is the distance between sites 2'.
and j. In the following the case of truncated Coulomb interaction, for which there is

interaction only between neighbors that are in some range, also will be considered.

11

The value of the background charge depends on the values of Qb, Q“, and on the

concentrations of those ions. The average charge per site due to the positive ions is
(f: J‘Qb+(1—$)Qw- (22)

The value of the charge that provides the uniform background with overall charge
neutrality is then given by q = ——q. Thus Ag, = (2;, + q = ((2,, — Q,,,)(1 — at) and
Aw = Qw + q = (Q, — Q,,.)(—.r). Introducing JU : VU-(Qb — Qw)2 we rewrite Eq.

(2. l) as:
H = 2.1.,(72, — .1:)(n, — .1:). (2.3)

K]
Where n,- = 1 if site '1'. is occupied by a black ion and n, = 0 if site '1' is occupied by
white ion. Thus the precise values of the positive charges do not cl1ange the general

Properties of the model. From Eq. (2.3) follows that at a fixed concentration 1‘

H = 2 1,112,711 — .172 Z Jij. (2.4)

-i<j i<j

ThLIS at fixed concentration the model is essentially a lattice gas model because in
this case the second term in Eq. (2.4) is the same for all configurations of ions. This
term is due to the neutralizing background and cancels the divergence of the energy
a‘SSOCiated with the Coulomb repulsion when the interactions are of infinite range.
It is easy to estimate the values of the constants Jij if the lattice constant and
Cllal‘ges Q, and Qw are known. For the aluminium substituted nickel layered double

h~y'dl‘oxide a. = 3.032 2M8], Qb = 3e and Q“, 2 26. For the nearest neighbors we have:

‘)

)-

J : l/nn(Qb _ Q10)? : _ ' (25)

EU.

Tllus if we use a dielectric constant appropriate for water 6 = 80, then J 2 60 meV

0 . .
1‘ about 500 K. Short range screening, due to the presence of water and other ions,

12

 

can effectively decrease the values of the J0.

In order to use the grand canonical ensemble to study the equilibrium properties
of this system of charges a chemical potential term flax, n, should be added to
the Hamiltonian. Now the concentration I can fluctuate and hence the background

Charge will also fluctuate. The grand canonical Hamiltonian can be written as:

H = Z J,J(n,- + Ql('71'j+ (I) — 111:": a (26)
<ziJ> i
where q describes the uniform backgrormd charge which is adjusted to be equal to
q = - < Tl," >-_— —.’17.

Although the primary goal is to study the two-dimensional triangular lattice with
infinite range Coulomb interactions, it is worth at first to consider the linear chain as
an example to gain a better understanding of the model since there are exact analyt-
ical methods that can be used in two limiting cases: finite range interactions can be
Studied exactly using transfer matrix methods [25] and infinite range Coulomb inter—
actions can be described in terms of a Devil’s staircase formalism[17]. The study of
the linear chain will also give an insight into the precision of the numerical techniques
that will be used for the triangular lattice with Coulomb interactions.

In the following sections the following questions are addressed. What. is the
equilibrium structure of the charges for a given concentration of the ions? How
does the equilibrium energy of the system depend on the concentration of the ions?
HOW does the chemical potential depend on the average concentrations of the ions
at different temperatures. Do phase transitions occur in the system and what is the

131—1886? diagram of the system at finite temperatures?

13

2.3 Monte Carlo Method for Simulations of the

System

In Monte Carlo simulations a l\=’Ietropolis algorithm[26] will be used to accept or
re 3' ect elementary moves that will be performed on the charges to bring them into an
equilibrium configuration. Two different types of moves will be used: A interchange
the positions of black and white ions in the lattice. B change the color of the ion at.
a particular site.

The simulations can be carried out with either constant values of the concen-
tration or with constant values of the chemical potential. In case of simulations at
a constant value of concentration only moves of type .4 were used. In case of simu-
lations at a constant value of the chemical potential both types of moves A and B
were used. The A—move does not change the value of the background charge since
it does not change the concentration. The B—move changes the concentration and
thus the background charge has to be changed at every site in the lattice. In order to
decide Whether to accept or reject the move it is necessary to calculate the energy of
the system before and after the move. A lattice of size. N in case of the linear chain
and N x M for the triangular lattice are considered. Standard periodic boundary
Conditions are applied with respect to this central zone. Thus the central zone occurs
in the center of bigger lattice surrounded by sucroumlt'ny zones. In order to calculate
the energy of the system we calculate the energy of interaction between all sites inside
the central zone and the energy of interaction between the sites in the central zone
With sites in all surrounding zones. Since the system is charge neutral, the contribu—

tion to the energy from surrounding zones that are far away from the central zone
are much Smaller than the contribution from surrounding zones that are close to the
central Zone. In fact, it was found that if N and M are of the order of 10, then it. is

e . . . . .
nough t0 Consrder the lattice of s1ze 5N x SM in order to calculate the energy With

14

 

sufficient precision for almost all crmcentrations. In other parts of this chapter we
refer to the size of the central zone as to the sample size. with the periodic boundary
conditions described above.

The total energy of the system given by Eq. (2.3) can be separated into three
parts which represent interaction between bla.ck-l_)lack, white-white and black-white
sites

E : Ebb + Eww + Ebw - (27)
It follows from Eq. (2.3) that Ebb can be written as
,, 2 ,,b,.b
Ebb : (I — .1.) E Jamil”. . (2.8)
K}
In the sum above, the index '27 runs over all sites in the central zone and the index j
over all sites in the central and su'l'romrdilng zones (1' 7é j). The quantity a? is unity if

a black ion occupies site '2'. and zero otherwise. Using the notation obb(j) = 2} Jun?

EQ- (2.8) can be rewritten as:

Ebb = (l — gr)2 271:)01)h('ll : (1 — M20“).

i

am

In the same manner we can write EM. = 1:20“, and Eb", = 2.1?(1 — I)ob,,,. Then the

energy Of the system per site can be written as:

I o
[(1 — @2051, + 2:1:(1 — .T)Obw + Tamw] , (2.10)

aNxAn

\Vhere the first term in square brackets comes from the interaction between black-

blaCk siteS, the second from black-white sites and the third from white—white sites.

It IS easy to see from Eq. (2.10) that the equilibrium energy of the system is a

Symmetric function, E(;17) = E(1 — :r), with respect to a: = 1 / 2. The ground states

COI‘I‘eS . . .
pondlng to concentrations :r and 1 — a? can be obtamed from each other by

15

 

 

 

 

 

 

 

j I I I I I I I I I I I I I I I I r I I I T I I IMI I T I I I F? I I I I I T? I I I I I I I I ..

l .

(ii/J) = 0.606 E

,

.R 0.9 1. . ’ w a;

'5‘ [l ' (T/J) = 0.08 i

‘5; :

a _

5 0 8 1

g (u/J) = 0.375 3 F3“-

0 ,r, . _ . .
U i i 1 l
0.7 __ .
(ii/J) = 0.315 i

0 6 I l l I I 1 L 1 1 l_ 1 1 1 1 1 l 1 14 I l l l l l l l l l l 1 1 1 1 1 1 I 1 1 l l l 1 1 1 J 1 l 1:

° 0 2000 4000 6000 8000 10000
MC sweep

Figure 2.2: Dependence of the concentration i' on MC sweeps at a fixed values of
the chemical potential ,u at temperature (T / J ) = 0.08 for the triangular lattice with
30 X 30 sites in the sample. The values of ii are shown near their corresponding curves.

changing all white sites into black sites and all black sites into white sites. In this
case a: H 1 — a: , obb <—+ o,,,,,., 0m H 0b”. and it follows from Eq. (2.10) that the energy
remains the same. In order to accept or reject the move we calculate the change in
AE — MAE. For calculation of AE it is necessary to take into account that if we turn
a. White site 2' into a black site then charges at the sites corresponding to the site i'.
but situated in the surrounding zones should be changed also. When an elementary
move is performed, then values of 01,5, 05“,. am can be updated by calculating the
Sums Ubb(2'), Ubw(‘l), ou,u,(i) for the particular site «i that participate in the move. This
Significantly reduces the calculation time since it is not necessary to recalculate the
energy 0f the whole lattice again after every move.
We say that one MC step was performed if one attempt to perform operation A
or B was made. We say that one MC sweep was made if as many MC steps were

ade as there are Sites in the sample. In Simulations at a constant p. we initially tried

16

to vary the frequency with which operations A and B were performed, but we found
that 1 : 1 ratio was close to the optimum value. For every value of the concentration
or the chemical potential, simulations start at a relatively high temperature T 2 J.
If simulations are to be performed at constant a? then in the initial configuration
black and white ions in amounts corresponding to .r. are randomly distributed over
the lattice sites. If simulations are to be performed at a constant value of ii. the
initial configuration is less ii‘nportant.

We used the following criteria to check the equilibration of the system at a
given temperature. Let B2 will be the average value of the energy in the last 10
IVIC sweeps and B1 will be the average value of the energy in the previous 10 MC
sweeps. Let 052 and orgl be the average fluctuations of energy in those two cycles. If
[E2 — E1] 3 % min(oE,,o'E.2) then we say that. the system is sufficiently equilibrated
in order to collect the data. If this condition is not fulfilled another 10 MC sweeps
are made until this condition is met and so on.

After the equilibration, in order to obtain statistics, we calculated and stored
the values of parameters of interest after every MC sweep. Their convergence to
equilibrium values was verified by plotting them versus the number of the MC sweep.
The number of required MC sweeps varied depending on the size of the system, type
Of interaction and temperature. When the necessary data at a temperature T were
Collected, the temperature was decreased by a small amount 0T. For smaller values
Of the temperature T, a smaller value of (ST was used.

It will be shown below that first order phase trai'isitions occur in these systems.
In other words ions on the lattice should separate into two parts with different concen~
trations of the black ions in each part. Parts with different concentrations of the black
ions ShOUId also have different values of the background charge. But it is assumed

in our InOde] that the value of the background charge is the same everywhere. Thus

the . . . . . .
background in our Simulations does not allow the systems to split into parts With

17

 

different concentrations and thus does not allow the phase separation to be observed
directly.

When we perform simulations at a. constant value of average concentraticm i. it
is possible that we may choose some. particular average value i: that can not occur
in a homogeneous systei‘ii. The energy curve E (5:) obtained in this case does not
really give the dependence of energy on concentration, but rather shows when phase
separation should occur. This will be demonstrated explicitly in the next section
using the linear chain as an example.

In simulations with a fixed value of the chemical potential the phase transitions
are more pronounced. At high temperatures (T/ J) _>_ 1 changes in the chemical
potential lead to the smooth changes in average concentration at. However, at low
temperatures (7/ J ) << 1 there are discontinuities in the 2(a) curve. We assume that
the borders of the discontinuity region are the borders of the phase separation region.

The situation at intermediate temperatures is more complicated. As an example,
F ig.2.2 shows the concentration as a function of the number of MC sweeps for a
30 X 30 triangular lattice for three values of ii at (I/ J) =2 0.08. The sharp jump in
concentration from a‘: 2 0.88 to i? = 1 that occurs at ii/ J = 0.606 shows that both

concentrations :5 r: 0.88 and it = 1 lead to the same minimum value of Helmholtz
free energy and thus are stable. Homogeneous equilibrium configurations in the range
0f concentrations 0.88 < i: < 1 have higher values of the Helmholtz free energy and
thus are unstable with respect to phase separation into two parts with concentrations
5 =0.88 audit-<1.

Fig-2.3 shows the histogram of the distribution of concentrations corresponding
to those in Fig.2.2. Peak positions give us the values of average concentrations. For

some Particular values of ,u. the system migrates between two significantly different
Concentraliions, as for example for [,l. = 0.606. We assume in this case that there

is h ~ . . . . .
p ase Separation and relative areas under two peaks give us the relative Sizes of

18

 

 

IIIIII IIITYYTIIYIITIYIIIYIIIIIIIIYIIIIIIT

(gm/J) = 0.315

 

6O

50 (TU) = 0.08

40

30 ., (ii/J) = 0.375

 

 

 

20
01/1) = 0.606 .

 

10

Distribution of Concentrations

l[III{11111111llYlIIImlIIIITIII

 

llllilllllllllillllillll[11111111

 

 

_llllilll 111111 .4144 it

0.8 _ 0.9 1
Concentration x

 

Figure 2.3: Histograms that show distributions of concentrations for the curves from
F ig.2.2. Data were collected over 10000 MC sweeps. The sharp peak corresponding
to (,u/J) = 0.315 has height 4000. The peak corresponding to (ii/J) = 0.606 has a
height 1700.

the two phases. Thus the appearance or disappearance of a peak tells us about the
appearance or disappearance of a phase. We will use the positions of the peaks when
they appear or disappear as the borders of the two phase coexistence region.

At low temperatures (T/ J) << 1 our MC procedure becomes less effective and the
System can become frozen in some configurations. One of the reasons for this is the
local character of the moves A and B that. we use to search for a new configurations
e.g, every MC move involves only one or two sites. To study the properties of the
System at very low temperatures and critical behavior of the model (when simulations
Of the large systems are required) it might be useful to implement other simulation

methods [27], but that is not the objective of the present work.

19

 

2.4 Linear Chain

There are two exact analytical methods to study the model in the one dimensional
Case of the linear chain. In particular, if the Coulomb interaction is truncated at some
distance then a transfer 7n.at7*i:r[25] technique can be applied to calculate the free
energy as a function of chemical potential and temperature. Then the dependence
1v( ’11., T) can be studied and the phase diagram can be obtained.

In case of the infinite range Coulomb interaction at T = 0, a devil ’s 501171105417]
forrrialism can be used to predict the equilibrium structure for any concentration and
calculate the ground state energy of the system. In both cases exact results will be
compared with the results of simulations to establish the precision of the numerical

methods.

2.4.1 Nearest Neighbors Interaction

If we restrict the range of interaction to nearest neighbors only, the Hamiltonian

becomes

H = J ZUQ + q)(’n,+1 + q) — #2711. (2.11)

Initially we assume that q is a constant and is not connected with concentration.

Then the grand partition function ZN for a cyclic chain of N sites can be expressed
In terms of the largest eigenvalue of the 2 x 2 transfer matrix as Zn; 2 A91“. with

8—13qu .
AW = ' 2 [1+ 7 + \/(i — 1)2 + 4765-1 , (2.12)

 

 

W'here ’7 = efifwwq‘” and [3 = 1/(kBT). The grand potential per site is given by

Q : ‘kBTln[/\mal.]. Using the fact that it 2 (n1) 2 ail/an, we can find it as a

 

A!“ a" ' .

function of the independent variables T and a. The resulting expression can then be

inve . . . _ . . . _
“Ed 11Sing the charge neutrality requirement, q = —.r, to obtain p, in terms of :1:

20

and T with the following result:

 

2 i —1 2 W2
)1 = J — 2.1.? + ,— sinh‘l (f / )6 (2.13)

3 .T.‘(l — (T)

Fig.2.4 shows the transfer matrix predictions for the u versus :7: curves for various
values of T/ J and also the results of numerical simulations. The upper left frame
shows the results of the transfer matrix calculations at various temperatures.

Note that p. : 0 corresponds to .i' = 1/2. At. high temperatures, )1 is a i‘noi‘ioton-
ically increasing function of i but at low temperatures, a has regions where the slope
is negative. This behavior is therniotlyi‘iai‘riically unstable and indicates that phase
separation occurs. The upper right frame shows transfer matrix predictions and the
results of simulations at (T/ J) = 0.50. Simulation points lie on top of the exact
curve. The lower left fran‘ie shows the results at. temperature (T/ J ) = 0.20 which is
just above the maximum temperature for which phase separation occurs. The regions
with a low density of simulation points indicate the appearance of the regions of phase
separation that occur at lower temperatures. The concentration :2 in this region is the
average over the two peaks that occurs at the intermediate temperatures as shown
in Fig.2.2 and in Fig.2.3. The lower right frame at. (T/.]) = 0.10 clearly shows the
sharp jumps in concentration that. occur at low temperatures.

In the simulations, phase separation manifests itself as a discontinuity in the
dependence 5:01.). In contrast, a second order transition would correspond to n iii-
creasing monotonically with a discontinuity in slope.

The first order transition can also be seen in plots of the grand potential per site,
9, versus chemical potential. Loops corresponding to the unstable branches appear
at low temperatures. However, we find it more convenient to plot the Helmholtz free

energy Der site, F = Q + in, or the internal energy per site, E, as a function of the

Concentration it.

21

 

F‘

0.2 0.4 0.6 0.8 l 0 0.2 0.4 0.6 0-8.

 

#—

 

    
   
  

 

] J'I‘I'U‘YIIIIIYUIUUIIUIIII'I'IIIYI IIIII IIIY'I YY‘TIIT'II'ITYTTT'ITTTYT‘I'I[IFIYYTIIII {VIII
_
h b
F p
b u b 4
F .1 F q
- d b .1
h h d
- 100 . ~
05- - - T - -05
n I- y- -- o q n
F - P d
. 050 - . .
b ' l- -1
P h d
- 020 - - <
n u q p u
h -l I- .1
I- i- d
0- - - _ s0
4 n d
I- - P- 4
- 005 « - -
b ‘ q p -
d F -
1 -. .- -
I 4 h “
F d h —1
d h u
-05 — - --05
C .1 h d 0
d b -
n F q
q ~ q
n b d
-1 - -1
j b d
p d
1 111111 [1111 lllllllllllllllllllllllllllllLLlLf 1111111 1111111111111lllllllllllllllllllllLJlll 1

 

 

 

 

 

 

L

I IIIIIII I IIIIIII FY11VIII‘IIYIIYYIITITWITYIYIIITL
q

.1

Chemical Potential (ii/J)

  

 

 

_
p
—
—
p
—
_
-
05 .—. _ .......................... 05
. '- — o u
p
- u
- —
'- —1
- q
— —I
p —
- q
0 — _0
'- —1
p m
p d
"’ —4
‘ d
h u
p d
-
q
-05 .......................... __05
. _ C
q
u
q
- u
" -1
.' h. _
p u
1 lllllllllllllllllll]lllllllllllllllllllllllllr l
— -

 

~

0 0.2 0.4 0.6 0.8 l 0 0.2 0.4 0.6 0.8

Concentration x

Figure 2.4: Dependence of the chemical potential ii on the average concentration :7“.
for the case of a linear chain with interaction limited to nearest neighbors only. The
plateaus seen in simulations for (T/ J) z 0.10 correspond to the regions of phase
separation.

A finite temperature phase transition is not expected in one dimension for finite
range interactions. However. the presence of the background charge effectively makes
this an infinite range problem and produces a first order phase transition at a finite
te{Tlperature This behavior is similar to that in the Van der Waals theory of liquids
where long range attractive interactions lead to condensation phenomena.

Figure 2.5 shows both energy E and Helmholtz free energy F = E — TS as a
furlction of 5: at (T/ J ) = 0.05. The difference between these two curves is due to the
entro13y of the system. Since the difference between these two curves is significant,
the entropy plays an important role even at quite low temperatures.

Since :7: is the independent coordinate, we apply the double tangent rule to the

22

 

 

 

 

 

0 ,_~_ I I I I I I I I l I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 4

Q -005 :— “I
\ I 2
LL} : :
V — ~
m : 1
a) -0.1 _— _:
H ,_ _.
'0-1 ,.. _
VJ : :
$— 2 I
Q) : I
C“ -0.15 r 1
>5 : :
on - . , .
H ’- - ., _
3 5 T/ 0 05 5
LL] -02 :— ( J) — - Energy -3

E " ------- Free Energy 3

- -1

3 :

-0.25—IlllIlLlll1LlllllllilllJlLlllilJiLLLlJlllllJlllll
0 0.2 0.4 0.6 0.8 1

Concentration E

Figure 2.5: Energy and Helmholtz free energy as functions of the average concen-
tration :1“: at temperature (T/ J) = 0.05 for a linear chain with interaction limited to
nearest neighbors only.

Helmholtz free energy F (or to the energy curve at T = 0) to determine the equi—
librium concentration. The slope of the tangent line gives the value of the chemical
potential H- This predicts that at T = 0 the system will separate into two phases with
c011centrations 0 and 1 / 2 if .i‘ is between those two concentrations. If 1/2 < j: < 1
then the system will separate into two phases with concentrations 1 / 2 and 1. Hence a.
plot of n versus it should display fiat horizontal sections at values of a corresponding
to t’he slopes of the double tangent lines. This process is equivalent to a Maxwell con-
str 1ICtion applied to the regions in Fig.2.4 where the ,u.(;1:) dependence has a negative
Slope.
I f in MC simulations we are trying to produce the E (‘1’) curve and we fix the con-
cerltration at a particular value that cannot exist homogeneously across the system,
then we create an internal stress in the system: the ions tend to separate into two

phases but the background charge, through Coulomb interaction, does not allow this

23

0.2

 

 

 

 

 

I I I I I I I I I I I I I I I I I I I I I I I I I I I I II I I I I I I I I I I I I I I TfiiT r
l- ’ \ u : —.
P- ”, \ \ 'A A! —-4
" I” ‘\\ T
’ N
u- II \‘ o -1
‘ \
I
0.15 — . x
’ \
Q _ ” \\ I). A O 4
\ , \
R '- l' ‘\ -4
v II \\ n .
._ ' ‘ .
a x 1.
a ,1 ‘1
a 0 l I, “ I. o
I \
g ; x 0 90 Sites
\ .
D- x '- A 240 Sites 0
5 ~ ' ‘1 . .
a) x 0 480 Sites
H I‘- : \\‘ . {D
OIOS _:' ‘\ _—
-- ' ‘I u €>
' \
I \
—' \ ..
I 1
I \
r-J \ -4
i \
J X _
1
or I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I J I I I I I I I I I I I I I

Concentration Tc

Figure 2.6: Phase diagram for a linear chain with nearest neighbor interactions. The
solid lines indicate the borders of the phase separation regions. The dashed curve
shown for :r < 1 / 2 only is the spinodal curve. Both curves were obtained with the

transfer matrix method. The points on the right. show the results of the simulations
performed on lattices with different sample sizes.

phase separation since it is constrained to be uniform. This stress should effectively
increase the value of the energy E (if) in the simulations. It. is possible sometimes
to see in structures obtained from simulations the tendency to phase separation. In
Fig. 2.7 one can clearly see this tendency for the phase separation: The average con-
centration in the top row is 1 / 2 while the average concentration in the bottom row
is 1/3.
Fig.2.6 shows the phase diagram in the temperature—concentration plane. The
SOlid curve is the transition temperature TC. This curve was obtained using the equal
area rule applied to the ,u versus if curves (see F ig.2.4) obtained from Eq. (2.13).
The dashed curve is the spinodal which corresponds to the locus of points for which

fag/<31?) = 0. The spinodal is only shown in the region 27 < 1 / 2 but it is syii‘in’ietric

24

with respect to :r = 1 / 2. The regions between solid and dashed curves correspond to
the regions of metastability.

In order to obtain the phase diagram in MC simulations at low temperatures

(T/ J ) S 0.12 we refer back to Fig.2.4. The plateaus in the Mr) dependence give
us the concentrations for which phase separation occurs. In other words, any c011-
centration on any plateau should split into the two concentrations on the edges of
the same plateau. Thus the edges of the plateaus give us the borders of the two
phase coexistence regions and the phase diagram of the system. This approach gives
good agreement. with exact results at low temperatures but. fails at. ii‘iterinediate tem—
peratures when fluctuations in the concentration are big or when system migrates
between two different concentrations. To obtain the phase diagram at» intermediate
terllperatures 0.14 g (T/ J) S 0.20 we used the observations and interpretation that
was discussed above in Fig.2.‘2 and Fig.2.3. In Fig.2.6 data for the sample with 90
sites were obtained from the borders of the jumps in .r(,u.) curves shown in Fig.2.4.
The remaining points are obtained from the distribution of concentrations at different.
T and p using larger sample sizes. The peaks in the distributions of concentrations
for the sample with 480 sites are narrower than those with 240 sites, as expected.
Tile full width at half of the maximum for the peaks at temperatures (T/ J) 2 0.14
is about 0.05 and at temperature (I/J) :3 0.18 it is about 0.1. Thus at temperature

(T/ J) 2:: 0.18 peaks corresponding to two different concentrations overlap.

The good agreement between the exact phase diagram and the one obtained from

Slrl'lulations indicates that our interpretation of the histogram of the concentration

frequencies is correct. It also indicates the precision of our technique.

2'4- 2 Finite Range of Interaction

The Hamiltonian of the model at fixed concentration .1: is given by Eq. (2.3).

In . . . . . .
tlle case of nearest neighbor interaction J1 and a second neighbor interaction J2

25

OOOOOOOOOOCOOOOOOOOO
OOOOOOOOOOOOOOOOOOOO
OOOOOOOOOOOOOOOOOOOO
OOOOOOOOOOOOOOOOOOOO
OOOOOOOOOOOOOOOOOOO.

OOOOOOOOOOOOOOOOOOOO

Figure 2.7: The linear chain sample of 120 sites with interactions up to fifth neighbors
included at (E = 3/ 8 = 0.375 and temperature (T / J ) = 0.01. The second row from
the top is the continuation of the first row and so on. The picture illustrates the
tendency of the system to phase separation: in the top row .1: z 1 / 2, while in the
bottom row a: a: 1 / 3. Separation cannot really occur due to the same backgrormd

charge at every site of the chain that corresponds to .7": = 3 / 8.

only, the ground state energy can be obtained directly using the following reasoning:
for concentrations in the range 0 < :17 < 1 / 3, the background charges (white charges)
corltribute an amount —(.]1 + J2);1'2 to the energy per site from the last term in
Eq- (2.4). The black charges can be placed on every third site so as to avoid the
reI)ulsion in the first term. However. in the concentration range 1/3 < :r < 1/2
the repulsive interactions contribute an additional amount (J1 + J2)(;r — 1 / 3) to the
energy. The energy in the range 1 / 2 < 1' < 1 is obtained using the symmetry property
EC 1') = E (1 ~- 1'). Similar reasoning can be used for larger ranges of the interactions
bu t the expressions become more complicated.

Fig.2.8 shows the ground state energies as a function of concentration if at T = 0,

for finite ranged interactions .1" = J/n, where n. = 1,2,3.4,5 and 6. The solid

curves were obtained using the arguments described above and were verified using
the numerical transfer matrix method results at very low temperatures. The MC

si . . . . . .
qulations were performed on lattice samples wrth 90 and 180 Sites and periodlc

b . .
0L1 Ildary conditions.

26

 

A\.\.VJ\ v

0.. 0. l. 0.
0. . 0. .

aJu .- mJ. -uaJnn \A.”N-.aJ-m.~

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Concentration x

O 0.25 005 0.75 0025 0.5 0.75 1
0 IttTI‘IY'IIII’YTIYYIIIIYIIIVIIIIIIUYIII IIIIIIIITITITIITITTIIIrrTTIYIITIII[11"
I- - -
f- —1 - -
y— -. -— —
h d F- _1
-0 05 — — — — -0 05
' — —1 r— 1 '
— —1 — -1
b —1 -— —1
i- -1 r- —<
-01~ — — — -01
° -- -1 l- —4 '
>- -4 u— .1
r— —1 v— —-.
— u +— _.
-0 15 ~ - — — -0 15
. >- d r— —4 '
P. -l ’- -‘
-4 >— —
_ h —1
o ’- _‘ __ __ o
L— —4 i— —
1— t _
r- —l
‘ ‘ LLLllllllllllllllllvlLlllLll11111111111 111111111ll1111111lililelllllllllllll
\ O ‘I'IIIIUII'ITIIII'IIIII'IIIIIIIYTrYII III'IUIIIIIITY1TTITIITIIITrIIIIIIIIIIYI 0
LL) — — - ..
_ -—1 h -1
v I— —4 — —1
r—- — — —1
-0 05 — a — — -0 05
Q.) - e ~ ~ - -
H —l r- -1
O“ K —1 r— —‘
i— .4 P- —I
' r— — — —1 ’ ‘
H — — .— -1
D _ — F d
r— — l— A
9-. 0 15 -— — — - 015
- -
o _ _‘ _ _ o
P“ —4 I—
>‘. ._ .1 _ 1
Q0 “ ‘ ' ‘
H ' p— —< - - '
0.) ~ 1 — —
h- D— -1
a _ _
m luJLlLllLLlllllllllll11111]lllllllll lllllll‘lllllllllllllllllllllllllllllLl
O I'ITIIIUI'UIIIIIII'IIIIIIIIIIIIIIYIIIII WITIIIIII'ITITTTTTITYYIIIIIIIIIIIIIIY O
l— —1 I— —1
_ - _ d
i— -4 I— -1
n— —1 r— -
’ - -l r- - '
_ 4 h i
n— —1 — —4
n- — .— -1
' >— _‘ p— _‘ t
P' —1 >— 1
>- -4 r— —1
i— -1 i— d
' v- -4 >— -. 0
I- -1 v— 1
h- —4 ~— _‘
u— —1 - q
- -
. h- -: r— —1 °
- c1 — d
_ —1 r— .4
I- cu — —1
llJuLllllllllllllllllllllLllllllll‘lll lllllIJJLIJIIIllllllllllllllllllllllL
O 0.25 005 0.75 0.25 005 0075 1

FigLIre 2.8: Dependence of the energy per site on concentration, at very low tempera—
tureS, for the linear chain with interaction extending up to first, second, third, fourth,
fifth and sixth neighbors. Points from simulations are plotted on top of the exact solid
_Curves obtained with transfer matrix technique. Dotted curve shows limiting case of
lnfirlite range coulomb interaction.

27

The panel 2 in Fig.2.8 shows the energy curve for the interaction between nearest
and second nearest neighbors only. It is easy see that for .i‘ : 1/2 it is energetically
favorable for the system to split into two parts with different average concentrations
f1 = 1/3 and i2 = 2/3.

If we apply the double tangent. construction to the ELI?) curves in Fig.2.8, cor-
responding to the interactions with larger range, we find additional transitions com-
P'dI‘Ed to the nearest neighbor case. For interaction up to second neighbors we have
first. order transitions from :r = 0 to :r = 1/3, I. = 1/3 to .r = 2/3 and :r = 2/3 to
17 = 1. For neighbors up to the third included, there are four transitions involving
513 =— O,1/4,1/2,3/4 and 1.

As the range of interaction increases there is an increasing amount of structure
in the energy curves but the double tangent construction does not always lead to an
increasing number of transitions. It will be shown in the next section that in the limit
Of Coulomb interactions, the energy curve approaches a form which predicts only four
transitions involving the values .1? = 0, 1/3, 1/2, 2/3, 1.

Phase separations predicted from the energy curves cannot occur in our MC sim-
UIations because this phase separation requires two different. values of the background
Charge, while in our simulations background charge is constrained to have the same
Value at every site.

However if concentrations of the two parts are close to each other then the ten-
CleIlcy to the phase separation can be seen even if the value of the background charge
is the same everywhere as shown in Fig.2]. This behavior can also be seen on the

fifth panel of F ig.2.8 that shows a disagreement between exact curve and results of
Si ITlulations for the case of i1’1teraction up to fifth neighbors. There is a region between

0‘ 325 g 5: 3 0.425 in which simulation points seem to follow a horizontal line while
thee , .fi r . . ‘ .fi
xact curve has a peak in this region. We believe that the Oilgm for this behav

1S the tendency for the phase separation descrlbed above. Thus every Simulated

28

-w-..F. ' mil-muff
U

 

 

 

 

 

0.00 I I I I I T I l I T Y I I T T I Y I 1
.1
t -
g -0.05 — —
V " -1
D >- ..
.‘.-2."
U) " "
8 "' -4
o. -O.10 — ~
>5 h -
an __ _
Q)
C " .
LL} _ _
-0.15 — —
' l .l L l l l 1 i l l L l l l l I l I l
0 0.2 0.4 0.6 0.8 1

Concentration x

Figure 2.9: Energy per site as a function of concentration it for a linear chain with
Conlomb interactions. The solid curve was obtained by knowing equilibrium struc-
tlll‘es predicted with the devil’s staircase formalism on the lattice sample with 1000
Sites. Circles give the results of the MC simulations at low temperature for samples
With 180 sites.

DOint in this region is a weighted average of two concentrations that correspond to
the edges of this horizontal region. The weight is proportional to the fraction of the
VVhole sample that is at a given edge concentration.
In order to apply the transfer matrix method to the case with range of interaction
T" 2 2 neighbors, we can group the sites into consecutive blocks of length n and the
illteractions are then only between nearest blocks. The transfer matrix formalism can
136‘ Used again, but the Hamiltonian given by Eq. (2.6) leads to the transfer matrix of
SiZe 2" x 2" and it is difficult. or impossible to solve the problem analytically. However,
i t is possible to calculate the largest. eigenvalue of the transfer matrix numerically and
O btain the thermodynamic properties. we have used this method to study the energy,
I\Ielfhholtz free energy and chemical potential as a function of concentration if. Our

1‘ e . . . . .
8111 ts are in agreement w1th the predictions made from the energy curves, shown in

P-
lg-2.8.

29

 

'fi 33.: _
4
I

 

0.8

IIIIIIYIVITITIIIIIIIIIIIIIIIITTIIT‘FITIIIYIIIIIIIIIW_,
I
I .
I

.I
I

 

 

 

 

 

 

 

 

 

 

 

: 60 sites, (771) = 0.40 j

; -.- 60 sites, (T/J) = 0.20 ;

: 60 sites, (T/J) = 0.15 3

E 3 60 sites, (T/J) = 0.10 , ’ _j
3:3 0.6 _— . 180 sites, (T/J) = 0.10 ,z ’ ‘5‘
_, : — Zero Temperature Limit , r’ " ,
9'8 C /’ i :
Q, 0-4 T ‘33:; 1
'5 E 'Ir' “la/530.5 ‘ i . . . 1 L . . . i . . . . l . . . . E
8 '- Ir! 7; ‘ '..,. . 1 4

. .- - I : f—H—t a
E _ 1,. .,// 0.25 1 m f .

a) : : -. ~- (T/J) = 0.015 t 2

'5" 0.2 —- .2“ 0 4 .. . - f7
U : ’11,» j 960 Sites i 1
i -025 —§ _ , . E I

: .JV -05 v . . . I v v , . I . v v . I . . . . j

_ d" 0 0.25 0.5 0.75 1 1

O _l i l 1 J J l l 1 l 1 l 1 l 1 1 i 1 l l l l l l l l l l l l l I .l l L L4 4 L4 I l l l 1 l 1 i L 1 H

0.5 0.6 0.7 0.8 0.9 1

Concentration 17

Figure 2.10: Dependence of the chemical potential ,u on concentration .'ir for the linear
chain with Coulomb interaction in the region :1’: > 1 / 2. The solid line shows the limit
of zero temperature obtained from the energy curve on F ig.2.9. The inset shows the
results from the lattice sample with 960 sites at very low temperatures in the whole
range 0 < 515 < 1.

2.4.3 Infinite Range Coulomb Interaction

If the Jij correspond to the bare Coulomb interaction, the interactions satisfy
the positivity and convexity condition [17] which allows the ground state to be found
for any rational value of i. If .7: = 1/n where n. is an integer, then the black charges
are equally spaced along the chain at a distance of n neighbors apart and form the
One dimensional analogue of the Wigner lattice[28]. In the more general case where
57 2 19/ q is the ratio of two integers, the ground state configuration is periodic with
period q and has p black charges in each cell. If there is a black charge at site 0 then
blaCk Charges are at the sites with numbers [nq/p] = [vi/2:] where n is any integer and

[A] denotes the integer part of A. Since the structure is known for any rational value

30

 

      

 

 

 

 

 

0.12 I I I I Ij I I I I I I ITI I I I I I—I I rI Ij I I I I I I I I I I I I I I I I I I I I I I I I I _

0.1 g— “2

S a 2

E 0.08 r -

v : 3
H L
a 0.06 :
H ..

O.) I I

Q. I I.

E —_ 4

Q 0 045 I

I" E a»

0.02 E— a

0 :v I 1 1 l 1 1 J 1 1 i 1 l 1 1 l l' 1 1 l 1 L 1 1 v 1 1 1 1 I 1; 1 1 1 1 1 l l I l 1 1 1 l l l I 1:,

0.5 0.6 0.7 0.8 0.9 1

Concentration x

Figure 2.11: Phase diagram for the linear chain with Coulomb interactions in the
region 1: > 1/2 obtained from MC simulations.

of :1? the ground state energy can be calculated using the same techniques that. was
used to calculate energies in MC simulations. But since in this case it is not necessary
to run a relaxation procedure much bigger samples can be considered.

Fig.2.9 shows the ground state energy as a function of .r. for the Coulomb po-
tential obtained using both the exact ground state configurations and our simulation
technique. The value of the energy at :1: = 1/2 corresponds to alternating black and
white Charges and is equal to —(ln 2) /4. For smaller values of .1: = 1 / n corresponding
to Period 72. Wigner lattices of equally spaced black charges, the energy per site is
given by E = 132111 .1:. For values of :1: between these values the energy is slightly
larger than if we use the same formula for all 1:. The energy curve has a sequence of
cuSPS located at all rational values of :1: (devil’s staircase).

Since at T = 0, the chemical potential is given by n = 813/ 0.17, we predict that

at zero temperature In. versus .i‘ will display a series of jumps as shown in Fig.2.10.

31

The double tangent rule applied to the energy curve predicts phase separations for
some values of 5:. For 0 < I < 1 / 3, the system should separate into phases with to
r = 0 and a: = 1/3 whereas, for 1/3 < .1; < 1/2, it should separate into phases with
:r 21/3 and .17 21/2.

Phase separation is due to the second term in Eq. (2.4) which is entirely due
to the background. It corresponds to a long range attractive interaction. In the
absence of the background, the model would not display phase separation and the
chemical potential versus concentration curve would be a devil‘s staircase with an
infinite number of jumps corresponding to the rational values of .17.

The Devil‘s staircase formalism can be used to obtain the E(1) curve at zero
temperature and thus predict #(l') dependence at zero temperature. In order to
obtain the p.(.r) dependence at non zero temperature we used MC simulations. The

p.(.“lf) curves at several temperatures obtained from MC simulations are shown in
Fig.2.10. Simulations were performed on samples with 60, 90, 120, 180 and 960 sites.
The role of size effects can be seen for the curve (T/ J) = 0.10 for which the results
from the samples with 60 and 180 sites are shown. For higher temperatures the
differences are less significant. \Ve found that up to very low temperatures there are.
almost no differences between the results on samples with 120 and 180 sites.

The inset in Fig.2.10 shows a rather large difference between the predictions
for 11(1) from the devil’s staircase forn’ialism and from the results of simulations at
low temperatures. There can be several reasons for this disagreeii’ient. One reason,
discussed earlier, is that simulations at low temperatures may not be reliable since
the SYStem can become frozen in some local minima that it cannot leave due to the
large transition barrier associated with the local character of moves .4 and B. On the
other hand the entropy contribution can make the “cusp” at a? = 1 / 3 and :17 = 2/ 3
in the E (:1?) dependence more “rom'ided” and this can lead to an extended range of

chemical potential near concentrations :r : 1 / 3 and :r = 2 / 3 on the ,u(.z:) curve.

32

In order to obtain the phase (:liagrain shown in Fig.2.11 we used the technique
already discussed above. As expected the size of the plateaus corresponding to the
regions of phase separation decreases as the temperature increases and the concen-
trations of both phases becomes the same at the critical tell'll‘)€1‘atlll'€. Only half of
the phase diagram is shown since it is symmetric with respect to r = 1 / 2. At tem—
peratures 0.04 g (T/J) S 0.08 the data were obtained using a sample with 480 sites.
At other temperatures the data were obtained using a sample with 960 sites. There
are three distinct regions of phase separation that correspond to the following tran-

si tions in concentration at zero temperature: 1/2 —> 2/ 3 with critical ten‘iperature
(T/J) 2 0.03, 2/3 —+ 3/4 with critical temperature (T/J) 2 0.0175 and 3/4 —+ 1

V’Vi th critical temperature (I/]) 2 0.11.

2 . 5 Triangular Lattice

The Hamiltonian for triangular lattice with Coulomb interaction is given by Eq.
( 2 - 6) :
H = Z J,J(71,+ (1)01] + q) -— In 271,. (2.14)

<ij>

Tile details of model and simulation techniques were described in section 2.2 and

2‘ 3 ~ In the 2D case with long range interaction we used MC simulations as the main

I’llet hod to study the system. Different sample sizes were used in the previous studyIl]
in C)I‘cler to obtain E(i7) dependence shown in Fig.2.13.

For the triangular lattice there are multiple “cusps” in the E (1?) curve that are
Sirllilar to the devil’s staircase behavior in case of the one dimensional chain. In
tile region 1/2 _<_ if S 1 significant cusps in E (.1:) curve occur at concentrations:
1/2 , 2/3, 3/5, 3/4, 6/7.

Ground state configurations for .r : 1/2, 3/5, 2/3 and 3/4 are shown in Fig.2.12.

I Q _
I l ('1 iStributions obtained on the 30 x 30 sample with MC simulations at three different

33

 

 

 

 

0_ I I j H I I l I l ..

c g 2
g -005 _— €
3 2 3
E71 : 1
8. -0-1 r 1
>5 C 1
OD - -
H " 4
0.) " _
I: I Z
”'l -015:— €
: l l l l l l l l l :

0 1/7 1/4 1/3 2/5 1/2 3/5 2/3 3/4 6/7 1

Concentration Tc

Figure 2.12: Energy per site as a function of concentration :1 for the triangular lattice
with Coulomb interactions. The solid curve connects the points obtained by simu-
lations in previous workIl]. Dashed vertical lines give the position of concentratioi‘is
that are stable in MC sin'iulations at very low temperatures.

val ues of In. at (T/J) = 0.01 are. shown in Fig.2.l4. In the ground state at. :r = 6/7
‘Vllite sites should form a triangular W’igner lattice with the spacing (ix/7 between
i 0118 I 1],

For values of .r = 1/(111‘2 + r12 + 11111.) S 1 / 3 or 1 — .1: 2 2/ 3 where m and n. are
illtengers, the ground state configurations are triangular Wigner crystals[28, 29] and
Elle ground state energy can be calculated exactly using the numerical formalisn'i of
B 01 IS all and Maradudii'iIBU]. For these triangular Wigner crystal structures the energy

is giVen exactly by E(.1:) = —2.10671262(;r3/2 — 11:2). Between :1: = 1/2 and .1: = 2/3,
trle ground state configurations have a rectangular rather than triangular geometry.
Reg ular rectangular structures can be formed at concentrations .1 = l/ (2(1 + 1)), where
Z -

1‘ ~ . .
3’ any integer[1]. Thus I = 4 corresponds to :r = 2 / 5 and is the same structure as for

’—

§

3/ 5 with the interchange of black and white sites. The same numerial method[30]

Q - -,
‘11 1 be used again to obtain E(3/5) = —0.1709803 and E(1/2) = —0.1755589512

34

00.0.0000...
000000000000
0.0.0.000...
000000000000
0.00.0000...
000000000000
00000000....
000000000000
0.00.0000...
000000000000
00.000.00.00
000000000000

000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000

000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000

000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000
000000000000

Figure 2.13: Ground state configurations for the triangular lattice with Coulomb
interaction for :1: = 1/2, 3/5, 2/3, and 3/4.

The first row of the table 2.1 shows the values of the energies obtained from
sirnulations at (T/ J ) = 0.002 with sample sizes 2 18 X 18. Per every concentration
the sample size was chosen to be commensurate with particular concentration. The

Second row shows the results of simulations from the lattice samples of size a: 30 x 30.
The third row of the table shows the exact values of the energies obtained with the
nletllod of Bonsall and Maradudin.

The major difference between the simulated and the exact values arises because

Our Simulation method can not relax the system to the exact ground state configu-
rations for concentrations 1/2, 3/5, 3/4, 6/7. On the other hand it is easy to obtain
t lle ground state configuration for 1: = 2 / 3 and the agreement between the simulated
an (1 the exact values of the energy is much better in this case.

Double tangent construction applied to the energy curve Fig.2.13 shows that

011 13’ concentrations 0, 1/4, 1/3, 1/2, 2/3, 3/4, 1 are stable with respect. to the phase

Se
IDaJI‘ation at zero temperature. Thus at zero temperature the system is always a

35

0000000000000000 0000000000000000 0000000000000000
00000000000000. 00000000000000. 00000000000000.
0000000000000000 0000000000000000 0000000000000000
00000000000000. 00000000000000. 00000000000000.
0000000000000000 0000000000000000 000.000.000.000.
00000000000000. 00000000000000. 00000000000000.

0 0 0 0 0 0 0..
0.00000000000000 0000000000000000 0000000000000000
00000000000000 00000000000000. 00000000000000.

Figure 2.14: Lattice structures obtained from simulations at three different vahies
of p at (T/J) = 0.01. Concentrations of black ions are close to 1? = 1/2 = 0.5, :r =
3/5 = 0.6, :1: = 6/7 9: 0.857 for the pictures from left to right.

mixture of parts with these concentrations. For example if :7: is between 0 and 1/4
tllen the system should split into two homogeneous parts with .’L‘ = 0 in one part and
1" = 1/4 in another.

The slope of the double tangent lines gives the values of chemical potential at
“'Ilich the transition from one concentration to another should occur. In the region
‘1‘ > 1/2 we have: (p./J)1/2_.2/3 = 0.0003. (p./J)2/3_s3/4 = 0.5087, (,u/J)3/4_.1 =
O ' 5 1 59. The corresponding dependence 11(1) at zero temperature shown in Fig.2.15

as the dotted curve. Direct simulations of the ,u.(:r) curve were performed on lattices

 

f 1/2 3/5 2/3 3/4 6/7
E3.m(18x18) 0.17066 0.16534 0.16953 0.12450 0.06794
Esim(30X30) 0.16964 0.16456 0.17070 -0.12671 -0.06807

Ema 0.17556 0.17098 0.17136 0.13167 0.07076

 

 

 

 

 

 

 

 

 

 

 

T .

Q a'l’ble 2.1: Comparison of energies at special concentrations obtained in simulations

at lattices with sample sizes 18 x 18, in the first row, ~ 30 x 30 in the second row
1 d exact values calculated with method of Bonsall and Maradudin in the third row.

36

 

 

  

 

  

 

 

 

E I I I I I I I I ll,’

0.6:— '

E

C: 0 4 E— E
33 02;— é
a a s
8 E 3
o :— a
a. OE ‘ ., 3
Te _..;.;,>' 9x9, (T/J)=0.20
-§-0-2 :— r,ng 9x9,(T/J)=0.10 g
a) E ,7; 9x9, (T/J)=0.08
"50.4 g— ' ; + 18x18, (T/J)=O.10 —;
o—o 30x30,(T/J)=0.03
-0.6<E, .71 : , - Zero Temperature Limit _f

E” /'I I I I I I I I I

0 1/6 1/4 1/3 2/5 1/2 3/5 2/3 3/4 5/6 1

Concentration )_c

Figure 2.15: Dependence of chemical potential ,u. on concentration j? for the triangular
lattice with Coulomb interaction. For the cases of temperatures (T/ J ) = 0.2 and
( T/ J) = 0.1 the size of the sample was 9 x 9 and no points are shown since their
density is high.

Of d ifferent sample sizes: 9 x 9, 18 x 18, 18 x 20, 24 x 24 and 30 x 30. Some results
are Shown in Fig.2.15. Crosses for the sample size 18 x 18 show that there is no big
di fference with the results obtained on the sample of size 9 x 9 at (T/ J ) = 0.10. For
hi gher temperatures the difference is even less significant.

When we consider long range interactions in 2D, the amount of time needed to
Calculate the energies before and after the MC step increases quickly with the size of
t he system. In order to save time we tried to perform calculations on smaller samples
\v 11611 possible. Thus at high temperatures, we performed simulations on small 9 x 9

SE3"11'1ples. Every point is the result of averaging over 20,000 MC sweeps. Even at

C T/ J) = 0.05 there is almost no difference between the results of simulations on the

37

lattices 18 x 18 and 30 x 30. However, at very low temperatures, there are some
structures that can be observed only on large samples. Fer example at a temperature
(T/ J) g 0.03 we can clearly see concentrations 3 / 5 and 6 / 7 on the 30 x 30 sample. but
we cannot see these concentrations on the 18 x 18 sample and they are not pronounced
in case of the 24 X 24 sample.
At low temperatures (I 2 0.01) there. is a significant difference between the
direct simulations of 11(17) dependence and the predictions from the E(.1>) curve. In the
simulations, the ranges of stability (with respect to the change in chemical potential)

at concentrations :1: = 1 /4, 3 / 4 and 1/ 2 are much larger then they should be according

to predictions from the energy curve, while at concentrations .1: = 1 / 3 and 2/ 3 the

ranges of stability are smaller. In simulations we also see stable concentrations 2 / 5

and 3/ 5 that should not appear according to the energy curve. Several effects can

lead to this disagreement.
First of all, the difference between the energy and the free energy curves can be

Significant even at very low temperatures. This difference can lead to higher stabilities
0f some concentrations. l\sIoreover some concentrations that should not appear, as

follows from energy curve, can appear due to entropic contributions to the free energy.
— 2/ 3 are highly ordered and have

For example, configurations at :1: = 1/3, and 1?

r elat ively small entropy. The concentrations between .1: = 1 / 3 and .r. = 2 / 3 are highly

(iisordered (large entropy) and have approximately the same energy as energy of the

System at concentrations 1 / 3 and 2/ 3. Thus one can expect that concentrations 1 / 3

and 2/3 will have smaller range of stability in ,11, while concentration 1/2 will be
more stable than follows from energy curve predictions. This is in agreement with
Fig - 2.15. This feature can also lead to the appearance of some configurations with
1 /3 < :1: <1/2 and 1/2 < :1: < 2/3: for example :1: = 2/5 and a: = 3/5. It also explains

1y we see these concentrations on larger samples and do not see them on smaller

8 -
(11 r1IDIes.

38

Another reason can be in the way simulations were performed at a constant
value of the concentration. For example, if the system has average concentration
1 / 2 < i: < 2 / 3 then it should separate into two parts with concentrations .1' = 1/2 and
.r = 2 / 3. But since at every site the background charge is the same, the system cannot.
separate. This tendency of the system to phase separation, that cannot occur, creates
internal stress and effectively increases the energies of intermediate concentrations.
When simulations are performed at a constant value of the chemical potential there
is no problem connected with background charge that. prohibits the phase separation
and there is no stress in the system. Thus energies at intermediate concentrations
effectively become smaller.

The last reason that we can mention is the local character of the A and B
operations that were used to introduce changes in the system. Their local character
can also become important at low ten‘iperatures.

All the reasons discussed above can lead to the larger regions of stability for
concentrations 1' = 1/4, 1/2, 3/4 and to the smaller regions of stability for the con-
centrations ;r. = 1/3 and .r = 2/ 3. They also explain appearance of concentrations
1‘ = 2/5 and :1: = 3/5.

One can obtain the phase diagram for the triangular lattice in the same way
as for the linear chain. In order to obtain the phase diagram of the system at low
temperature we had to perform up to 100,000 MC sweeps on the lattice with the
88‘I'nple size equal to 30 x 30. The phase diagram of the triangular lattice is shown in
Fig- 2. 16. Only half of the phase diagram is shown since it is symmetric with respect
to 513 = 1 / 2. There are five distinct regions of phase separation that correspond to the
transitions in concentrations 1/2 —> 3/5, 3/5 —-> 2/3. 2/3 -—> 3/4. 3/4 ——> 6/7, 6/7 —+ 1
at Zero temperature. Some regions of phase separation exist only at low temperature

and it is difficult to obtain an accurate phase boundary for them.

39

 

   
 

 

0.12- I I I I I:

E a

E i

0.1: ‘2

S E 2
R. L r
a E o 0%
a : :
a 0.06:— i’ ‘ “E
8 : :
.— 41 0 WI

CL 5 5
E 0.04:— ., I .5
cu ; . . 1. 1»:
E— E U ' E
1 I

0.02 . \ —:

.. g

 

 

 

 

 

 

but

3 __ - I
(l/2 3/5 2/3 3/4 .. 6/7
Concentration x

Figure 2.16: Phase diagram for the triangular lattice with Coulomb interaction in the
region at > 1 / 2 obtained with .\IC simulations.

2.6 Conclusion

We have introduced and studied a model to describe the possible ions ordering
in layered double hydroxides. In the model ions situated at the sites of the triangular
lattice interact through long range Coulomb interaction. The exactly solvable example
of the liner chain was used get insight into the model properties and to demonstrate
the precision of the MC simulation methods employed. The model predicts multiple
phase transitions and phase separation regions.

Our results are in agreement with experimental measurements in a sense that.
concentrations :1: = 1/ 4 and a: = 1/ 3 are special [3, 4, 5, 6, 7]. However a large number
of predicted possible phases with different fractions of metal ions can lead to the
situation when effects are hard to see in experiments. If we assume that some phase
Separation really occurs, then we should allow for the fact that there are different

amounts of [CO3]2‘ anions in the different regions of the same gallery. Regions with

40

a large number of [CO;312_ anions should have a larger interlayer spacing than the
regions with a small number of [C03]2“ anions. This should lead to mechanical
stress in the system and that will also resist the tendency to phase separation. This
competition between the Coulomb tendt—aicy to phase separation and mechanical stress
can lead to the disintegration of the compounds and limit. the composition range
over which laYered hydroxides can be symhesized. However the effect of nonzero
temperature can move the system over the phase separation boundary and bring the
system to a uniform distribution of charges. In this case there should not be internal
stress in the compounds. Thus Al substituted layered Ni l’Iydroxides that are stable
at higher temperatures may become unstable and disintegrate at lower temperatures.

It would be interesting to see detailed experimental measuremeiits that can sup—
port evidence for the charge ordering and the presence or absence of phase separation

in these systems.

41

Chapter 3

Absence of Decay in the
Amplitude of Pair Distribution

Functions at Large Distances

3. 1 Introduction

Pair Distribution Function (PDF) or radial distribution function C(r) has been
used to study atomic structures of materials since 1927l32, 33]. It is a real function of
a single real variable: radius 1". Peaks in PDF centered at values of r =< f1) > that
correspond to the average distances between atomic pairs. Thus in principal PDF
contains rather limited information about the structure of the materials. Because of
it PDF is used mostly to characterize local atomic environment in materials with the
absence of long range order like amorphous materials or liquids. Importance of PDF
function is caused by the fact that PDF is related in a simple way to the scattering
intensity in X-ray or neutron diffraction experimentsl34].

If a material consists of atoms of only one type then experimental PDF Gex(T)

can be obtained from reduced scattering intensity F (q) E q[S(q) — 1] via the sin

42

Fourier transformation:

. 2 0° ,
Cer(r) = —/ F(q)s1n(qr) dq, (31)
7T 0
where S q _=. 55%, I q is scatterin ‘ intensitv avera ‘ed over the directions of momen—
1» f 8 .. 8

tum transfer vector q, f is the (q dependent) atomic form factor in x-ray scattering.
and is the (constant) scattering length in neutron scattering. The total number of
scattering centers in the sample N in practice is defined in such a. way that S ( q) goes
to unity as q goes to infinity-(34].

As it was already mentioned PDF is related to structure of a material. If material

consist of only one type of atoms then:
C(r) = 47rr[p(r) — [)0], (3.2)

where the average number of atoms in a unit volume (average number density) is
p0 and the number of particles in the spherical annulus of thickness d'r is given by
dN(r) = 47rr2p(r)dr. It is clear that radial density p(r) should oscillate around p0.
In order to compare modeled structure with real structure of a material experiment
modeled p(7') usually is calculated as follows.

It is a widely used approximation justified by the Debye-VValler theorem, that in
solids, if '1",- and F]- are equilibrium positions of atoms 1: and j, the density with respect
to atom 1'. created by atom j is given by[34]:

2

1260") = 7—3 exvl ], (3.3)
where fij = 7"]- —- 1‘; and 0,2)- =< (fij- < fij >)2 > is the mean square deviation of 171:}
from its equilibrium value < 17,]- > due to atomic vibrations. In order to obtain p(r)

it is necessary to perform angular averaging of pij(F). This can easily be done in 3d

43

leading to:

1 (r — r, )2 (r + '13.»)2 .
0110‘) = , {exp [—7240 ‘“ EXP l—TJ—li' (3'4)
4mm, / 27mg. 01]- 011'

 

Since in solids deviations of atoms from their et‘iuilibrium positions are much smaller
than interatomic distances 0,]- << 73,-, the second term in the numerator can be ignored
and 73:,- in the denominator can be substituted with r. The total radial density with

respect to the atom 17 is obtained as a sum of contributions from different atoms j:

A; 1 l (7' — r,~)2
p10") = 1—3 X —_ e3’<Pl———2—2]—l- (3-5)
‘ 7” fig 200?] 011)

Finally in order to obtain p(r) it is necessary to average p.,-,( '1‘) over all sites N in

the sample:

1 N ,
710‘) = N EMF)- (36)

PDF C(r) is obtained by substitution of p(-r) into (3.2)[34, 58, 36. 37, 38]. From
formulas (32,3536) f(_)llows that peaks in C(r) occur at values of 1' that correspond
to the average interatomic distances in a material and that since p(r) oscillates around
pO PDF should oscillate around zero.

By comparing PDF calculated for some model structure with experimental PDF
obtained from scattering data it is possible to verify how close the model structure
is to the real structure of the material. PDFs for amorphous materials and liquids
contain only a few peaks, whose amplitudes quickly decrease as r increases, as it. is
shown in panel C of Fig.3.4. Because of this, little attention was paid to the behavior
of PDF at large distances.

In the last 20 years PDF has also been used to extract information about in-
teratomic interactions in crystalline materials[39]. This can be done because 0,3 is

interaction-dependent[37, 38]. Experimental PDF for crystalline materials decays

44

 

 

 

 

 

 

 

’- I I I I I I I I I f I I I I I I I I I _,
200 _ —
A 7 Experimental F(q) from Ni at 15 K -
cr - _
E: 100 *- _,
0 l l L l L l l l l l i l J l l J l i

0 10 20 30 40

Momentum Transfer q, (131")

    
  

 

I I I I I I I I I I I I I I

Experimental G(r) from Ni at 15 K

(3(r)

 

'Zililili

a _
l I l

 

0 50 100 150
Distance r, (A)

Figure 3.1: Experimental F (q) and 0(7) of N i at 15K measured by neutron scattering.

significantly slower than for amorphous materials, as shown in Fig.3.1 that shows re-
duced scattering intensity F (q) obtained in the neutron diffraction experiment at the
NPDF diffractometer at Los Alamos National Laboratory and corresponding PDF,
G(r), obtained from F (q) by Fourier transformation (3.1).

However traditionally there still was interest in behavior of PDF at small dis-
tances because that is the region from which structural information can be easily
extracted. Thus it seems to us, to the best of our knowledge, that although PDF
was used to characterize materials for almost a 100 years, behavior of PDF at large
distances was never carefully discussed. For example, it is unclear why PDF decays
at large distances. In fact, there can be three different reasons that can lead to the
decay of PDF with an increase in distance.

The First and the most natural reason is the following. Note that from (3.5)

follows that the value of p,(r) at 'r is determined by those atoms that are in a-vicinity

45

 

Figure 3.2: Only those sites contribute to the value of 0,-(7‘) that are in a vicinity of
the distance r with respect to site i.

of the spherical shell of the radius 1'. 111 other words, in order for the site j to contribute
to the value of p,(r) the following condition should hold: |r — 'rU-I ~ 0' (see Fig.3.2).
The volume that these atoms occupy can be estimated as (fl/W. ~ 477720, while the
number of atoms in this volume can be estimated as dNW ~ 41rr2p(r)a. As r increases
both de and dVW increase and it is natural to expect that p(r) ~ d/ngr/dl/W
converges to p0 leading to the decay of peak amplitude in PDF (see (3.2)). However
the question of how {)(1‘) converges to p0 was never properly addressed and we will
see that prefactor 7' in (3.2) plays a very non-trivial role.

The second reason that can lead to the decay of PDF with increase of 'r is
averaging (see (3.6)) of p,(r) over different central atoms 1' each of which can have a
different atomic environment.

Finally the third reason that can lead to the decay of experimental PDF is finite

instrumental resolution that increasingly broadens the width of the peaks with an in—

46

crease of the distance. In order to achieve better agreement with an experiment, PDF
calculated according to (3.2,3.5,3.6) is usually convoluted with so-called instrumental
resolution function that is r dependent-(37]. This convolution increases the width of
the peaks at large values of 1' more significantly than the width of the peaks at small
values of 7‘.

It is clear that finite iI'Istrumental resolution leads to the decay of PDF, but do
both the first and second reasons also lead to the decay of PDF? Here we are trying
to discuss this issue carefully. This study was initiated by the observation (a puzzle)
that PDF calculated for the fee structure of crystalline Ni, in assumption of infinite
instrumental resolution (no convolution was made), does not decay at large distances
at all. The average amplitude of peaks in PDF function persists up to the very large
distances as shown in Fig.3.3.

All atomic positions in the fee structure of Ni are equivalent. Because of it
there is no need to perform averaging over different atomic sites (3.6) Thus it turns
out, surprisingly, that increase in the volume of the spherical annulus with r does
not lead to the decay of PDF (first reason mentioned above). That is something
counterintuitive. However it basically means that the amplitude of fluctuations in
p(-r) — p0 calculated with respect to one particular atom in fee structure decays as
1/1" (see (3.2)). The question is why it (_lecays as 1/1‘. The answer to this question
for fee lattice and other lattices is not completely known. This problem is related to
the area of mathematics called “lattice point theory” [40, 41, 42]. The problem that is
basically equivalent to the one under consideration was originally formulated by C .F .
Gauss[43]. It was studied for more than a hundred years and it is still not completely
solved.

In our study we took a somewhat different approach. Since in nature, besides
crystalline there are also amorphous materials we decided to study at first structures

With disordered distribution of atoms in them. We a found simple exact solution for

47

 

I I I I I

3:1 “. I 11‘ 111 T
0155.-..-.14'14‘0.) 11'ka] Alma.1«1111111714111111111111111111111111de”1110111

 

 

 

 

-3

I 1 a 1 1 L I 1 1 1 1 1 1 1 1 l l 1 1

0 5 10 15 20
5 I I r I I I I I f I I I I r 1 I T I I .0
O F i I Y]
'3 — 1 l l 1 1 I l 1 1 1 1 l I 1 1 1 I I 1—

 

 

 

100 105 110 115 120

 

(30) = 490) - p,I

 

 

3_ I I I I I I I I I I I I T I I I F I I _—
F— —1
0 '1.
I- ~—
_ I-I
-3_1 1 l l I l 1 l 1 L L 1 1 l 1 1 l l l.—

 

1000 1005 1010 1015 1020

 

 

 

 

3 _ I I I I I I I I I I I I I F I I I I I _
0 '
-3 —m 1 1 1 1 1 1 l 1 1 l l ,1 1 1 1 I l 1 —
10000 10005 10010 10015 10020

Distance r, A

Figure 3.3: PDF of Ni at 15K at (-lill‘crcnt distances. The blue solid lines are results
of calculations while the red dotted lines on the two panels on top show the result of
neutron scattering experiment taken from Fig. 3.1

 

 

 

 

 

 

 

 

 

 

DJ

IIIIIIIIIIIIIIIIIITIIIIIIIIrIrrIII
3— A
F- _-
y—
A0
L“; -
0 r
-3—
g 1L41111111111111111lilillJliilllJ
0
'5
o IIIIIIIIIIIIIIIIIIIIIIIIITfifTIfiIII
C h _
:33_ B —
h F —1
t: - a
.20 ‘
H _
s l-
.-° “
E-3_11L1111111111111111111111111lilllL—]
(D
'5
8—1
'v-l
CU
Du

 

 

114111111111111111111

20 30
Distance r, (A)

 

 

Figure 3.4: PDF for the amorphous Si: (1,, 2 2.4 A, a = 0.1 A. A: with respect to 1
atom, B: averaged over 10 atoms. C: averaged over 20000 atoms.

the case when atoms are distributed completely randomly. Randomly means that
position of two atoms can even coincide. The case of amorphous materials is different
from the random case, in particular. because there is excluded volume around every
atom. We found a way to take into account the role of the excluded volume. These
two cases allow consideration in d«dimensional space and we follow this line. The case
of crystals is the most complicated and we discuss some aspects of PDF behavior on
crystals at the end.

Thus experimental PDF obtained on a crystal decays due to finite instrumental
resolution and imperfections of the crystalline structure. In amorphous materials the
decay of G(r) occurs primarily due to the averaging (3.6) of PDFS over many atoms

that have different atomic environment as shown in Fig.3.4.

49

For simplicity we discuss everywhere the case of materials consisting atoms of
one type. It is evident how to make a generalization for the case of materials that are
formed by different atoms.

Mean square deviations of the distance between a pair of atoms from its equilib-
rium value 0,2,- depends on positions of the atoms in the sample. For a given structure
the values of of, can be calculated from a particular force model[58, 36], as it was
made for the top two panels of Fig.3.3. However if the distance between atoms 1 and

.2

j is large then their motion is uncorrelated and it can be assumed that 0",}- is the same
for all atoms that are far away from each other 0,2)» = 02, as we do in this paper.

The chapter is organized as follows. At first we define PDF in d—dimensional
space. Then we consider the case of completely random atomic distribution. After
that we discuss the case of the random distril;>ution of atoms, with randomness limited
by excluded volume. Behavior of PDF in crystals is discussed after that. We conclude
by discussing how reduced scattering intensity F (q) would look like on a perfect
crystal if there would be infinite instrumental resolution. In the main body of the
chapter there are only simple ex-aluations with results of numerical calculations (on

square and triangular lattices for 2d and simple cubic. orthorhombic and fcc lattices

in 3d), while more complex derivations are presented in the appendices.

3.2 PDF in d-Dimensional Space

3.2.1 Continuous Definition of PDF

It is possible to generalize the definition of PDF in 3d to the general case of d-
dimensions. As it will be clear from the following PDF in d—dimensional space should
be defined as:

(2.0) = ,r—.— [p.101 — pa]. (3.7)

50

This definition differs by the constant factor (4w) from the definition of G(r) in 3d
given by (3.2). It is not important for us whether we define PDF with this prefactor
or without. For the purpose of generality everywhere below we will use definition of
PDF given by (3.7).

The angular averaging of (3.3) leading to (3.4) cannot. be performed in 2d or
generally in d—dimensional space in a closed form. However we will assume that

forms similar to (3.5) with correct prefactors are valid for any dimension:

'2

1 1 (7' — I'L'J')
szf') — W Z —)—: ekl3l*—Ej—;—], (3.8)
‘ 17H drag ‘1

where Cd is the total solid angle (91 = 1. (22 = 2%, Q] = 4a). In order to obtain
pd(r) from pd,(r) averaging over different points (3.6) should be performed.
In the following we will also discuss PDF calculated with respect to a particular

Site:

([1

(;(Ii(r) : "—Tipdif") _ pUi' (39)

We will see in the following that the period and amplitude of oscillations in G'Mr) at
large values of r are almost independent of r for the given type of the structure, as
can be seen from Fig.3.3. The period of oscillations in G150) is determined basically

by 0. Thus one can expect that average value of the integral:

‘0 1 ”2 .,
< '“l- r >2 ———_— GP; 1' (17‘. 3.10
do firm .<> < )

calculated over a big enough interval (R1; R2) is independent of R1 if R1 is big enough.
This approach will be used further for all types of materials.

For the case of amorphous materials. instead of integration of (15,0) over some
range one may want to average Cid-r) (wer the different sites 2'. For the case when

sites are distributed in space completely randomly it is also possible to do averaging

over different distributions. All these methods should, as it. will be shown, lead to the
same result.

For a particular type of lattice, < (132.0) > can depend only on the lattice
parameter (1. (i.e. density p0 : 1//\d) and the value of 0. Thus the dimensionless

combination of < 65,-(7') >, po and (I which is:

ate/A) =< (”Ii-(r) > 3 (3.11)

0
can depend only on dimensionless combination of po and 0 i.e. ad/po or on a/A. The
same is true for a particular type of amorphous material or for a random distribution.
In the following, instead of yaw/A) we will use the notations 95(L) for lattices.
g3,(R) for random distributions of points, {Ii-(F) for random distributions of points
with excluded volume and 931(0) for distributions of points that model glasses. Some-

times, when it is obvious what is meant, we will drop indices (1 and 2' or both without,

notice.

3.2.2 Definition of PDF Through Bins

Definition of PDF given in the previous section is a continuous definition of PDF.
This definition is useful because it. allows us to compare. model and experimental
PDFS. For our purposes, it is convenient. to give a. different definition of PDF also.

Let assume that we want to calculate PDF with respect to site 2'. Then we can
define radial density with respect to site i. at distance 7‘ by counting the number of
points in the interval (bin) [1' — (5/2, 7' +072). If there are N,((5, 7‘) points in the interval
then radial density at 'r can be defined, assuming (5 << 7‘, as:

IV) ((5, '7’)

n (50 = -——-..
p1( I) Qde—lo

(3.12)

It is clear that P<11(6~r) z p()- Then PDF is defined in the same way as before (3.9)

52

 

—I

I I I I I I I I I I I I I ]

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 b s
i: :
0 ELI-id" U1 LIILJITUUr-i LiL‘J “trial ill—[1Tb ULU‘IHJFLJ‘JT ULIJFLI “NJ-[111155
A '4 :— 1 1 1 1 l 1 1 1 1 l 1 l l 1 l 1 1 1 I —3
CL 0 5 10 15 20
CD _ I 1 . 1 , I I d
5 4 j
‘5 0 in] “fills ”11*:th “9 ull-oil ”INJWWWJHI 1115111111
2:3 ,4._ _j
I-I-u 1 1 1 1 1 1 L 1 1 1 1 1 _
g 100 105 liO 1i5 120
:5 0: flu” ”W fta‘m a firm with a) w. m W51“ -
e -4 i:— I 1 1 1 l 1 1 1 L 1 1 _
'3 1000 1005 1010 1020
Q I I I I
4; —1
0 Mb “[1”; till militia-H rifled] Vwflmd i ' _
-4 :- 1 1 1 g I 1 1 l 1 I 1 1 1 1 .5
10000 10005 10010 10015 10020

Distance (r/a)

Figure 3.5: PDF calculated for a square lattice a — 0.1(1. according to definititms
through bins (blue lines) and Gaussians (red lines). Fig.3.1

with ail-(6, 7") instead of {1111(7‘).

In order to make a con'tparison between continuous (glefinition of PDF and defi-
nition of PDF through bins it is convenient to establish some relationship between a
and 6. We find it convenient to define this relationship in assumption that in both

approaches the weight. of the peak associated with every point is unity:

exp —

. 2
1r _ r“
20“

3° 1
1 : f _—

where the height of the peak in definition through bins~h is forced to be the same as
the height. of the peak in continuous definition, i.e. h : 1/\/ 27:02. Thus (5 : 27:0 .

It was already shown (l’ig.3.3) that amplitudes of peaks in PDF cal(_‘-ulated ac-
cording to (38,39) for fee (3d) structure of Ni persisted up to very high distances.
Figure 3.5 shows PDFs calculated according to both (continuous and through bins)
definitions of PDF for square lattice in 211. Thus we see that. both definitions are
equivalent in some sense. We also see that. the amplitude and frequency of oscilla-

tions are basi tally the same at all distances.

3.3 Random Case

It is rather difficult. to understand behavior of PDF in cases when atoms or points
form ordered structures lattices. We found that the case when points are distributed
randomly, besides being easily understandalfle, also provides a great insight into the
origin of non-decaying behavior of PDF.

The value of < G5,(r) > (see (3.10)) in the random case can be easily estimated,
if it is assumed that averaging is done over different distributions of sites 2. If follows

from the definition of PDF through bins (38,39) that:

 

< [erd‘ldpd,(6, r) — erd"1(5/)0]2 >

95741—152 (313)

< Giff) >:

54

mgr-r

 

c>

U)
l

'D
u
'o
q
n
.0

C

leIII

-03

[Ill

 

 

 

 

1

.O
b)

1 1.11

I
.O
U.)

 

 

l

Scaled PDF, r (p(r) - p0) (o/po)]/2
.9
D) O
I I I I I I ’ I I l

b

i» o
ITIIII
llll

 

 

llJJIlQllIlIgllllll

100 150 200 250 300

 

r/o

Figure 3.6: Sealed PDF function g3,(r) E ng/U/IOO vs. reduced distance r/a

O“
CH

As follows from the definition of pd,(6,r) (3.12) the nmnber of particles inside the
spherical annulus of radius r and thickness (5 is given by N,(6, r) : erd‘ldpd,(6, r),
while the average number is given by 1’\",(6,r) = le‘d—ldpo. The well known result of

statistical physics can be employedHS]:
< [N — QT]? >= N (3.11)

Thus we get:

erd—ldpo 1 ,00
< C2, 1‘ >= ——.—— = ——.—. 3.15
d1< ) gird—102 d o ( )
Exact calculations for the continuous definition of PDF perforn‘ied (see text

(3. 10)) in appendices (All) in ensemble averaging and in (A.1.2) in the integral

approaches lead to the following result:

0

_ 1
2 Z [12' . _ = _
gdi(R) < C dt(r) > p0 Qfifld

= const. (3.16)
Frorn (3.15.316) it follows that in the random case < G3,(r) > does not depend on
73 This result means that. an increase in the number of particles inside the spherical
anTllllus of thickness 0 with 7' does not lead to the decay of PDF. It also means that
in tlle random case /),('r) converges to p0 as r 2 , e.1. as 1/\/7—‘ in 2D and as 1/7‘ 111

3D

Note that for the random case the value of 93, does not depend on the value of

0 / /\ .

Fig.3.6 shows dependence of scalet'l PDF g3,-(7') E G31(r) \/(7—/p: 011 scaled distance
'r/ (7- It shows that the amplitude of oscillations in the sealed PDF g3,(r) does not
depend on p0 or a. It also shows that. 0 determines the period of oscillations. Note

that the period T ~ 40 is of the order of 0.

56

 

‘PM' ' VEIMIL-n—A v- rug
I

 

 

Figure 3.7: Distributions of points for two different excluded distances. The left figure
shoxvs distribution that was obtained by choosing excluded distance in such a way
that excluded area covers 1/4 of the total area. In the right figure, excluded area
covers 1/2 of the total area.

3.4 Random Case with Excluded Volume

3.4 - 1 General Discussion

There is an important. difference between arrangement of atoms in real materials
and arrangement of points in the random case. In real materials two atoms cannot
be t‘00 close to each other. In other words, there is excluded volume around every
atofn in which there can not be another atom. Thus, in the random case, for a given
finite density, the number particles inside the spherical annulus can vary, in principle,
be’tVVeen zero and infinity. In the amorphous case, in contrast, the maximum number
of particles inside the spherical annulus is limited by excluded volume.

In real materials, the excluded volume around every particle and the atomic num-

er density of the material are closely related. However one can imagine a situation

57

 

Figure 3.8: Only those points that. are in a vicinity of the distance 7‘ with respect to
site i contribute to the value of G,(r). Excluded volume around every point reduces
the volume accessible to the other points, that may want to enter the annulus due to
fluctuations, and thus reduces the size size of fluctuations in the number of points.

When there is no direct relation between these two quantities. For example, we can
Place particles in some box randomly, but every new particle can be placed at some
DOint only if there is no any other particle within some distance. Suppose that we
alvVays want to place the same number of particles in the box, so that density would
be the same, but we want to vary the excluded distance. If the excluded distance
is 8111311 enough we would be able to achieve our desirable density. However, as we
WOUICI increase the excluded distance, it would be more and more difficult to find
the Configuration of sites that would satisfy both requirements: given density and
95in11 excluded volume. But if the excluded distance is too big we may not be able
to arChieve desired density. This problem is related to the so called problem of the
ralldom close packing of spheres[46, 47, 48].

Figure 3.7 shows an example of two distributions of sites with the same density
in each of them. However, the excluded distance in the left half of the figure was

ChOSen in such a way that excluded area covers 1 / 4 of the total area. The right half

58

 

of the figure shows distribution in which excluded area. covers 1/2 of the total area.
As it was pointed out. the value of < Gib) > is determined by the average size of
fluctuations in the number of points inside the spherical annulus, i.e. < (N — NV >.
If there is excluded volume around every point, then for any point that wants to
enter the annulus due to some fluctuation, not the whole volume of the annulus is
accessible: it is reduced by excluded volume of those points that are already in the
annulus. Below we show that excluded volume decreases the value of < 93,( E) >
con‘ipared with the random case (316).

If the excluded volume around one site is V1 f ~ Ad and there are IV 2 poV sites
inside the volume V, then the volume available for the placement of a new site is
V — V1 f1? (in the random case it. still would be V).

Since, in order for the site j to contribute to the value (1,1,:(1‘), its distance from
the site i should be ~ 7‘ i a, we can say that site j should be within the volume
I” N erd’la (see Fig.3.8). If inside this volume there is already one site, then
accessible volume is reduced by V1 f ~ yad'la, where 7 can be different for different
structures. If there are T = pOV sites inside the annulus then the accessible volume

/* - —, . _ 1
I IS reduced by ”VA compared w1th V, so that. we have:

1/* —_ v — Vlfiv ~ v — «ad-la . pOV (3.17)
1 0
r A, d—l _ ,d—l _ ,_
~I/[1—,/\ 0-—/\d]—er a[1 7A].

In COInparison with the random case for which we had (314): for the amorphous case

We Write:

< [N _ T]? >= T" = pow.

 

so that (compare with (3.15)):

 

Qd~7'd—10p0[l _ 7%] _ i&[1_ ,0]

2
< C 7‘ >= . .. —-
( ) Qj'rdTIUZ 9d a

(ii

(3.13)

Thus the presence of excluded volume in the amorphous case leads to the linear decay
of < 0,21,.(1') > with an increase of a. A more rigorous derivation made in Appendix
(A- 2) leads to (compare with (3.16)):

933512 1 ll-‘r‘il- (3-19)

2 \//,I:|S2l1

 

Results (3.18.319) were (.)bt.aine(.l in assumptions that — 157)?) << IV: and that.

 

peaks corresponding to the different sites do not. overlap, i.e. (7//\ << 1.

3.4-2 Results of Simulations in 2d

In order to verify numerically the role of excluded volume in 2d, Np = 220 x
220 points were placed in the box of size [—110,110] in :1: and y directions. The
trial coordinates (.1:,. 1),) of a particular point. were generated using a random number
generator. Point was placed at. (Jr,,y,-) if there was no any other point within the
excluded distance 2g. Otherwise a new attempt to generate coordinates of the point
Was made. This procedure was repeated until all N,D points were placed in the box. In
the 1') eginning, when the box is almost empty, points can be placed basically anywhere.
AS the box fills up it. becomes more and more difficult to find suitable coordinates
(errlpty space) in the box. The fraction of excluded volume to the total volume of
the box can be found, if the concentration c = 1//\2, as 3; /S101 = 7r§2//\2. In our
Sin'lulations /\ = const = 1.

Then, for the particular choice of 5 (distribution), for all sites that were inside
the central boa: of size [—50, 50] in .1: and y directions the PDF 02,0") was calculated

for different values of peak width 0. The average value of gift. =< (ii-(I) > a/po was

60

 

 

 

 

 

 

 

I I I I I I I I I l I IW—I I I I I I I I I I I I I I I I IT I I I I l I I I l l I I I I I I I I
l . :__ i : J; ...... A ...................... a ....................
73k ‘ . -
r— \ u i .-
l- “ g‘ . ”.‘.. l _. « : - t . - - O 05 —
0.8 — “\ i i —’
>- \\' ‘ ..
A \ " " 0.10
N " __
a 0.6 ~— \ . .
A _ . . _
I- ‘ \ 3'... l, .n
0.2 — -. __
t- " . ‘ _‘
h “.‘. .\ T h.“ 0.525 -
O 1 1 LI 1 I I l I I I I I I I I I I I I I I [III I l I I I II-J I I I I 14 I I I II. I I I l I I

 

0 0.1 0.2 0.3 0.4 0.5
Peak Width 60»

Figure 3.9: Dependence of 9%,(E) and g§,(G) on a for configurations with different.
excluded volumes. For the curves (with circles) shown, excluded volume covers 0.05 ,
0. 10, 0.20 and 0.525 of the total area of the box. The curve with squares corresponds
to the amorphous distribution of sites. It basically coincides with the curve that
has the fraction of excluded volume 0.525. The curve with triangles corresponds to
Completely random distribution.

found by averaging over the points that lay inside the central box. Our simulations
Show that g§i(E) does not depend on T, which is in agreement with (3.19). At larger
Values of r > 20 some size effects can be observed similar to those that occur in 3d
(see description in 3d section).

Figure (3.9) shows the dependencies of the ratio of ggz-(E) to the constant exact
(3~ 16) value 934R) for the random case on a/A for values of excluded volume Wég/AQ
equal to 0.05, 0.10, 0.20, 0.525.

It is interesting to compare results obtained from the random case with excluded
VOlllrne with results obtained on some modeled amorphous structure. The curve

with squares represents the results that were obtained on 3-fold coordinated network

61

 

 

 

 

 

 

 

 

 

0.2 _rrrrrrrrrlrrrirrrrIIrrwrrrrirInwrrlIIIIIHIIrrlrllri1_

I 2

: 1

_ 1

l- -1

p- —

015’— 2d —‘

. _ _.

I 2

_ _

>. - _
fl

" lI—I—r TII #1 III III "

Q) — A I l—r I l V -

8‘01? 9‘ : 3 t

m t N -- " ‘ ‘

7- N "’ '1 -4

I- on : :1 -4

I >\ 05: 1 I

I— _ -4

— e E - .

0.05:- m 0.25_ —‘ “i

_ a : : q

r- ,_ .4 —1

: DI) D’Aiiliiiliiiliiiliiii“ :

r 0 0.1 0.2 0.3 0.4 0.5 2

0:riiiiiiiilliiiiiiiiliiiauiiirliiiiirauuliiiiiiiiiLL11;—

0 0.1 0.2 0.3 0.4 0.5

Fraction of Excluded Volume néZ/kz

Figure 3.10: Dependence of the slope (linear part) of the curves 9%,. versus a (see
Fig.3.9) on excluded volume. The value of excluded volume 0.525 is very close to
the limiting accessible value. The inset shows how the ratio of the asymptotic value
(large 0//\) of 9.3,» to the value of 93,- for the random case depends on excluded volume.
Note the scale on the y-axis.

Corlstructed from the honeycomb lattice by amorphizing it with an amorphization
Procedure similar to WWW method[49, 50]. The curve corresponding to the amor-
PhO us case basically coincides with the curve obtained on the random structure with
the fraction of excluded volume 0.525. In order to obtain the data for amorphous case
to the central box of size [-25.0250] in 17 and y directions that contained 800 points
periOdic boundary corulitions were applied. The averaging of 03,0") in the interval
Of 7‘ 6 (5.0; 15.0) was made over the 800 points that are inside of the central box.
Although this approach has obvious shortcomings (see their discussion in 3d section),
it Seems, that obtained results are insensitive to them.

The curve with triangles shows the data with respect to the point at the origin
While the others 10002 points were randomly distributed (zero excluded volume) in

the box of size [—1000.0; 1000.0] in .r and y directions. The PDFS were calculated in

62

mt

 

the interval of 'r E (3.0; 7.0). The averaging of 03,0") was made over 150,000 different
configurations/distrilnitions in order to calculate g§,(R).

As excluded volume increases, the absolute value of the curves slop c," at small

values of a/A also increases, while the asymptotic value of gig-(E) reachable at large
val ues of a/A decreases. We can see that g§,(E) decreases linearly with a/A at small
values of a/A, in agreement, with (3.19). Figure 3.10 shows dependence of the slop 7."
and asymptotic (large (.7) values of 9.3,(E) 011 forbidden volume.

The precision with which we determine the slope is about 5% of the value of the
slope itself. Precision with which asymptotic value is determined is significantly lower
~ 20%, since for the large values of 0//\ there occur only few oscillations in G210")
0n the interval of the studied distances. We used as an asymptotic value of g§,(E—) its
value at a/A = 0.75.

The fact that 9.3,(E) does not (_lecay to zero at large values of a/A indicates that
even in the case when there is excluded volume. the number of particles inside the
Spherical annulus can fluctuate, although fluctuations are limited. Due to the limited
size of fluctuations, the value of 9.3, (E) in cases when there is excluded volume is
SHlaller than for the completely random case. The bigger excluded volume leads to
Slnaller fluctuations and thus it results in a smaller value of g§,(E). Behavior of
9‘3: ( E) at large values of 0//\ corresponds to the situation when peaks that originate
from} different points overlap. That is the limit when amplitude of atomic vibrations

beCOnres comparable with interatomic distances.

3-4-3 Results of Simulations in 3d

Numerical simulations similar to those made in 2d can also be made in 3d. It
1 ~ - . . . . . . .
S ll’lteresting to compare the results obtained for the random distribution of points
1' - . .
11111ted by excluded volume w1th results obtained on some structure that represents

8 .
(”he real amorphous material.

63

— u 1 H245!" “Wt—ME:
‘
i.

 

 

   

 

 

 

 

IIIIIII—rIrIIWIIIfiIWIIIIIIIrl—IFFIIVIIVIIIIIIIITI]
~ - . _
~ V... " ' p 0.02 —
0.8 - ‘ “ .. _
95 0.05
NV - . ‘
p; 0.6 — -
Lu 0.10
v ' = '
N -~ 0.4 - —«
M _ _
on _ 0.30 _
0.2 +- 0.20 —*
Glass 1
0 IIIIJILIIIILIIIIIJII‘I‘Ili-IIIIIIIIIIIIIIIIIIIJII‘I‘IIIIII
0 0.1 0.2 0.3 0.4 0.5

Peak Width, (5/).

Figure 3.11: Dependence of g§,(E) and 931(0) on 0 for configurations with different
excluded volumes. For the curves (with circles) shown, excluded volume covers 0.02,
0-05 , 0.10, 0.20 and 0.30 of the total volume of the box. The curve with squares
corresponds to the distribution of sites in amorphous Si. It basically coincides with
the curve that has the fraction of excluded volume 0.35 (not shown). The curve with
triangles corresponds to completely random distribution.

As a structure that represents a real material, we used the modeled structure
0f amorphous Si[51, 52]. In amorphous Si the average distance between the nearest
atOIns is ~ 2.5 A. Most of the simulations were made on the sample containing 20,000
atOIns that occupy the cubic box of size [—36.23; 36.23] A in :r, y and 2 directions.
PeI‘iodic boundary conditions were applied to this box. The square of PDF G3i(7‘) was
Ca*lculated for all atoms in the original box. The average value 934E) was found by
averaging of 03,0“) over 20,000 atoms. In order to find the dependence of 9%,(E) on a
Cal(llllations were made for different values of a. The described method of calculation
is not quite correct, because there are some predefined correlations in atomic positions

due to periodic boundary conditions and the fact that contributions of some pairs of

64

.9. .' ‘I'J

 

 

 

 

 

 

 

 

 

5 rfI I I I I I I l I I I I I I I I I l l I I I I I I I I l I I II
I
4* 3D “
>— 3— —
0" A l... r I I I I I I l I I I r1 _,
8‘ ‘5 Z :
m 2— '53 : : —
£1) 05; g
_ A : _I _
E: ~ 3
N ._ 0.25:— -:
l_ m >— -4 _
CD 0‘ i I l I i i i l i l I i
’ 0 0.1 0.2 0.3 “
04’: l i IIJ__LII l i i I i IIII l l I i i I i l I l i l i 1 ii
0.1 0.2 0.3

Fraction of Forbidden Volume, (4/3)1tl';3/?t3

Figure 3.12: Dependence of the slope (linear part) of the curves gi- versus 0 (see
Fig.3.11) on excluded volume. The value of excluded volume 0.525 is very close to
the limiting accessible value. The inset shows how the ratio of the asymptotic value
(large a/A) of 93,. to the value of gig, for the random case depends on excluded volume.

atoms to {13,(E) are counted more than once. However, we repeated some of our
calculations on a bigger sample containing 100,000 atoms without periodic boui‘idai‘y
conditions. From that we conclude that the obtained results are not. sensitive to the
shortcomings of the used calculation procedure.

Coordinates of the points for the random case with excluded volume were gen-
erated using a random number generator by distributing 903 = 729, 000 points in the
box of size [-45.0;45.0] in .1:, y and z directions. A point was placed at generated
coordinates if there was not any other point within excluded distance. Otherwise
another set of trial coordinates was generated. Then from those points that were
in the smaller box of size [—20.0; 20.0] (approximately 403 = 61.000) 10,000 points
were randomly chosen. With respect to them G3,(‘r) were calculated. We found that,
as in 2d, the amplitude of oscillations in 03,0) is distance independent, if size ef-

fects are ignored. PDF G3,(r) was averaged over these 10,000 points in order to find

93,—(E). The value of g§,(E) is constant at values of r < 10.0. At larger 7‘ values
g§,(E) slowly increases, which is the size effect. caused by the finite size of the original
box (finite number of points also limits the size of fluctuations). The last statement
was verified by calculating g§,(E) from the bigger box [-90.0;90.0] in 1 y and z with
1803 = 5, 832. 000 points in it with the size of the small box still [-20.0120.0]in .r, y
and 2 directions. For this bigger box the ii‘icrease in the value of g§,(E) was smaller.

Figure 3.11 shows the dependencies of the ratio g§,(E)/gf,(l{) (where, in the ratio,
the numerator represents the value that was obtained from the random (_listributioii
with excluded volume and the denominator represents the value that was obtained for
the purely random case) on 0//\ for different. values of excluded volume. The curve
that corresponds to the amorphous Si is also shown.

Figure 3.12 shows how the slopes at small values of 0//\ and asymptotic val-
ues at large values of 0//\ of the curves g§,(E)/g§,-(R) depend on excluded volume

(47763)/(3/\3)-

3.5 Crystals

3.5.1 General Discussion

Behavior of the PDF at. large distances in crystals is the most intriguing. This
case is also the most important. because only in crystals can non-decaying behavior
of the PDF be observed experimentally. Thus, if a crystal is formed by atoms of
one type and all lattice sites are equivalent (in assumption that crystal structure is
perfect) then there is no need to perform averaging (3.6) over different lattice sites
since 0,-(7‘) is the same for all of them, i.e G,(‘r) = G(r). The case of crystals also
turns out to be the most complhtatetl.

In this case we cannot say that points inside the annulus are randomly dis-

tributed, since points are fixed on the grid. Thus we can not use the same philosophy

66

 

 

 

Figure 3.13: Illustration for the Gauss’s circle problem. Sites that lie inside of the
circle of radius 7‘ shown as black filled circles, while sites that lie outside of the circle
of radius r as Open circles. Only those sites contribute to the value of (1,-(7') that are
in a vicinity of the distance r with respect to site i.

that was used for the random distribution of points to describe the case of crystals.
However, as it will be shown, the geometry of the circle leads to some similarities
between these two cases.

Behavior of PDF in crystals at large distance is related to complex and exten-
sively studied mathematical area — so called lattice points theory. The problem is
the most studied in 2d. More than a hundred years ago C.F. Gauss formulated the
following problem[43, 44]: how many lattice points of the square lattice are there
inside the of circle of radius 7‘ with the center at one particular point?

On Fig.3.13 lattice points that lie inside the circle of radius 7‘ are shown as small
filled circles, while sites that lie outside are shown as small open circles.

Assume that every site of the lattice is in the center of the square with side a,

67

as shown on the Fig.3.13. Thus, there is a correspondence between the area that lies
inside of the circle of radius r and the number of lattice points inside of the circle.
It is clear that the number of points inside of the circle can be written in the form
N('r) = NO) + h(r), where T0) = 7r(7"2/a.2), and }l('l‘) is the error term.

It is easy to see that. Gauss‘s circle problem is related to the behavior of PDF.
Let us use definition of PDF through bins (3.12). Then from the definition of PDF
in 2d (3.7), if 6 << 7‘, it. follows that:

G(r) = [27rrr5p(r) — 27r'r6p0]/(27r\/F(5), (3.20)

while from definitions (3.12) of p(r) and p0 it follows that:

27r7‘6p(r) : [N0- + (5/2) — N(r — (5/2)] (3.21)
and
27%)), = [Mr + 6/2) — We — 5/2)]. (3.22)
Thus:
G(r) : [he + (5/2) — h(r —— 6/2)]/(27r\/7_'6). (3.23)

The simplest estimation of h(1') made by C.F. Gauss[43, 44] follows from Fig.3.13.
All squares that correspond to the sites that are inside of the circle of radius r lie
completely inside of the circle with radius r + d/ 2, where d. = \/2(I is the diagonal
of the little square. From the other side, the circle of radius r — d/ 2 is completely

covered by the squares that correspond to the sites that lie inside of the circle of

68

radius 7". Thus, we conclude that:

l
—,7r(r —

u . <
(12 fl (1. (I2

i.e:
(3.24)

Formula (3.24) suggests that the value of the error term h(r) is determined by the
length of the circumference: 2m. From (3.23.3211) it follows that:

— V215; g G(r) g (5:? (3.2.5)
(I.

(I!

\-\'e see that the obtained bounds for h(/') are very weak. If 12(1') would depend
on r as 12(7‘) ~ 7' then G(r) would increase as 7' while numerical simulations suggest,
as follows from Fig.3.3 and as we will see further, that G(r) ~ cons! meaning that in
fact Mr) ~ VI.

A number of works [5351, 55] were devoted to the more precise estimation of

Mr). In 1990 .\'I.l\'. Huxley showed that [53]:
|h(r)l g Cr", (3.26)

wl’iere C is a constant. and the bounds on 9 are

L.

'6
< 6 g : z 0.63 . (3.27)

{I

mIr—i
w

The lower limit was obtained independently by OH. Hardy and E. Landau in 1915[54,
55].

If 0 would be 1/2 it would be natural to expect that amplitudes of oscillations in

69

G(r.) on average. are distance-independent. Since 6 > 1 / 2 it indicates that extreme
amplitudes of peaks in G(r) grow with the increase of distance 7'. One may argue
that. (3.26.3.2?) basically put limitations on the extreme valucs of [2(7').

Thus (326,327) suggest that [2(1') can exhibit, in principle, rather complicated
distance-dependent behavior. However. results of our mimerical simulations suggest
that, 011 average, amplitudes of peaks in G(r) are almost distance independent. It
means that 6 is almost 1/2. Being more precise, it seems that, 9 being just a little
bit bigger than 1/2 converges to 1 / 2 from above, as 'r ——> 00.

This observation is in agreement with result that. was obtained in 1992 by P.
Blelier who showed (we are not. giving here precise formulation) that. there exists a

liiiiit:

 

1 R l- 2 .
[[2210 5/0 I 'fr)’[)(I‘/R)(IT‘ = (:2, (3.28)

where p(.r) is an arbitrary probability density on 1' 6 [0,1] and (r2 is a positive
constant. This result suggests that there is also another limit that. is important for

us (see (3.23)):

1 Jr) 2 (
—.—— ; (1‘: *2. 3.29
13.2 — It, /,,, (r) ' a ( )

In the following sections, after the. discussion of the exactly solvable 1d crystal, we
present some puzzling results of our numerical simulations in 2d and 3d. Then we

return to the general discussion again.

3.5.2 Exact Solution for 1d Crystal

According to definition of PDF in d-dimensions (37,3839), in 1d:

G(r) = ()(r') - pm (330)

 

IIIIIFIIWIIIIIIIIIIJIIIITIIll IIIIIIIII

 

 

 

 

l :— r
T i
r- ‘1
— 4
l——- —
f :
_ _
>— _.
I :
[—— _—
[_ —.
g , s
)— _.

0 I l l 14 l l l l l l I II I I l l I II I I I l l IT‘LI I I 1
O 0.1 0.2 0.3 0.4

Peak Width o/a

Figure 3.14: Figure shows the dependence of the ratio gfL(a) / 9123(0) on 0. At small
values of a, when different peaks in G(r) do not overlap, this ratio linearly decreases
as or increases. As peaks start. to overlap (large values of a), the rate of decay in the
average amplitude of PDF decreases.

where

 

p<r>= 1 Semi—W] (3.31)

and )00 = Me.

It can be shown (see appendix (A.3)) that for any value of a:

 

 

Girl”) 0 00 71202 .
2 ) = [1 — 2\/ffg] + :exp [— 02 ]. (3.32)

This curve is plotted on F ig.3.14.

The last term in (3.32) can be ignored at small values of 0. As 0 increases it
changes the the rate of decay in the average amplitude of PDF. Thus, behavior of
the scaled amplitude gfL(0)/gf(R) for 1d crystal is similar to the its behavior of in

glasses (39,311).

71

3.5.3 Results of Simulations in 2d

In order to demonstrate that on average the amplitude of oscillations in 12(1')

- '0
indeed scales as 7"“

we calculated explicitly for the square lattice the number of
lattice points inside the circle of radius r E (0; 100. 000)u with step 0.010.. where a is

the lattice spacing. \Ve define function F(r) as [56]:

 

F(r)=" ' ‘ (3.33)

If h(r) ~ N” then the amplitude of oscillations in F(r) on average should be constant.

Insets on the Fig.3.15 show F(r) vs. T dependence in three different ii’itervals
when r E (0; 100, 000).

The values of F2(r) were averaged over 100,000 different values in a. few ii'itervals
of length 1000a. The average value of < 172(7) > as a function of r is shown on
F ig.3.15 as a dashed curve. Circles are plotted in the beginning and in the end of
the corresponding averaging interval. From the vertical scale it follows that while r
changes from 0 to 100, 000 the average value of 172(1) changes by less than 1%. From
this we conclude that basically h(r) ~ 7“”. There is also a slight increase in the
value of < F2(r) > with r. This increase becomes slower as 1‘ increases. From that
we conclude that. 9 converges from above to 1 / 2 as 7“ increases.

An example of PDF calculated for the square lattice in a continuous approach is
shown on figure (3.5) as a dashed line. PDF calculated for triangular lattice exhibits
similar behavior. In order to calculate g§( L) (see (310,311) for the square and trian-
gular lattices as a function of 0, the squares of corresponding PDFS were integrated
(see (3.10)) in the interval of r E (1000)\;2000)\). We should point out here again
that the values of g§(L) are very insensitive to the position of the averaging interval.
Calculation in the interval r 6 (300A, 600A) leads to almost the same results.

Figure 3.16 shows how, for the particular type of the lattice, the ratio of g§(L) to

72

 

 

 

 

 

     

 

 

  
 

 

 

 

 

 

Plllllllll]IIIIITTIIIIIIIIIIIITTIIllllllllllllllllq’
: .0 ........... my ..
_ .' """"""" 116.
2.575 _ ,. —
Z ..... -°°'. Z
_ .0. :
._ 3 " .-°° _
- es .. :
2.57 -— .3? —
_ a _
A : I :
A :

NV 1' § ‘
LLI _ as a
V - 0 200 400 600 800 1000 -

2565 __ on 5 - I I I l I I __

_ 3 O _

i 5 Le , . . I

_ _5 11’ l a 1 l I f ' l r _

_ 1 l l I l a

”‘36 9000 9200 9400 9600 9800 l 0000 ‘

-E 5 I l I I I I 1

2.56 j , . , P —

I; 0 I

:3 _51' I‘ll ‘ H l . ”l ‘ I l .1 1w... :

t-: 1 L l 1 l 1 l l n

: 99000 99200 99400 99600 99800 le+05 :

2 555E1 L11111111111111L11111111111111 1111111111111111 1-
' 0 20000 40000 60000 80000 1 e+05

Distance r/a

Figure 3.15: Inset shows how F(r) depends on r for r E (0;1000)a, r E
(9000; 10, 000)a and r 6 (99,000; 100,000)a. It shows that on average amplitude
of oscillations in F(r) does not depend on 1'. At least, it is impossible to see it by
eye. The dashed curve on the main figure shows how < F 2 (r) > averaged over few
intervals of r of length 1000a depends on the position of the averaging interval. At
the beginning and at the end of every averaging interval circles are plotted.

73

 

I
. _ . I . . . . ..
—\ . 0 . _\y . v I .r. m1 . .- «a... l
..)l... lam .Q WI. I. u. III. 1
l ...c . I . l L t .0
YD} h... lniih rl no: I s
It . .ylg. . , Ill.‘ Fit . A I
1 -Nu . MI 1 ... rI In .\.. 9. P. .I. . ..
. \. .w ..I . ,. .1 l f. LI“ .I . .w
A N.— v nm\A n vi u _.:.t ,fw‘. .h.. s-.. . Pl Us h. a? 3L
. NU AI. PL . .M .. s L rolx Iév v 11194 1.. n

 

 

 

 

g22<L)/ng(R>

 

 

 

 

Square Lattice ~
0 i " l [U I [D I T
0 0.05 0.1 0.15 0.2 0.25 0.3

(50»

Figure 3.16: Dependencies of gfi(L)/g§t(R) on a/A for triangular and square lat—
tices. Triangles and squares represent the results of numerical calculations for
triangular and square lattice respectively. The solid fitting curves are given by
formulas g§(t'r2')/g§i(R) = 2.483(a/A)'1/3 — 18.264(0//\) for triangular lattice and
g§(sqr)/g§i(R) = 1.849(0/A)'1/314.478(0//\) for square lattice. Dashed lines high-
light the effect of excluded volume. Expressions for them are given by the first.
divergent term in the formulas for the solid curves. They do not contain the second
term (that might be) caused by excluded volume. From the figure, it follows that the
effect of the excluded volume is bigger for triangular lattice than for square lattice.

74

the g§,(R), where 9.3,(R) is the exact constant value in completely random case (see
(3.16)) depends on a//\.

As it was already pointed out {13(L) being dimensionless can depend, for the
particular type of the lattice, only on a dimensionless combination of po and a i.e. on
p002 or on a/A (see (3.11) with related text). Thus the curves presented on Fig.3.16
are universal for triangular and square lattices.

From (316,319,332), Fig.3.9,3.1:1 and Fig.3.16 it. follows that. behavior of g§(L)/g§,(li)
vs. a in case of 2d crystals is very different. from its behavior in the completely ran-
dom case or in the random case with excluded volume (the amorphous case too),
or in the case of 1d crystal. In the random case this ratio is unity by definition
(93,0?) = const). In the random case with excluded volume it decays from unity (at
small values of a) to some finite smaller value as 0 increases, as it does for 1d crystal
also. However, in the case of 1d crystal this ratio decays to zero. In 2d crystals
93(L)/g§,-(R) diverges at. small values of 0 and it seems that it is decaying to zero at
large values of a.

We tried to fit numerical data with analytical curves in the form:

)"’ — c-)( ), (3.34)

where 01,02 and T] are positive constants. The first term was chosen in the simple
form that can provide divergence of g§(L)/g§,(R) at small values of a. The second
term originates from the notion that in the crystalline cases, in some sense, there is
also excluded volume around every lattice point. Thus it was chosen in the same form
as it is in the random case with excluded volume or in the amorphous case (3.19)

In fact the form (3.34) provides a rather poor fit in 2d. However, the same form
provides an extremely good fit in 3d case, as we will see.

For triangular lattice the values of coefficients that can provide the best fit in

75

Iit ll'lailt' ‘ 5

.‘ v
v . . ~ o
.IIH‘JIII? I1 -
t

 

l ‘
l‘ (.54.... 4
It: $144.th I"

 

:14-
-‘I
.
.n .
(hf-“is.
' ‘RI
F1
l.
.
.
, .
-.

the whole studied range of a are: cl : 2.483, r] = 1/3 and (:2 : 18.261. However, at
small values of sigma, different values of coefficients provide a better fit: c1 = 3.712,
T] = 0.245 and (1'2 : 413.433.

The situation with the square lattice is similar. The values of (.‘.oefficients that
provide the best fit in the whole range of a are: c1 = 1.819, T] = 1/3 and c2 = 14.478.
The values of coefficients that provide a better fit at small values of a are: (:1 = 2.806.
I] = 0.215 and c3 = 33.4110.

Note that the value of the coefficients c2 for triangular lattice is bigger than for
the square lattice. That is in agreement with the observation that triangular lattice
is more densely packed than square lattice. 111 other words, the ratio of excluded
volume to the total volume of the sample for triangular lattice is bigger than for
square lattice. Thus from F ig.3.9 it follows that (:2 for triangular lattice. should be
bigger than for square lattice, as it is. This observation can be considered to be an

indirect indication that. the form of the second term in (331) is correct.

3.5.4 Results of Simulations in 3d

An example of PDF calculated for fcc lattice is shown on Fig.3.3. In order to
obtain the values of 93(1.) as a function of 0 the square of PDF was integrated for
fee and orthorhombic (b/a = 2.0, c/a. : 3.0) lattices in the interval 1' E (300A; 300A +
20000), where [)0 = 1//\3 ((‘Ibtained results are insensitive to the position of the
averaging interval as in 2d).

Figure 3.17 shows how the ratio g§(L)/g§,(R,) depends on 0 for fee and or-
thorhombic lattices. This ratio in 3d also diverges as a ——> 0, as it does in 2d.
We tried to fit the data obtained from numerical calculations with the. fitting curves

in the form similar to the one used in 2d (3.34). We found that. c1.1rves:

M = 0.78(

93le

)‘1 — 200(3) (3.35)

a
/\

>’

76

 

 

   

 

 

 

 

 

 

 

 

 

 

90 I I I I I I I I I I I I I I I I
_ .. A 15 3 _
80 _— “ fcc : 7
_ A " ._
70— 10 ‘— —
_ 5 _________ ‘5 -
60— II _ __
ax ' ..
NV 50— —
3:.
OD ' ._
>\
a 40 ._ - —
NV
of)" I . ‘
30 _ -—
2O _ ‘ 2‘ '—
10 ’— i ‘ i 7““ ......... A
orthorhombic ’ ' " = -
0 t- _
l 1 l l 1 1 1 l L 1 1 1 1 l l 1
0 0.05 0 1 0.15

Peak Width o/k

Figure 3.17: Results of numerical simulations for FCC and orthorhombic lattices. For
rectangular lattice b/a = 2 and c/a 2 3. Points represent results of simulations. Solid
fitting curves are given by formulas described in the text. Dashed lines highlight the
effect of excluded volume. Expressions for them are given by the first divergent term
in the formulas for the solid curves. Thus they do not contain the second term (that.
might be) caused by excluded volume. It can be seen from the figure that the effect
of excluded volume is bigger for FCC lattice than for orthorhombic lattice.

77

I" ’,__ . . o o .
ILJ ELI. idx . .

 

Emttm:.

i I_...
.'--I 11.1 ldi

'."lI"'
n~haru

II}: III

p
.. .

J x

‘1‘ '7‘! v
a 2...:

90
DU
0'

_LIN

 

for fee lattice and

(1.3% (W)
UIIIIH)

 

) (3.30)

for orthorhombic lattice provide a. good fit to the numerically calculated data.

Note that the value of the coefficient. c2 for fee lattice is bigger than for or-
thorhombic lattice. It is. as in 2d. in agreement with conception of excluded volume.
Since fcc lattice is packed more densely than orthorhombic lattice, excluded volume
for fee lattice is bigger than for orthorhombic lattice. Thus according to Fig.3.11,3.12
it should be 02(fcc) > c2(0rf).

We would like to highlight again that the curves plotted on Fig.3.17 are universal
for a given lattice, i.e. the ratio g;f(L)/g§,(}?) indeed can depend only on the ratio

0//\.

3.5.5 Speculations on the Origin of a-Divergence.

In order to find 513(1.) it is necessary to calculate an integral (3.10):

'7 1 H2 4),
< Cflr) > : m /n Gg(r)dr (3.37)
,

= -———R9 _ RI R {[Q()I'd"lop(r) — 9.17"!"la/JOIQ/[Qird‘la‘z]}(1r .

Since the value of the expression above (almost) does not depend on position
of the averaging interval, the value of the last integrand should also be on average
r-independent. We know that, in the case of random distribution of points, the
fluctuations in the number of points inside of annulus is determined by the volume of

the annulus or by the number of points in it (see (3.13,3.14,3.15)):

[erd‘lap(r) — Qd'r“l“10p(,]2 ~ Qd'l‘d‘lapo. (3.38)

I ,~ I .-
II... I

gif'lll‘i‘b {II

 

 

Thus
[IflMHS
FIII‘III
a: 13.35.1i
.1. Ari 3

Have

 

Thus for the random case we would get:

0 .3 1
QWF<%M>3~—=mm. mm)
.00 ad

In (.'rystalline case, when lattice points are fixed on the grid, it is natural to assume
that the size of fluctuations is determined by position of the surfaces enclosing the
annulus and thus basically by the surface area of the annulus. In principle, fluctuation
on internal and external surfaces can be correlated. The size of correlations can be

a-dependent. Thus. for the case of crystals, we write:
[Qiil‘d710/Jfl‘) — 9.11"!"10/IOIQ ~ Qd'l‘d—l/\[)07](U/)\). (3.410)

So that for crystals instead of (3.39) we would get:

f, a 1 0'
fieb<QM>—~—q

[)0 (2d )‘17]d(a//\). (3.41)

Thus, if for the random case we had 93(H) ~ coast, for the ordered distribution
of points, we can get. divergence if, for example. 77(a//\) ~ co-nst.

Fitting curves for fee and orthorhombic lattices in 3d (see Fig.3.35 with formu-
las (335,336)) suggest that ‘7]3((7/)\) ~ const. 111 2d in order to obtain divergence
(0//\)‘1/3 (see Fig.3.16) we should assume that. '7]2(0//\) ~ (U/A)2/3.

However, everything that is written in this section should be considered a hypoth-
esis (we are not specialists in lattice point theory) and further investigation of this
rather complicated problem, in connection with Gausss circle problem, is obviously

Il€C€SSEtI§C

79

3.5.6

IIr r:;4"|
5, ‘\A41

Q

-..,.. -r
.‘..lilull
Tl
‘

I !
'wa‘] ..
Lu“ 14(1

.H. . .
bx, ._ ~
~. unto“

vav‘ I
.‘.
. Q ‘
.
'v‘..‘.
-ALI¥,‘1-§

tyrant-i, ‘
A4 1'

ski
‘ .

 

3.5.6 Order and Disorder. Similarities and Differences.

III behaviors of PDF on crystals. on completely random distribution of points,
on random distributicm of points with excluded volume and on glasses there are some
sin'iilarities. as well as differences.

Thus, in all of these structures amplitudes of oscillations in PDF that were
calculated with respect to one particular site persist (or almost persist for the case of
crystals) with an increase of distance. We understand the origin of this behavior in the
case of completely random distribution of points and in case of random distribution
with excluded volume. The case of glasses is very similar to the case of random
distribution with excluded volume. We do not really understand the origin of this
behavior in the case of crystals.

We also understand the dependence of the ratio gi(E)/gj,(R) on 0 in the random
case and in random case with excluded volume (amorphous case too). The case of
crystals again represents a puzzle.

As it was already pointed out. behavior of PDF in crystals is closely related to the
behavior of the error term in Causs‘s circle problem that remains under investigation
for more than 150 years. It is interesting that this behavior that was (almost[57])
an abstract. mathematical problem can be measured, in principle, in a real physical
experiment.

It is clear that the size of oscillations in 0.1.0“) is determined by the size of
fluctuations in the number of points near the surface of d—dimensional sphere of radius
1‘. That is also true with respect to the behavior of h(r) in crystals, since points deep
inside of the sphere are not. sub ject to fluctuations. However, it is hard to say that.
points that are fixed on a grid participate in any kind of fluctuations. Thus the origin
of these fluctuations, in case the crystals, lies in geometrical compatibility of the
surface of the circle/sphere with its non-zero curvature and geometry of the lattice

whose edges / faces have zero curvature.

80

-.‘l‘ 1

 

. , “
infinite

' r. .
'r F‘ I! .
“Auk-Atf"

l
4. I g‘
1 " I
>- I‘l
It. “LII; I"
V
The.
A.)A
‘P .v, 1.2.7
[Afl“)‘l‘il::
“P“? r_, “
Ava-..l1\ ‘\‘
I
1’1»

u ,‘
mpg, «.1
TI
Ib,1‘l' ’1
1‘" 3 ~h.t
FW-

«91

3.5.7 The Role of the Spherical Geometry.

Spherical geometry is basically imbedded into the definition of PDF due to its
connection with diffraction experiments. However. in order to illustrate the role of
spherical geometry, it is also possible to define tritmgular density /),,.,(r) or square
density psqr(r) (and so on ...) by counting the number of points inside triangular or
square annuluses that are shown on Fig.3.18. Here we assume that PDF is defined
through bins. The “re-Idius” of the triangular density or square density can be defined
as the radius of the circumference to which this square or triangle is inscribed.

It is almost obvious from Fig.3.18 that square density on the square lattice ba-
sically turns 2d problem into 1d problem. Thus PDF defined according to (3.7), as
G(r) : fiIper) — p0], will diverge as 1‘ increases. In order to avoid divergence
PDF should be defined, as in Id. as G(r) = [p,,.qr('r) — p0]. It. is not so evident. but it.
is also true with respect to triangular-pdf. Figure 3.19 shows that the amplitude of
oscillations in [pm-(r) — p.,] persist. as 7‘ increases.

It is easy to see that if points are distril'nited randomly then PDF should always
be defined (for any polygon) according to (3.7): G(r) = fiI/IU) — p0]. On the
contrary, if points form a lattice then for any polygon with finite number of edges
PDF should be defined as G(r) : [/),,.q,.(r) — p0]. But if the number of edges becomes
infinite (circle) then PDF for the lattices should be defined in the same way as it is
defined for the random distribution of points.

Thus it is the transition from a finite number of edges to the infinite one—
transition from zero curvature of the edges to non-zero curvature of the circles circum-
ference that makes random and ordered distributions of points somewhat equivalent.
Appearance of this equivalence changes the way in which PDF should be defined.
Thus the role of the spherical geometry (non-zero curvature) is extremely important.

for the properties of PDF.

81

 

 

 

 

 

 

 

 

Figure 3.18: Triangular or square PDFS can be defined through triangular or square
densities, i.e. by counting the number of points inside triangular or square annuluses
and dividing this number by the area of the annulus.

82

 

3 Ifi IHI ‘91 I I T I I I I I I I l I I I I I I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.. - _ _ , L
J I I I I
1...; -r‘
0... -i
H mm MWWWWWW
1 11111111111111 1111 1 1 1 1
0 5 10 15 20
3’- I I I I I I I l I I I I I I I I I 1 I“
O 2" I F r I“ I' F H I‘ I -‘
Q. I- a
g 1— _
it ~ _
a? 0— WI —
-lillllllllllllllllllllllllllllllllllllllllll llllllllllllllllllll
1 1 1 1 11 l 1 l1 1 1 1 l 1 1 11 1
80 85 90 95 100
3~I 1 T I I T I T I I I I I l I I I I 1‘
2_ I II n Ir
1*- IF
I— .
O— ...—
_11IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII llllllllli
1 1 1 1 1 1 1 1 1 1 11 14 141 1 1
480 485 490 495 500

Distance r/a

Figure 3.19: Dependence of pm — p., on distance r/ a for triangular PDF on square
lattice. Since the amplitudes of oscillation persist as distance increases the figure
demonstrates that the use of the polygons (any polygon can be divided into triangles)
instead of the circle changes the way in which PDF should be defined.

83

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

A) L=100A. 0:0.05A;
LL. 1 l 1 1 l 1 l
.g —1
g B) L=200A. 0:0.05A_:
ED 1 1 l 1 I; J L -4
'c
2
8 C) L=100A,o=0.1A “
m i
'U
0.)
g l l 1 l 1 . l *
'0
Q)
m h —1
ZOOFJH [J I ll 1 D) L=200A, 0:0.1A fl
0 +- T Win [wil'milir‘lfl‘f'i‘A‘JL ‘
-200 1 l 1 l 1 l 1 l l l 1 l
0 IO 20 3O 40 50 (,0

-1
Momentum Transfer q, A

Figure 3.20: Reduced scattering intensity obtained by Fourier transformation from
PDF calculated for the values of lattice parameters similar to those for the Ni crystal
at 20K and 300K: (1. :~_» 3.53A, ago = 0.05Aand 0300 = 0.1A.

3.6 Perfect S (q)

The quality of the measured scattering intensity and, as follows from (3.1), of
experimental PDF is limited by the finite instrumental resolution. However, one
may want to calculate “perfect” scattering intensity, as it would be if instrumental
resolution would be infinitely good. It. can be done by performing inverse Fourier
transformation of perfect PDF that can be calculated if lattice parameters and prop—

erties of atomic vibrations are known.
[{rnar.
F(q) = / G(r)sin(qr)dr . (3.42)
0

But, in the case of crystals, PDF does not decay at large distances and it becomes
unclear at what maximum value of anx integration should be terminated.

Figure 3.20 shows the reduced scattering intensity F (q) calculated from PDF

84

’mv'q ‘
. .‘.

o 1.4.

.
' tr _
‘ II

(A -V

v ,? ‘v

‘I rm' 4

AA.. ~‘(A

91. _7

Lglt‘ \Yv
. ..

,

tnll

"Al. “I.\

I 1‘
ll. ( I
,,‘

V ,_
.

‘,

'i$.\“!r
_ ‘.

function for the lattice parameters of the Ni crystal. It was assumed that atoms of
Ni arranged in the perfect FCC lattice and that a is the same for all atomic pairs.

The panel A) shows F(q) calculated for the value of 0 = 0.05 A that is close to
the value of a for Ni at. 20K. Integration was terminated at the value Rm” = 100
A. The panel B) shows F(q) calculated for the same value of 0, but integration was
terminated at Rum. = 200 A. Panels C) and D) both were calculated for the value
of a = 0.1 .A which is close to the value of 0 for Ni at room temperature. The
intri‘gration was terminated at Hm“. =2 100 A for panel C) and at Rm“ 2 200 .A for
panel D).

Thus we see that although PDF does not decay. F (q) always decays. Moreover
(2mmr at which occurs decay of F(q) is determined by (.7 and not by Rum)... However,
Rm“. affects the amplitude of the peaks: compare vertical scales on the. panels A)
with B) and C) with D). This behavior could be easily understood. The width of
the smallest feature in PDF is given by 0. Thus if q is that big, that sin (qr) makes
full oscillation on the length ~ ‘20, contribution of any feature in PDF to F(q) would
almost vanish. Thus (2,,m.20 ~ ‘27r and Qmmr ~ 7r/0. Thus increase in Rm“ does not
lead to the increase of QM”, but develops the structure of F(q) in the fixed range of
q between zero and Qmar.

This note is important. in the sense that sometimes it is assumed that an increase
of (2mm, in experimental measurements would lead to a better quality of PDF. This
assumption, in general, is wrong. If QM! ~ 7r/0 then in order to obtain high quality
PDF it. is more important to imprtwe instrumental resolution in the range q < Qnm.

than increase QM”.

85

v
., v-Y'

Emil-“l ll 1
All? Ur i2.
1' f ‘ 1 ‘
.11 cc.‘ 1

«.~ .'. 7'

u.....; as

 

V
0. Y‘
.l.‘“Pl
U
I
vv .

‘c

3.7 Conclusion

The PDF with respect to particular atom can not be measured by any presently
known technique. In any experimental measurement. the averaging which places all at
the origin is automatically performed. Thus, the pair distribution function (’ralculated
in assumption of infinite instrumental resolution only decays at. large distances if
different atoms at. the origin have different atomic environments. Therefore for perfect.
crystals the pair distribution function does not decay. However even for perfect
crystals the experimental PDF obtained from the diffracted scattering intensity will
eventually decay at large distances due to the finite instrumental resolution. This
means that the measured rate of decay of the measured pair distribution function
for very good crystals could be. used to test instrumental resolution. The decay of
the PDF at large distances can also be used to investigate the sizes of nano-crystals
and strained regions. The 1i)ersistence of fluctuations in the pair distribution function
has a. connection to a value of the rms fluctuations for the Gauss circle problem in
d-dimensions, thus providing a link between ordered and disordered distributions of
sites. From my point. of view it is important to have better understanding of the
origin of this C(mnection. Thus. this problem requires further investigation. However.
I do not believe that we (by ourselves) would able to move significantly forward with
this problem due to the absence on sufficient knowledge in the area of the lattice

point theory.

86

C]

i3

N

V.

s;

4‘

N.
-.- ‘J‘.

x

I
A

Chapter 4

Quantum correction to the Pair
Distribution Function calculated

classically.

4. 1 Introduction

Recently, due to the subtleness of the protein folding problem especially, there
is a demand for different techniques that will allow an accurate study of molecular
structures and their dynamics, yet will not. require time-consuming calculations, at
least for relatively small molecules.

The pair distribution function (PDF), that can be obtained by Fourier trans-
form of powder diffraction data, traditionally has been used to describe short-range
correlations in atomic positions.

In recent years a technique has been developed that allows one to achieve an
extremely good agreement between calculated and experin‘iental PDFS for crystalline
materials. [58, 36, 37].

Molecular Dynamics or Monte Carlo techniques are usually implemented in order

87

II ‘
-'v I
, ,,
.‘.H ‘i

.
O
.. -\
'L .
¢
r‘\ H—
...c
.
‘
y .1 1
- :.
‘Aaa A
"u; . ‘
»'.l‘ 4
.
V
'..." l.
.
., r '
“-1 11

 

 

"F F
“ 4‘ 51
‘V
4'
' §‘-
‘JA.‘
1!“

 

 

to study molecular structures and their dynamics. Thus. quantum effects are. usuallv
ignored. Other approaches. like the Car—Parinello method [59], that are l.)ased on
quantum mechanics. are more time consuming than classical approaches and only
feasible for a small number of atoms.

Accurate experimental PDFs can be obtained from X-ray or neutron diffraction
experiments. Thus. in order to use PDFS to compare modeled molecular structures
and their dynamics with real molecular structures and their motion, it is necessary to
calculate molecular PDFS accurately. It is possible that PDFS obtained in classical
calculations, in ignorance of QM effects, do not provide sufficient accuracy any more.

Thus, it is important to estimate. the role of QM effects in calculations of PDF
for molecules. Further, we will see that the accuracy of PDFS calculated using clas-
sical methods can be significantly improved without switching to QM calculatirms
completely, but by implementing some quantum corrections to classically calculated
PDFs

It is known, that even at zero temperature, if it would be accessible, atoms would
vibrate near their equilibrium positions due to the conmletely quantum effect that. is
called zero—point motion [(52, 63]. This effect is completely missing in classical MC or
MD calculations.

AS was already discussed in the previous chapters, the width of the peaks in the
PDF of solid materials is determined by the mean square deviation of the distance
between a pair of atoms from its equilibrium value. Thus, at low temperatures es-
pecially, ignorance of zero point motion in classical calculations should result in a
peak-width that is (significantly) smaller than the one that occurs in reality.

In this chapter the role of atomic zero point motion is discussed. The size of the
effect is estimated by comparing the dependencies of the mean square displacements
of a particle in the Morse potential on ten‘iperature in QM and CM calculations. W)

used a Morse potential because it was developed and proved to be useful in modelling

88

 

 

“IAN.

A

‘-_

the properties of diatomic molecules. This potential is also convenient because it
allows an exact QM solution. After that we develop a technique that. allows to take
into account the effect of zero point motion in PDFS calculated classically for more.
complex molecules. The same method is used to correct the classically calculated
heat capacitance.

Usually, even in very simple molecules, there are different energy scales associ—
ated with the borid-stretching, angle-bending, dihedral angles rotations and so on.
Although our correction can be applied to all pairs of atoms, it. is the most important
for those pairs that strongly interact with each other so that the energy associated
with this pair (at a given temperature) is still almost the ground state energy of the
corresponding vibrational mode. Thus, our correction is the most important for those
pairs of atoms that bring sharp peaks into the PDF.

Molecules are usually in some enviromnent that also contributes to the value of
PDF at a particular distance. Often there are many molecules that form a liquid.
Intermolecular interactions are usually much weaker than intramolecular interactions
and thus they can be. treated mostly classically, i.e. with MC or MD technique. In
order to calculate accurately the PDF for a. real sample it is necessary first of all to
be able to calculate accurately the PDF for a single isolated molecule. That is the
question that we address in this chapter.

Thus, our technique suggests the following algorithm for more accurate calcula-
tions of PDFS. The initial PDF can be extracted from classical MC or MD simulations
that provide coordinates of atoms in a molecule as a functions of time, i.e. molec-
ular trajectory. This modeled PDF already can be compared with the experimental
PDF. Then, in order to achieve a better agreement between modeled and experimen-
tal PDFS, quantum correction can be applied to the PDF obtained from classical
simulations.

In this chapter, we at first. consider the case of exactly solvable Morse potential

89

‘. I
l‘ .'a' ', -l
,4-‘ ..u-
[fail (”a
r} ll 1‘
x I »1\‘
f9
grrfy
nip-,4

,.‘u

(27 ~
..IA

.

If? “'V.

.‘ -AA|

9.ng

for a diatomic molecule. We use this example to demonstrate the precision of classical
calculations of PDF and the importance of the quantum correction. Another quantity
that can be corrected with our method is the heat capacitance that was calculated
by classical methods.

After that, in an attempt to use our correction for more complex molecules, we
apply our method to the C6111.) molecule. At first we study the role of quantum
corrections at very low temperatures when the molecule is in one of the equilibrium
configurations and it can be assumed that the atoms only slightly vibrate about their
equilibrium positions. If the an‘iplitudes of atomic vibrations are small enough, it can
be assumed that their motion is harmonic. In this case the problem could be solved
exactly using either classical or quantum approaches. Thus, the role of quantum
effects can be estimated at different temperatures. Finally, we calculate the PDF for
the 06H 14 molecule at higher temperatures when the potential cannot be considered
as harmonic any more and at even higher temperatures, when the molecule becomes
flexible.

In our calculations we used the “TINKER” molecular dynamics simulation pack-
age that can search for equilibrium configurations and calculate eigenfrequencies and
eigenvectors of molecular vibrations in these configurations. It was also used to run
molecular dynamics simulations at a particular temperature. The obtained molecular
trajectories (coordinates of atoms as a function of time) were used in calculations of

PDF.

90

 

'llflun. «mind; r , .

..r.
\rnI
J» .' .

v
71‘? ,
”beam.

'1
Q

~
‘I‘IV,

.a‘rj‘l

h ' n-.
'J'». .9"

4.2 Pair Distribution Function for a single molecule

The experimental PDF function (},._,.(r) is obtained from the powder diffraction

data via a sine Fourier transform of the normalized scattering intensity 9(IQ)
Ix.
cum = — / ewe) — 11<Qr>de (4.1)
o

where Q is the magnitude of the scattering vector. For elastic scattering Q =
:17.’ sin B/A with ‘26 being the scattering angle and A the wavelength of the radiation
used.

On the other hand, the PDF is related to the structure of the material. The
PDF is simply the bond length distribution of the material weighted by the respective
scattering powers of the contributing atoms. It can be calculated from the structural

model using the relation:

 

4m? < b2 >

I f),'f)j ‘
cw) = .im-{[ Z Z ———6(r — 731)] — p0}, (4.2)

where the sum goes over all pairs of atoms i and j swarated by "iij within the model
sample. The scattering power of the atom 2' is b, and < b > is the average scattering
power of the sample. In the case of neutron scattering b, is simply the scattering
length, whereas in the case of X-rays it is the atomic form factor evaluated at. a given
value of Q. For Q = 0 the value b,- is simply the number of electrons of atom i.

The first term in square brackets in (4.2) represents the radial density [)(r) nor-
malized in such a way that the number of atoms in the spherical annulus of thickness
(11' is given by 47r7‘2p('r)(1r. In solid materials, where atoms vibrate near their equi-
librium positions 6(r — '1‘,J-)—-functions in (4.2) should be substituted with Gaussians,
whose width 0?]- is determined by the mean square deviation of the 7“,}- from its equilib-

2

rium value: at) =< (7",1- — f”)2 >. For our purposes it is convenient to define, instead

91

Hf]

-'.»

W? l' til

‘. c .
L‘A'. a)"
.

Y!“ ‘ y
...I’I'. 'I

.'. .;
Jail
l-1I.'l

.:‘l‘.l
A _ ‘*,l

of Gc(r), a different function g(r) to which we will refer to as the PDF. Sunnnarizing

we can write:

 

(‘)-1r='2 o—z: “’1' 1 ,.. [ """f'JVH (43)
g! _- Mp, _‘ j <b>2\/‘2:t2mp 20'?) i
1 H U

In order to accurately calculate PDF it. is necessary, in particular, accurately calculate
the peak width 0?]- for the atoms that belong to the same molecule. That is the issue
that we discuss in this work.

Substitution of (5-functions with Gaussians occurs due to the time averaging of
the distances between atomic pairs. If vibrations would not be small, for example if
a molecule becomes flexible, (5—functions should be substituted not by Gaussians, but
by normalized probabilities to atoms in the pair at a given distance. For a particular
pair, the shape of the probability function can be rather complicated. For example, it

can contain several peaks that. correspond to the different equilibrium configurations.

4.3 An example of Morse potential. The idea of

the method.

Assume, for simplicity, that we are interested in an accurate calculation of the
PDF for a diatomic molecule. Also assume that the only variable on which potential
energy U depends is the deviation 1: = r — 1'0 of the distance 7‘ between the atoms
from its equilibrium value To. If there is only one minimum in U = U (:r) then the
molecule has only one vibrational mode.

In the center of mass the Hamiltonian of the system can be written as:

.2

[I = £— + f/(.I:), (4.4)
Zn

92

’xl'lt If? I

I
... -.v
III", Ar
lau;z-1..\

‘1’ II”;
M.

' ' .
Lox ‘. ‘A

' I
o, i: 1‘

‘7 ta . . [A ,

~ ,.
nr- 5
Li“ l

A -
.

L)’;

Y

’11

' 1

where ,u = m / ‘2 is the reduced mass.

In the limits of low and high temperatures. this problem could be solved approx-
imately. At low temperatures. when the deviation of the distance between the pair of
atoms from its equililn'ium value is very small, instead of the potential U (.r) we use

the harmonic approximation:

 

A) a)
1‘ ‘“ .77 — .l‘ “
[w ( U) a (.1 I')

L’H(.I‘) 2 I + 9

SO that we can write:

I] _ [)2 + /I.w2(.r -— .1:“)2 (4 6)
— ‘2)1. ‘2 l '

 

(It. can be assumed. without changes in the results that U0 = 0) The solution of the
problem in the harmonic approxirnation is well known. The energy levels are given

by :

v

I
E” : hw‘I71+§).(~l.7)

In the harmonic potential the average values of kinetic and potential energies are the

same: < Ekm >=< EM >= F],,,,/2. Thus

 

p") [lec'2(;l' — .170).2 l , _
< — >:< >= 4...; <4 7 > —. 4.8
2/1. 2 21 l n( ) +2] ( )
From this we get:
‘ l h r 1 1
0f1c2.\1(T) 3< (l‘ ’51'vlz >: El< "(1) > +§l (4-9)

where abbreviation HQM stands for Harmonic Quantum Mechanics. The Bose-

93

*9

r...

l

5.4.aL.

.lll'dl

E

a

“‘4‘,

I

...a
«\g

Einstein function < 71(7’) > is given by:

l
< n(’1‘) >= n j . (41.10)
exp lfi] — l

'L

 

Expanding this at high temperatures (hw)/(ka) << 1 we get:

ka

-'1. l l
flu) ( )

 

tow-e

< 72(T) >2“
Thus, at high teinperatures:

< (.r — 1:0)“ >= —— = i ., (‘1-12)

 

If the spacing between the energy levels AE = free is much smaller than ka (i.e.
fw/ka << 1) the problem could be solverl classically. In the last case. the probalnlity

to find the particle at some coordinate .r is given by the Boltzman probability:

 

 

 

Pcatfl‘s T) = Z37,” exp l— [::;)]’ (413)
where
ZCM : / exp [—i:;)]d.r. (4.14)
So that:
< .1' >= /.L'PCM(;I:,T)(1.I:, (4.15)
and
< .172 >= /1‘2PC',\,(.1:,T)(1.13. (4.16)

94

.- -o

Thus:

a§.,,,(T) :< (:‘r— < :r >)2 >=< 3:2 > — <1? >2 (4.17)

could be found. In the harmonic case. (U(.r) = pw2(;r — .r0)2/2): < 1' >= .170,

 

27rka/ha'2 and U§,C,,,('F) = ka/pcu2, which is the same result as (4.12).

At higl’ier temperatures, when the amplitudes of atomic vibrations are larger, it
may not be appropriate to use the harmonic approximation for U(.r). Instead. it is
necessary to solve the problem for the original potential U(:r). Again, as before, if
the spacing between the energy levels is much bigger than ka the problem should be
solved quantum mechanically. In the opposite limit AE << 1.“th the classical solution
could be employed.

There is another potential L(z) that allows the exact QM solution. That is
the Morse potentialIG-l], which is widely used to model the properties of diatomic.

molecules:

UM(.r) : —(..r",,,,,, + Uo{l — exp [—a:r]}2 (4.18)

= I—UHnn + U0) + (.f,,[cxp [—20-1‘] — 2 exp —[(r.r]].

In this potential there is a finite number of the energy levels in the discrete spectrum.
The values of the energies are given by:

(if), 1

En : _Lfmm + L’fo — (jo 1'_ _— I + — 2~ 419
< > [ ml )1 . < >

2

where 72. runs over the positive integer values from zero to the maximum value at

which the term in the square brackets is still positive.

95

A>Uv >1.._Q:h.~

-2936

-2938

-2940

-2942

Energy (eV)

-2944

—2946

 

0.5 l 1.5 2 2.5 3
Distance r (A)

Figure 4.1: The blue curve shows the Morse potential for the N2 molecule. The red
curve shows the harmonic approximation to the Morse potential at small values of :r.
The blue and red horizontal lines are the energy levels for the Morse and harmonic
potentials respectively.

96

 

H"

.. ”l?“

Ts

The non-normalized wave functions are given by:

, 52 3
Ulnar) : 6X1) i—E—lg wn(£)s (4.20)
where
2V? (2,
6: —"—Iexp(—a.r), (4.21)
oh
and

V 2 (“in 1 '
s = [—n — (71+ —). (1.22.)
(if! 2

Finally u.'n(£) = F(—n,. ‘28 + 1,5) is the confluent hypergeoinetric. function that. can

be found by the series:

 

 

‘ , ~ _ Oi (f)((l+l) ;_'2_ 9‘
(itn+1)(a+‘2) 3_3 a(n+1)(o+2)(u+3) 5:
+ ‘;'(1v+1)(<r+2) 3! + 7(w'+1)t7+2)(*z'+55) 4! +
In our case, when a = —n. the confluent hypergeonietric function is a polynomial of

degree n.
If the energy levels and the wave functicms are known, the probe-ibility to find the

particle at temperature T at position I is given by:

 

 

 

1 "In!!! t) E
P “(1.3T = ,7 'd’n :1: “exp — :1, , 4.24
(2 ) ZQAI n; l ( )l [ [Cb] ] ( )
where
"m” En
r _ _ _. . r
ZQM — ”2:; exp[ MT]. (1.2a)

97

llif‘Il <’

n: '.
tumui

 

'Hr ,
i WY”-
'4. ‘9;

 

Then < .r > and < 3‘2 > could be found as:

< (103/(T) >= /;I‘PQJ\[(.1‘)((;I',

(4.26)
and
< .13).,”(7‘) >= /.172PQM(.1')d;r. (4.27)
Finally from (4.26.4.2?) we get:
055,0) =< (.r — .7)? >=< :lthU") > — < .2:QM(T) >2 (4.28)

In these calculations the continuous spectrum was ignored. Thus, being exact at

low ten'iperatures these QM calculations should fail at very high ten’iperatures.

For the harmonic potential. the sums (4.24, 4.25) could be explicitly calculated
[63] with the following result:

 

})HQM(J“ T) = (3 )% 0X1) [—01.2]. (”4.29)
7r
where
mu; he , .
a — TL tanh [kaT]. (4.30)

If the temperature is high enough, instead the Boltzman probability (4.13, 4.14,
4.17) with U (.r) = UM(;r) could be used to calculate 02. MC and MD simulations,

being classical, produce 0'2 that correspond to the Boltzman distribution at all tem-

peratures. Thus, their results are not valid at very low temperatures.

98

 

 

. _F xd\ i.! . l I a-
.. ~ I‘I-IIV Ii
. . - n .-

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 l 1 l 1 11L11 11l1111l1111J1 8
105 0 g: 10 :
E 51 :4
I Z d _
O I I I l IF IIIIIIIIIIIIIIIIIIr 0
L J 1 l 1 11l1111l1111L111ml1 8
10.31 i: 20 _'_
E E: —4
~_ E- —
A a _ " -
X 0 I I I I I — III III I I I I I Uri O
Vt:
a L l 1 l 11l4111l1111l1111l1 8
2910—: 2 i: 30 r
5 3 E: :‘4
2% a 5. —
“g 0 I T I I IF II IIII I I I I I I I II 0
Q-a 1 l 1 l 11lra11l 1 1 11l1111l1 8
103 3 i: 40 ;
E E": :4
0‘! I I I— I —O
1 1 l 1 11l1111l1111 l 1141 l1 8
2 I .. l-
105 7 5.: 50 :
E E: :4
0‘1 I I I l- "'"0
1.2
Distancer(A)

Figure 4.2: The squares of the wave functions of the Morse potential for the N2
molecule. The little figures on the left from the top to the bottom show normalized
wave functions for n : 0, 1, 2, 3, 7. The figures on the right from the top to the bottom
show normalized wave functions for n : 10, 20, 30, 40, 50. There are 61 energy levels
for the chosen values of parameters.

99

\'- lii

V

'0 U.
L‘ ‘53-:
.. . ”1
Hui. \1

‘.r. "Vl‘

49" A.1.
.

~-
ll..'lAI.
1" r'lII
“.,“H
II ‘
ill". 'I
\._A
\
.-.‘ l
"l 1 .

v
'V
;a», rm.

4.3. 1 Molecule N2

As an example. we consider the Morse potential UM(;I?) parameterized for the
1‘32 molecule. The blue curve on Fig.4.1 shows the corresponding Morse potential.
It was obtained by pI-n'forming calculations with the Gaussia1198 [(30, 61] program
that was approximately solving the Shrodinger equation for the electrons when the
distance between the Inicleus of N atoms was fixed. The equilibrium distance r0 :
1.0828 A and the equilibrium energy U,,,,,, = 29470188 eV were found in a separate
optimization run. In the same way. the energies EN of the individual N atoms were
calculated. Double of this energy is 2EN = —Um,-n + UO = —2936.6608 eV. This
information is sufficient to construct Morse potential (4.18). The red curve is the
harmonic approximation for the Morse potential at small values of :r = 'r — Tmm- The
blue horizontal lines are the quantum energy levels of the Morse potential. while the
red dashed lines are the energy levels of the corresponding harmonic potential. For
the used values of I;)ara.1neters, there are 61 energy levels in the Morse potential.

The squares of the normalized wave functions of the Morse potential for the N2
molecule are shown in Fig. 4.2. The wave finictions corresponding to n = 0. 1, 2, 3. 7
are on the left and for n. = 10. 20. 30. 40, 50 are on the right. Thus, for small 12. the
particle could be found only near the minimum of the potential r0. As n. increases,
due to the asymmetry of the potential, the probability density, as well as < r >
shifts to the right, and thus the particle spends more and more time away from the
minimum.

Figure 4.3 shows probabilities to find the particle at a. given distance in the Morse
potential and its harmonic approximation at different temperatures. For a diatomic
molecule these probabilities represent the pair distribution function. At very low
temperatures the CM results for Morse potential and its harmonic approximation
almost coincide. However, the QM results are different. In QM, the ground state

energies for the Morse and Harmonic potentials are separated from the bottom of the

100

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

60 I I l I I l I l I I I I I I4 AL I I 20
- 100K — _ 1000K _
50— ()\l\l\“1\c _— — I6
40— IIQMHannonh _ 5.3 -
a ELNI Harmon _ — — 12
30—1 Corrected — ~ ' -
_ ; Morse
a — —8
20— — i _
10‘ I l. —_ \ —4
a» - i - ‘ '
E 0 I I I I I I I I I I I I TI I I I 0
c3 1 l l 12 l 11 12
'8 l l l l l
E 15 — l l I I J 1 l l _ 1 I I I I I I I I I- 6
‘ 2000K : 10,000K t;
_ E f :‘5
10— e—j :—4
Z I € E3
5— —{ E2
: ' : -_ :—l
0 I I I I I I I I I I I - I I I I I I I I I II I I - 0
1 1.1 1.2 0.8 l 1.2 1.4
Distancer(A)

Figure 4.3: The l’air Distribution Functions or the probabilities to find the particle
at a given distance in the Morse potential or its harmonic approximation at temper—
atures 100K, 1000K, 2000K and 10.000K. The red and orange curves represent. QM
(4.24.425) and ("M (4.13.411) results for the Morse potential. The blue and green
curves represent QM and ( ‘M results for the harmonic approximation to the Morse
potential. The black curves represent C‘M results for the Morse potential corrected
by the convolution (4.34) with (73,,” obtained from the harmonic approximation to
the Morse potential.

101

'i
.
..
~71
\I"
A.

.L

.‘I

13‘1
.41.»
v

I
5 .

. . _ .
a . x. 1.. .2... . I I
1. I ‘ 1| A I I I Irut
”H1 w u .41.. . il a I '11“ O ..
1r u I1 v.. I 4; I .1 \I .v .IU 111.
I L 1'4 I I (II .IIIV IHU

potential by finite (slightly different) energy values. Because of it. the anharmonicity
of the Morse potential is reflected even in the ground state of the particle. This
leads to the slightly different PDFS that were obtained using QM on the Morse and
harmonic potentials. At temperature 10,000K the non—Gaussian line shape of the
PDFS obtained on the Morse potential are well pronounced. The depth of the l\=lorse
potential corresponds to z 100. 000 K

Figure 4.4 shows dependtmcies of 02(T) on temperature T that were calculated
in four different ways. The red curve shows the result of the QM calculation for the
Morse potential, where the smnmation in (424,425) were performed (Immerically)
over the 61 energy levels. At the temperatures plotted T < 10,000 K this result
is not any different from the result that could be obtained by the sunnnation over
the first 35 energy levels. Thus, it can be assumed that the range of integration in
(4.28) is (—oc,oc), while for the calculation the range (05,25) Awas used. The
orange curve shows the result. of CM calculations of 02(T) for the Morse potential
(4.13,4.14.4.17). There is a subtle point concerning the range of the integration in
(4.17). This question is discussed in Appendix B. The blue and green curves show
QM (4.9) and CM (4.12) results for the harmonic potential.

Note that 02.,‘,(T) calculated classically for the Morse and Harmonic potentials,
converges to zero as T —+ 0. In contrast, 032M(T) calculated by QM methods for
both potentials converge to a finite non-zero value as T ——+ 0. This is the effect called
zero—point motion. In other words, atoms in QM are not. motionless even at zero
temperature.

Classical MD or MC simulations performed on the .\~Iorse and Harmonic poten-
tials would lead to the orange and green curves respectively. Thus, classical methods
lead, especially at low temperature, to significantly incorrect results. In order to
obtain the correct results, it is necessary to use the QM approach.

Potentials that are used to model molecular motion for complicated molecules

102

 

$5 ...b .255 is;

 

0.006 — -
_ QM Morse :
_ CM Morse _
0'005 - — QM Harmonic _

CM Harmonic
------- CM Morse Corrected

 

 

 

 

3’3 0.004 — —
”b i i
.45
"U - r-
; 0.003 — +—
x _ -
g _ r
a. _ _
0.002 — —
-1 I-
0.001 — —
O I I I I I I I I I
0 2000 4000 6000 8000 10000

Temperature T (K)

Figure 4.4: The dependencies of (72 on temperature 'I' that were calculated in four
different ways. The red and orange curves show QM (4.24, 4.25, 4.28) and CM (4.13.
4.14, 4.17) results for the Morse potential. The blue and green curves show QM (4.0)
and CM (4.12) results for the harmonic approximation to the Morse potential. The
black curve shows the C M results on the Morse potential with correction ((1.3l) that.
comes from the harmonic approximation to the Morse potential.

103

 

D

‘ \
.L.

‘ .

VI;

 

.\x

AIL

[.1 A

R... I.-.

 

 

 

do not allow exact QM calculations. The harmonic a})proximatitm still could be used
to solve the problem at low temperatures. It can be employed if equilil'n‘imn configu-
rations, vibraticmal frequencies and vibrational eigenvectors are known. In this case.
it is possible to obtain both QM and CM solutions. But the harmonic approximation
is not valid in the most interesting regime when more complex molecular motions
appear due to the anharinonicity of the real potential. On the other hand, MD and
MC simulations allow to take into account the anharmonicity of the real potential.
Thus it would be valuable to find a way to correct the results of the CM calculation in
such a. way that at very low temperatures the corrected results would reproduce the
results of low temperature QM calculations, while at high temperatures they would
track the anharmonicity of the real potential.

Consider the following expression:
020‘) = 02"...1'1‘) + [vamp — 02mm], (4.31)

where under 03.,”(T) we mean the CM results on some potential U (:I') that could be
obtained from the Monte Carlo or l\Iolecular Dynamics simulations and thus sensitive
to the anharmonicity of the real potential. Under 0?,QM(T) we mean the QM solution
for the harmonic approximation to U (.r) in the vicinity of r0. Finally Uflcuule
stands for the CM solution for the harmonic approximation to U (1?). At very low
temperatures, 02(T) = 0'12,QM(T), since 02011le 2 0 for the both potentials. As the
temperature increases, see Fig.4.4, the difference of the two terms in square brackets
in (4.31) decreases and thus becomes: 02(T) 2 UEMU'). Thus, the form (4.31)
reproduces the correct behavior of the real 02(T). It follows from Fig.4.4 that results
that were obtained using CM on the Morse potential, and corrected by (4.31). the
corrections that comes from the harmonic approximation to the Morse potential,

reproduce the exact QM solution to the Morse potential in the whole range of the

104

“.1:.

'Ivvi“
\Lab“

\ .

VI

,.

'V
.11 ‘5

1A

.‘.'1‘.

F‘II
...as

lg’r.‘

VAL

‘0

;,.

i

... .‘.

371+

.4-

‘1‘.

reasonable temperatures rather well.
In general the line shape of the PDF for a pair of atoms obtained in MD or MC
simulations is not. Gaussian. Thus the correction (4.31) could not be applied directly

. , . .) .
1n order to correct (4.3) smce the (7,“!- are not. known. It is known that:

 

oo 1 I (7‘ _ I,I)2] 1 [ (1‘0 _ 7.I)2]1 I
—— ex .. ex ) — p-
.3, 27m? p 20f) 27.77% I 20:22 ( I
1 _ do 2
= exp [— 7 r ) I. (4.32)

 

 

Thus, the convolution (41.32) of one Gaussian with another increases the width of the
original Gaussian. If instead of one Gaussian another function that has some peak
in it would be used. then the convolution still would increase the width of the peak,
while the height would decrease.

In MC or MD simulations, the distance between a pair of atoms would oscillate
around some average value, at low enough temperatures. Thus, the PDF obtained
from MC or MD simulations would represent some distribution with a peak. Thus

the correction (4.31) still can be applied by the convohition:

1 (I _ .’./)2
G r = G, 7" —eX)———dr' . 4.33
( ) / III)( )W I 20.30” I )
where the (:07‘1'(::I:ti(m width is given by:
030”le = l”f{Q.iI(T) — UTIC‘AI(TII' ("I-3‘1)

In (4.33, 4.34) CMDU') stands for the PDF (I'listribution of distances) of the atomic
pair obtained from the MC or MD simulations and G(r) stands for the corrected
PDF. The weighted sum over the different pairs leads to the corrected PDF (4.3).
The black dashed curves on Fig.4.3 Show the CM results on the Morse potential

(413,414) corrected by the convolution (433,434).

105

k’.~ ‘AL 5

_ ‘.
:4 ‘v
.1: 1,1

I I
hit
I: ‘

4“At

1
h.

.. -.
\ 1
'rf‘T
c
,.

4.3.2 Correction to the heat capacity

The same philosophy that was used to correct 02(T) could be used to correct
the classically (.'alc‘éulated heat capacity. According to QM, the average value of the
total energy of an oscillator (a diatomic molecule) in a potential (Morse or Harmonic)

could be found as:

"7 mar

2 E,, epr—

1120

 

 

< EQ\[(T) >= 13%,] (4.35)

ZQM

where ZQM is given by (4.25).
In the classical approach, the average value of the total energy could be found

as the sum of the average values of kinetic and potential energies:
< Egg/(T) >=< AIT) > + < U(T) > . (4.36)
The value of the kinetic energy at temperature T is given by [63]:
, , . 1
< 11’ (1) >2 gka, (4.37)

while the average value of the potential energy could be found as [63]:

 

 

. 1 U(;IT) .
< U T >2 /U 1: ex) — r, (11:, 4.38
< ) 2.... t > I i all < >

where ZCM is given by (4.14). See Appendix (B) on the integration range.
If the average value of the total energy for a diatomic molecule is known as a

function of temperature, then the vibrational heat capacity could be found as:

 

C(T) = f . (4.39)

106

 

 

Again, as we did in (4.31), we can consider the correction of the CM results for
the Morse potential by using the results from the harmonic approximation to the

Morse potential:

E(T) = EMC‘i/(T) + IEHQMU) — EHC‘.\I(T)I- (4.40)

Figure 4.5 shows the dcpemlence of the total average energy for a particle in the
Morse and Harmonic potentials. The values of parameters are the sai'ne that were used
to obtain Fig.4.4. The red and orange curves were calculated via the QM (435,425)
and CM (4.36.4.37,-1.38) for the I\lorse potential. The blue (QM) and green (CM)
curves were calculated for the harmonic approximation to the Morse potential. The
fact that the corrected curve is very close to the exact QM solution obtained on the
Morse potential shows that this correction method also works well for the average
energy in a wide range of temperatures.

Note that in the limit of zero temperature there is a small difference between
the curves that represent QM solutions for the Morse and Harmonic potentials. This
difference is caused by the fact that there is a very small difference between the first
few energy levels of the Morse and Harmonic potentials, i.e. the energy of the ground
state of the Morse potential is not exactly the energy of the ground state of the
Harmonic potential.

The curves that represent CM solutions for the Morse and Harmonic potentials
coincide exactly in the limit of zero temperature as they should. As temperature
increases the CM Morse curve deviates from the CM Harmonic curve. In the limit of
very high temperature this classical solution for the Morse potential should coincide
with the QM solution (exact) for the Morse potential.

Thus, the corrected curve (black dots) coincides with the QM harmonic curve

at very low temperatures (limit of zero temperature), while at high temperatures it

107

 

_1. “'2‘;

 

 

10000 a _,__
QM Morse
CM Morse /
— QM Harmonic -'

 

 

 

 

”To - CM Harmonic
<4 8000 ‘ ------- CM Morse Corrected ,='_.* —
35 — -
5
Cl)
5 6000— _
C3
33
Q .1 >—
E
B
.E __ __
>3 4000
CD
5
c - 1—
Ln

2000— L

O I I I I I I I 17 T—
0 2000 4000 6000 8000 l 0000

Temperature T (K)

Figure 4.5: Average energies for a particle in the Morse potential and its harmonic
approximation as a function of temperature. The red and orange curves represent
the QM and CM solution for the Morse potential. The blue and green (dashed)
curves show the QM and CM solutions for the harmonic approximation to the Morse
potential. The black dashed curve shows the CM solution for the Morse potential
with correction (4.31) that comes from the harmonic approximation to the Morse
potential. It is assumed that U“ : 0.

108

 

coincides with the QM solution for the Morse potential. Thus the corrected curve, as
temperature increases. performs transition from the blue (QM Harmonic curve with
the smaller values of energy at the same temperature compare to the red QM Morse
curve) to the red curve (larger energy values). Thus, the (.lerivative (heat capacity) of
the corrected curve should be slightly bigger than the derivative from the blue (QM
Morse) curve in some range of temperatures (between 1000K and 5000K in Fig.4.6).

The heat capacity can be obtained by the differentiation (4.39) of the energy
curves. Figure 4.6 shows how the heat capacities found by the differentiation of the
energy curves on Fig-4.5 depend on temperature. Since in the harmonic potential,
the average values of the kinetic and potential energies are the same and due to (4.37)
the classical heat capacity for the particle in the harmonic potential is just A}, or unity
in units used on Fig.4.6. The behavior of the QM solution for the harmonic potential
is also well known [63]. Thus, we can see that our correction method again leads
to a very good agreement. between the approximate (corrected) and the exact. (QM
Morse) results.

Note that there is a small difference between the corrected curve (black dots)
and the curve that. represents the exact solution for the Morse potential in the range
of temperature between 1000 K and 5000 K. This difference was already discussed
above.

In QM approach the (,lept.‘ndence of the heat capacity on temperature can be
found not only by the differentiation of the energy curves but. also directly from the
knowledge of the energy levels. Differentiation of (4.35) with respect to temperature

leads to:

l
< CQM(T) >= T[< E5M(T) > — < EQM(T) >2], (111)

109

 

t
."

 

———-— QM Morse
CM Morse
— QM Harmonic ”

----- CM Harmonic '
------- CM Morse Corrected ~

 

0.4 —

Heat Capacitance (Dimensionless)

 

 

 

I I l I I I l T I
0 2000 4000 6000 8000 10000

Temperature T (K)

Figure 4.6: The dependencies of the heat capacities for the particle in the Morse po-
tential and its harmonic approximation on temperature. The red and orange curves
represent QM and CM solutions for the Morse potential. The blue and green (dashed)
curves show the QM and CM solutions for the harmonic approximation to the Morse
potential. The black dashed curve shows the CM solution for the Morse potential
with the correction (4.31) that comes from the harmonic approximation to the Morse
potential. All these curves were obtained by the differentiation (4.39) of the corre-
sponding energy curves on Fig.4.5.

110

 

where < EQM(T) > is given by (4.35) and < EQM(T) >‘2 is:

’1 m 11'

Z E,:exp[ —A:T]. (4.42)

7120

< [QM/(T

 

 

22: 11

These summations can be performed for the Morse and Harmonic potentials. More-
over, in case of Harmonic potential, when energy levels are given by (4.7): E,, =
132.1;[11 + 1/2], the sunnnation (4. 35) can be made in a closed form leading to a well

known result:

 

 

1
< E11Q,\[(T) >2 hw-‘(< 72(7) > 4:3), (4.43)
where, as before:
. . 1
< 11(1) >2 f (4.44)
ex1)[A"“T] — 1

is the Bose—Einstein function. Expression (4.44) was previously used in (4.8). Differ-

entiation of (4.43) with respect to temperature leads to:

d < EHQM(T)>f1w hm

— —k()(— )Cxpf— Ab T

< (T >2. 4.45
(H. ka —-l "4 ) f )

CHQi/(T) =

 

In the classical approach for the harmonic potential. since the value of the total
energy is given by MT, the classical harmonic heat capacity is just 11:1,. Thus, the
value of the COIIO( tion te11n [C '11QM(7 )— CHC 11] can be easily found.

In the classical approach for the general potential we have to differentiate ex-
pression (4.36) (using expressions 4.37 and 4.38). Thus, we get:

, 1 1 , , .
< C'C,\[(r[)> >Z Eli), +— l(< Ucz 111(1) > — < UC‘JtICI) >2), (4.46)

111

 

where < UCM(T) > is given by (4.38) and < UQM(T) >2 is:

 

 

1 1) U(;I7)
< U T >= , /U‘ 11' ex — , (II. 4.47
( ) ,CM ( ) pl k1. Fl ( 1

Last formulas can be employed for the Morse potential. Thus we write:
CVQM = CMCM + [Chen/(T) — CHC'M] (4.48)

Both approaches: numerical differentiation of the energy curves and calculation of
the heat capacity using formulas above should lead to the same result, as they do:
the differences between the curve obtained by the numerical differentiation of the
corrected energy curves and the curve obtained using the formulas above would be
invisible 011 the scale of Fig.4.6.

In case when there are several inrlerwndent (orthogonal) vibrational modes the
correction should be applied to each mode independently. This can be used to correct

the heat capacity for the 111anyatomic molecules that. have several vibrational modes.

4.4 Vibrations of large molecules

Let us assume initially that atoms in the molecule vibrate near their equilibrium
positions with amplitudes of vibrations that are much smaller than any interatomic
distances. Thus for the moment we ignore the possibility of the large-scale atomic
motion that can be caused by the flexibility of the molecule. In order to calculate the
PDF it is necessary to calculate the mean square deviations of the distances between

the atoms from their equilibrium values: 0,2 Let the coordinates of atoms 1' and j

j.

be F, = 1"? + '17,- and 'F- = I“? + 17, where 7’“? and 17‘? are the ec uilibrium coordinates
1 J J J l J l

of atoms 1' and j, while their instantaneous displacements are 17,- and 17,-. Let also:

thj=7z}*7:i,7‘~—T "7:9

,J J ,, and 11,-]- = ”11]- — 11,-. Then 1t IS easy to show, under the

112

 

assumption 110 << 7‘8, that:

05:<(1'i1—I‘?J-)2 >§< [7.317012 > (4.49)

_ E )200 1013 , . , . _

a .13

where r3 = 17::- / rj’J. The indexes (1. .3 stand for the Cartesian coordinates of the
vectors.

We consider here the vibrations of a. molecule that consists of N atoms with
masses 111,. The potential energy of the molecule (7(7) is a function of the coordinates
.1" of all the atoms. There are 3N comprments .1:, of the vector 7". Thus the index '1'
runs over all the atoms and over all their Cartesian coordinates.

Let us assume that the atoms vibrate near their equilibrium positions .170 and

that their instantaneous coordinates are .17 = .70 + 17. If vector 17 is small enough the

potential energy could be expanded near the equilibrium (fa):

1
(17(17) 31 U0 + 3 E 0011,11, (4.50)
ii
where
I),~ : —— 4.51
J 0.171815 ( d )

is a real symmetric 111atrix. Thus for the Hamiltonian we have:
2
P1 1
H = —— + — E Dw’u'uu. 4.52)
, 2"” 2 .. U 1 J (
1 1}
It is convenient to introduce new variables

P1 = "11771, 'Uvz' = Cir/W (4-53)

113

 

that will transform (4.52) into:

71'? 1 ~
H = Z 7 + i Z D..q.q.. (4.54)
U

i
Where Di} = [)ij/1 /111,111J~.
Since B!) is a. real and synnnetric matrix it could be brought to the diagonal
form by the linear t1'a11sfo1'1nat101is of q,- :
- — A 155
(ll _ 61 QA: ( ' )
A

where 5" could be chosen as real, orthogonal and normalized:

2614614, 2 (SAN- (4.56)

In this case the transformation:

71, = Z 831), (4.57)
A

will also bring to the diagonal form the first sum in (4.54). Thus we can rewrite

Hamiltonian as:

1 .
H : Z E[103+ wiei]. (4.58)
A
The second quantization:
1 f1. + 73 r— +
Q) =5 LAO—)0” ‘f‘ CIA), ,PA :5 hLJJ)‘(C1/\ — CIA), (44.59)

114

 

with [R Q] = 171 so that. (1:11),, — axe: = (5M: transforms (4.58) into:
H = 2 mafia, + 1). (1.00)
A 2 i

as it should. Thus < 'drlfll'u'i >2 2, 1101.401 + 1/2).

From (453,4 55 4.59) follows that:

71
. /\ + ,
‘ '1' : ,3~ . . 4.61
11 EA 6, l( 2'111.,w1(a’\ + (1)) ( )

Thus, separating the atomic and Cartesian coordinates, we get:

 

 

 

<1,(11 11 II, >= Zi- (11 +—) 60” e)», (462)
m J3 2.1.3 A ‘2 run/111}. '

It is easy to show, using (4.62), that (4.50) transforms into:

Ui=Z—h—1n.+ —1{W———)—w12 (463)

,J A 24A (fit—1 m , .

where 11), is the Bose-Einstein function (4.10) that at high temperatures can be sub-

stituted by (4.11). It is easy to see that. (4.63) in case of the diatomic molecules

(m1 = 1112 = 111.) reduces to ( 4.9), if to take into account that 11 = 111/2 and
(517:0) — (5172-11) = 2-

MD or MC simulations at. low temperatures, when the anharmonicity of the real
potential can be ignored, would lead to (4.63, 4.11), while the correct result is (4.63,
4.10). It is clear from (4.7, 4.9, 4.10, 4.11, 4.12) and from (4.60, 4.63) that the formulas
(4.31, 4.40) used to correct the CM results in case of the diatomic molecule can be

modified in an obvious way in order to correct the CM results that were obtained on

115

 

Tm

a Iiiaiiy-atomic molecule, i.e.:

E(T) = (4.64)

5011le + ZlEHQ1\I(A3T) - £7110th, TN-
A

The width of the peaks can be corrected via

02(7‘) : 02'1” (T) + 0301'r(T)' (465)
where
030.17“) = Zlotmw T) — 02.1.1111] (4.66)
A

and the index /\ runs over different vibrational modes.
From the result. of MD or MC calculation only the total (not for the particular

vibrational mode) average values of the ECM(T) and 037M(T) could be found.

4.5 TIN KER package for molecular modelling and
the CSHM molecule

In order to demonstrate the role of the quantum corrections (4.64, 4.65) for a more
complicated molecule, it is necessary to find the eigenfrequencies and eigenvectors of
molecular vibrations as well as the results of the classical MD or MC simulations.
There are many different potentials and programs developed for the different kinds
of molecules. One of the most developed potentials are those developed for hydrocar-
bons. One of them is the so called “MM3” potential [65, 66, 67, 68, 69, 70, 71, 72].
One of the programs (a package that includes several different programs) that can

use this potential to optimize the molecular structure, calculate the eigenfrequencies

116

 

 

Om

 

 

 

 

 

‘23
‘9
N‘ O
O

 

 

 

 

Figure 4.7: The sketch of the (76H14 molecule. The carbon atoms are shown as black
filled circles. The hydrogen atoms are shown as open circles. The numeration of the
atoms coincides with those used in the Fig.4.9 and Fig.4.10,4.1l,4.12

(1.11),) and eigenvectors (5*) in a particular equilibrium configuration and to run MD

at given temperature is “TINKER” [73, 74, 75, 76, 77, 78, 79]. As an output of MD
simulations (for a single molecule) at given temperature, TINKER provides the coor-
dinates of the atoms. the total energy, the kinetic and potential energies as functions
of time.

We use the (761714 - molecule as a test. example. Figure 4.7 shows the sketch of
the CfiH 14 molecule, while the F ig.4.8 shows the geometry of the C'GH 1.4 molecule in
Long conformation. At room temperature the 05H 14 molecule is relatively flexible,
as we will see. At normal pressure these molecules form a non-toxic and not very
flammable liquid and thus its PDF can be easily measured.

At low temperatures the C611” molecule can be in two different conformations,
as follows from MD simulations and structures optimizations with TIN KER, which
we call Long and Short. The distances between some atoms in these molecular confor-

mations are shown in the table 4.1. The ground state tension energies (obtained from

117

 

Figure 4.8: The geometry of the CGH 14 molecule in Long conformation.

 

Pair ofAtomsi—j 1—4 1-—5 1—6
Distance between i and j (Long). (A) 3.94 4.56 5.07
Distance between 1 and j (Short), (A) 3.17 4.54 3.81

 

 

 

 

 

 

 

 

Table 4.1: The distances between some atoms in Long and Short conformations.

TINKER) of the molecule in Long and Short conformations are slightly different:

EL = 6.3479 (kcal/mol), (4.67)

E5 = 7.1293 (heal/moi),

where L stands for the Long and S for the Short. The difference in the energies

between these two conformations per molecule in temperature units is:

AE = E; — EL z 393 (K). (4.68)

118

4.6 Calculation of PDF for CGHM molecule

4.6. 1 Low temperatures

The blue curves on F ig.4.9 show the probabilities for finding the pairs of atoms
1-2, 1-3, 1-4. 1-5, 1-6 at a given distance obtained from MD simulations at temper-
atures 150 K and 500 K. In the beginning of MD runs, the molecule was always in
an equilibrium Long conformation. At temperature 150 K the molecule remains in
Long conformations during the MD run. At temperature 500K the molecule con-
tinuously switches between different conformations (it is flexible). The time step in
MD simulations was A! = 1 femtoseconds. The coordinates of the atoms were saved
after every 1000 MD time steps. The blue curves represent the distributions obtained
from 10.000 different molecular configurations. Thus the total run time was 10,000
picoseconds. The red curves were obtained by the convolution (433,434) of the blue
curves with the Gaussian whose correction width am“. that corresponds to the given
pair of atoms for the molecule is in Long conformation.

At low temperature, when the molecule is frozen in the Long or Short configura-
tion, atoms only vibrate slightly near their equilibrium positions. As the temperature
increases, there can occur transitions between Long and Short. configurations. This
behavior could be seen from the change in the peak shape for the pairs of atoms 1-4,
1-5, 1-6 when temperature changed from 150 K to 500 K. At 150 K those peaks are
relatively narrow and symmetric, while at 500 K the peaks become much broader and
their shape indicates that. there are two conftn'mations in which the molecule spends
most of its time.

If the molecule during the MD run remains all the time in Long or Short con-
formation it is clear what set of the correction om, should be used. If the molecule
changes its conformation during the MD run, it becomes unclear what own set should

be used: from Long or from Short conformation. We will discuss this issue later.

119

'1‘ = 150 (K) ’12—— 500 (K)

 

 

 

 

 

 

 

 

 

 

20_ I I I _, ‘20— I I I _4
:1-2 : :1-2 :
10_— — 10— {‘1 —-
_ [1L _ : ‘3‘. :
ol— 1 mjl l — Oi— 1 Jl Ll T
1 1.5 2 1 1.5 2
10 I . 10 . I .

 

 

 

 

Ur
IIIIIIIII
_ \
LJJllllll

kl}
IIIIIIIT
lllllllll

 

 

1]
o

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

: 2 25 3 2 2.5 3
o 10_ j 1 I ~ 10_ I I I _
b 51-4 II E E1-4 3
5E 5;— : 5r —:
2% : II : g A :
'8 0— 1 .'1_- o Dawes/1L-
a... 3 4 3 4
a4
»- I I I I ..4 I— I I I I d
541-5 II— 5—1-5 —
Oi- l l 1 Jlk ,— L-Mr l —
3 4 5 O3 4 5
_ I I I I _ I I I
5—1-6 5_1-

-I
6 _..
—4
—q
I—

I
lLllll
I

III
III

 

 

 

 

 

 

J ‘\ r _*M‘
0 1 l l l 0 l

5 3 4 5
Distance r (A)

L»
.1:.

Figure 4.9: The pr(_)bability for finding a. particular pair of atoms at a given distance
at temperatures 150K on the left and at 500K on the right. The blue curves were ob-
tained from the classical MD trajectories. The red curves result from the convolution
of the blue curves with the corresImmling correction Gaussian. Every pair of atoms
has its own (72 (T).

' ( (11‘ I‘

 

 

I
E
1
I

[ll

lriil

 

 

 

 

   

_i.
1-$
.-
G
v
'
D
D

L
F

 

 

‘ 111 III III‘III
l. I I I!
q

 

 

 

 

 

 

 

 

 

Average Distance (A)

 

 

 

A

llllllLJllll

*7

 

 

 

 

 

rTIIIIIjIfiT

l l l L i l l 1

i I
500 1000 0 500 1000

 

 

0

Temperature T (K)

Figure 4.10: The dependencies of the average distances between some pairs of atoms
on temperature obtained from MD simulations. In the beginning of every MD run
(at a particular temperature) the molecule was always in equilibrium Long conforma-
tion. At low temperatures these distances only slightly depend on temperature, as
molecule remains in the Long conformation. As temperature increases instantaneous
distance between a pair of atoms can have significantly different values that occur
when the molecule can be in different. conformations (or identical atoms interchange
their positions). Thus, as temperature increases, there occur sharp changes in the
average distances for some pairs of atoms.

121

 

 

 

 

 

 

 

 

 

   

 

O 500 1000 O 500 1000

 

0.0] I I T I Vi T I I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_ 1 I I r I I I 1 £1! ,1
N" I 1' 4 I I 1’ 4 ‘ .. . _ 01
_ a ._ a _
6 - — - -
~ — - : :
b O 005 M '- 4‘} _l
r r r- m 0.05
_ / ’- 9, -
_ .///’ I *- f‘ -
_ .4 _ 'd ;.-—--':
1 1 l 1 1 1 l 1 “iii—”T’1 l 1 1 1 1
0
0 100 200 O 500 100%
0.01 I I I I 1 1 r _ r 1 1 1 I r 1 1 1 4
h 1, 5 ‘ — 1, 5 e 4
l. ._ «c a _l
:7—//: :— 9 —-‘ 0.2
0.005 —- /— : .° 2
: / : :— :90 _: 0.1
/ ‘ : °« :
0 1 1 1 l 1 1 1 (W 0
0 100 200 O 500 1000

Temperature T (K)

Figure 4.11: The dependencies of (12 on temperature T for some pairs of atoms.
The orange triangles show the results obtained from MD simulations, the blue and
green curves show the QM and CM results obtained from the eigenfrequencies and
eigenvectors of molecular vibrations. The red circles show the MD results corrected
by adding to the orange triangles the difference between the blue and green curves.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

A )— .-
(:l< — . a ' — 0 8
V _ - ‘ -
Nb 0.04 M _ Effie _ 0.6
— — 0.4
0.02 — — _ 9‘ _
_ _ — 5 — 0.2
I I I I I I- i I I I I I
0
0 100 0 500 10%
0.02 I I I I I I I I I I I I I I g 0
: 8'” 6 c : —- 8’11000 0 ° 0 :3 jiOI
0.01 —- — - -
_ — — . — 0.05
0 1 III-’1 1 1 1 1 1 0
0 100 200 0 500 1000

Temperature T (K)

Figure 4.12: The dependencies of (12 on temperature T for some pairs of atoms.
The orange triangles show the results obtained from MD simulations. the blue and
green curves show the QM and CM results obtained from the eigenfrequencies and
eigenvectors of molecular vibrations. The red circles show the MD results corrected
by adding to the orange triangles the difference between the blue and green curves.

123

The orange triangles on the Fig.4.11,4.12 show how the 0‘5.“ obtained in MD
runs depend on temI’)erature (See also F ig.4.10). The blue and green curves show
02(T) that were obtained in the harmonic approximation from the eigenfrequeiicies
and eigeii\-'ectors of molecular vibrations in the frame of QM (1.63.111) and CM
(4.63.410) respectively. The red circles were obtained by applying the correction
(4.65) to the MD results. Thus at very low temperatures the results of MD simulations
agree with the CM results ol.)tained in the harmonic approximation, as it should be.
However the correct results at these temperatures are given by the QM results on
the harmonic potential (the blue curve). Thus our correction. if applied, shifts the
orange triangles into the red circles that fit QM results at low tci'nI')eratures very
well. As temperature increases the anharmonicity of the potential usually makes
potential softer, thus increasing 02(7') compared with the harmonic case. The sharp
increase in the 03(7‘) with increase of temperature that occurs for some pairs of atoms
corresponds to the appearance of flexibility. Thus the molecule becomes flexible
between 200 and 300 K. It is easy to see, from the numbers of the atoms. that it is
the dihedral angles that become flexible.

If the molecular trajectory is known the PDF (4.3) could be found if the distri-
butions of lengths (the blue curves on F ig.4.9) for every pair of atoms will be used
instead of the Gaussians in (43). The correction could be made by the convolution
of the distribution of distances for every pair of atoms with the gaussains (433,434).
The brown curves on the Fig.4.13 show the PDFS obtained in MD simulations when
the molecule remains in Long conforn‘iation before the correction, while the blue
curves show the corrected PDF s. Thus at low temperatures the correction signifi-
cantly changes the shape of the PDF. The role of the correction remains important.
at small distances at all temperatures, while at large distances its role decreases as
temperature increases because the widths of peaks that correspmid to the atoms that

are well separated become significantly larger in MD simulations.

124

2000 llllIlllllIlllIlllllllllIllllllI1J_IJ_llllllll 1111111||l1111

 

50 K
1000

lllllllllIlllllllll
IIIIIIIIIIIIIIIIIII

 

 

O IIIIIIIIIIIIIIIIIIIlllllIIIIIi-III'IIIIIIIIIIIIIIIIIIIIIIIIIIIIII

 

2000 IlllIlllllllllIlllllllllIlllllllllIlllllllllIlIlllllllIlllI

‘ l

100 K
1000

 

IIIIIIIIIIIIIIIIIII

J_Llllllllllllllllll

 

0 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII—IIIIIII

 

2000 llllIllLllllllIl11IlllllIlIlliLlllIIllllllllIlLIlllIllIllll

Pair Distribution Function g(r)

1 200 K
1000 ‘

IllllllLJIlllllllll
IIIIIIIIIIIIIIIIIII

 

 

 

 

' A .3 A
0 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIllllllll

1 2 3 4 5 6
Distance r (A)

Figure 4.13: PDF for the (7.317)., molecule in Long conformation at low temperatures.
The brown curves show the results of MD simulations before the correction is applied.
The blue curves show the results of MD simulations corrected by convolution (4.33,
4.34).

It is not. known in general how many molecules are in a particular conformation.
The relative contributions of different conformations to the total PDF depend on the
shape of the free energy surface. Since the line shape of the free energy is not known
the simplest way to combine the results obtained on different conformations is to
assign the Bolzman weights to the PDFs in different confornrations. Thus the weight
for the molecule in Long conformation:

E L
hibT

 

'11),,('[‘) or exp [— ]. (4.69)

Thus for the ('~Il . molecule at low tem )eratures the total PDF can be calculated
0 11
as a linear combination of the PDFs from Long and Short conformations with the

correspcmding Boltzman weights.

G(I‘, T) = ’I_1.1'L(rll)(fL(l‘, T) ‘l‘ U-‘§;(T)Gs(l‘, 11), (4.70)

where CL(1‘, T) and (15(11'1’) are the corrected PDFS from Long and Short conforma-

tions. The Boltzman weights 117(1) and 1113(T) are given by:

 

 

1
111(7' = ,. . (‘1-71)
1 + exp [—35]
and
A}?
exp "'77: .
11‘5(7l) [ MI] (4.72)

 

—1+epr—%]i

so that: 111L(T) + 015(7) 2 1. Table 4.2 shows 111L(T) and 1115(T) for temperatm'es
50, 100, 150, 200 (K).
The corrected PDFS for the C61! 1.; molecule in Long and Short conformations

at different temperatures are shown in Fig.4.l4 as the blue and green curves respec-

126

 

15(K) lllllllllllllllllllllllllllllllllllllllllllllIIIIllllllllll
— I

— .l

1":

1000-5 1 l 50 K

-0"

1
t
I
1
1
I

1
I
l
_1 1" '1
I
1

I
' '(\\ (r‘
I \ g I
' I I I O ' I
—I I. i ' I. \I I
1‘ I ‘ 1' ‘\ ‘ - I.
l ' \ / \ ' \ I
—I " I l I. I ‘J T ‘\ I \ I, \
I \ I I \ I, -\ I s 0’ ‘\ I \
__‘ 1 t 1 \fi '1 s“ '1 \ ,1 x‘
o W—

I
I
0
I
t
I
I
I
I
u
1
t
t
|
I
I
1

I II I II II I II II I

 

 

 

 

 

 

 

 

 

r-x , ~, ’1 ‘. ____ ,1" ..--. ---m...
b O ITI‘I'TIITIIIIII 111111111 111111111 111111111 1111111771177
00 1 1 1 1 1 1
C 1111l111111111l111111111l111111111l111111111l111111111l1111
.9 1500 _, 1. _
‘5 1 1': _
E; ‘ 11 100K -
LL. 1000'“ t _
c: 4 5]
£3 : 1'1 :
:1 ._ 11 _
_D 500 d A ,, A _
o‘_‘ 1‘ I! 1 1
b 1 /\ 1!. /\ ”
052 11 l‘. 1" l\ ’[l \J / I _
‘ / \1‘ , ‘
Q 0 1IIT’IHI‘TIIITTIHHTIIHI11IIIIIHIIITIIIIIIIIIrIrHHI—TTT?
H
”-1
CG 111111111111111111111111l111111111l111111111l111111111l1111
a, 1500 _ 1. _
1 II -
1 ’1 1 K i
1000—_- 1I 50 ,—
—1 l _
1
—1 11 .—
_ ‘i ._
500— I1 ~—
‘ A l 10 0 I
—1 I l\ 1—
4 \ l / _
0 1 L] -

 

 

IIIIIITIITTITWIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIITIIIIIIII

l 2 3 4 5 6
Distance r (A)

Figure 4.14: The blue and green solid curves show total corrected PDFs for the
molecule in Long and Short conformations. The red dashed curve shows Boltzman
combination (4.70) of PDFS from Long and Short conformations. At all temperatures
when the molecules remains in one or another conformation the combined PDF ba-
sically coincides with PDF in Long conformation due to the values of the Boltzman
weights.

 

Temperature T (K) 50 100 150 200
uiL(T) 0.9996 0.98 0.93 0.88
'u.'S(T) 0.0004 0.02 0.07 0.12

 

 

 

 

 

 

 

 

 

Table 4.2: Boltzman weights for CGH 1.; molecule in Long and Short conformations at
different temperatures.

tivelx. The red curves in Fig/1.14 show PDFS calculated according to (4.70) with the

I
\J

Boltzman weights taken from the table 4.2.

4.6.2 High temperatures

In order to use the Boltzman combination of the PDFS in different conforn‘iations
to obtain the total combined PDF, the molecule should remain in one or the other
conformation during the MD run, or it is necessary to separate those times (parts of
the molecular trajectory) when the molecule is in one or in the other conformation.
It is necessary to do this because every conformation has its own set of crew... that was
obtained using the eigenfrequencies and eigenvectors of the molecular vibrations in
the corresponding conformation.

In practice. however. it is impossible to say when the molecule is in one or in
the other conformation. For example even when the molecule remains in Long or
Short conformation the positions of atoms 8 and 9 can interchange during the MD
run. When it happens the molecule is already in a. different conformation since atoms
8 and 9 have different. sets of own. with respect to the other atoms. Thus although
the geometry of the molecule did not. change during the interchange of positions of
atoms 8 and 9, the sets of own should be modified. It seems to be an impossible
task to track all such little changes. Another reason that does not allow separation
of different conformations at high temperatures comes from the observation that
the changes in the distances between some pairs of atoms, due to the vibrations,

become comparable with the changes in the distances caused by the changes in the

128

conformations. Tlms basically it becomes impossible to say when the molecule is in
one or in the other conformation. Thus, at. high temperatures, the approach that
con'ibines PDFs in different conformations with the corresponding Boltzman weights
becomes unacceptable.

Instead another approach could be used. In order to obtain the corrected PDF
in Long conformation it is necessary to convolute the PDF obtained from the MD
trajectory in Long conformation with the set of am,” that. was also obtained on the
Long conformation. Assume that. instead of convolution performed with the set of
aw" from the Long conformation we will use the set of (7mm. from the Short confor-
mation. How different would the corrected PDFS be when obtained with these two
different sets of am" from the same .\‘ID trajectory?

The brown curves on Fig.1.15 show PDFS obtained from the MD trajectory at
temperature 50 K when the molecule remains in the Long (top) or in the Short (bot-
tom) conformation depending on initial configuration. The blue solid curve originates
from the convolutions of the M D results from the Long/Short conformations with the
correct sets of 0am from the Long/Short conformations. The red dashed curves were
obtained by the convolution of the brown curves with the incorrect sets of am": the
MD results from the Long/512,071 conformations were convoluted with the am". from
the Short/Long conformations.

Thus, although there are very significant differences between the non-convoluted
and convoluted (corrected) PDFS the differences between the PDFS obtained by the
convolutions with different sets of 0mm. are rather small.

The small differences between the MD results corrected with the different sets of
0am. could be easily understood. Our correction is the most significant (large 060".)
for those pairs of atoms that are tightly connected, i.e. by a bond or by an angle. Thus
the correction is significant. for those atoms that are close to each other. The mutual

vibrations of such pairs (distance distribution) are slightly affected by the general

129

 

 

 

 

 

 

1500 LllllllJLLLllllLLLlllllIlllllllllllllllilllllllllllllllllil

t Long conformation atSOK -

t _

1000~ a

c _

A l— _
i—

v _.
00

c: 500— _

.9. — -' —

H 1

g a _

u. i i

’2 0 lillllllllllllllllllllllllllllTIITlTlll‘lllllllkfflll‘llll
Q—D

'8 1500 lllllllllllllllllllllllllIlllllllllllllllllllllllllllllllll

E , Short conformation at 50 K —

.: f -

a _ _

‘3‘ 1000— -

500— “ —

t *1 "

’- c '5. T

f _

. J .. v P.
0 Illlfllllllllllllllllllllllllll[TTTTTTTIIITTIIIIllllllll
1 2 3 4 5 6
Distance 1' (A)

Figure 4.15: The brown curves on both figures show MD results obtained on the Long
(top) and Short (bottom) conformations. The figure on top also shows the corrected
l’Dli‘s obtained by the correction of M 1) results from the Long conformation with the
set. of am" from the Long (the blue solid curve) and Short. (the red dashed curve)
conformations. The bottom figure also shows the PDFS obtained by the correction of
MD results from the Short conformation with the sets of am" from Short (blue solid
curve) and Long (red dashed curve) conformations.

130

conformation of the molecule-they are almost the same for all of them. On the other
hand, if atoms in a pair are far away from each other then the corresponding own. are
small in general, independently from the molecular configuration. Thus the correction
is not that important for those pairs in principle. That means that in order to apply
the correction we can convolute the PDF obtained from the MD simulations with a
set. of 0am from any particular configuration. At higher temperatures the differences
between the curves obtained by convolutions with different sets of oc0.r,.( T) will become
even smaller since 0“,”.(T) decrease with increase of temperature (see F ig.4.4,??).

Figure 4.16 show PDFS calculated for the C6H14 molecule at different tempera-
tures. The PDFS at temperatures 50 (K) and 150 (K) were obtained by combining
the PDFS from the Long and Short conformations with the corresponding Boltzman
weights. The PDF curves at temperatures 298(K), 500(K) and 800(K) were obtained
by the convolution of PDFS from the MD simulations with the Gaussians with the
correction widths from the Long conformations. It follows from the figure that the
PDFS in the region of r < 3 (A) only slightly change in the considered interval of
temperatures. The inset. that shows the region of 3 < r < 6.5(A) on a bigger scale
shows how the flexibility of the molecule develops as temperature increases.

The region 3 < r < 6.5 (A) is the most interesting if the large scale molecular
motion is under consideration. In this region our correction is not very significant
for a single molecule. However, in order to compare the results of calculations with
the experimental results. it is necessary to calculate the combined PDF from many
different molecules. The PDF from the different molecules can overlap in the region
r < 3 (A). Thus, in order to extract accurate information about the mutual orientation
of different molecules. (that in addition can put some constraints on a single molecular
motion), it will be necessary to accurately calculate the PDF at small distances, i.e.

in the region where our correction is significant.

131

1500 llllllllllllllllllllllllllllllllllllllllllllllllllllllllll
TlrrIITT—rrTlijldlITIIIWTTIT‘fIIITjI

 

 

400“ — 50K —
_ - 150K

l — -- 298l( _
- l 500K _
-" 800K

I

1000

 

 

 

 

\

U1

o

o
l
I

 

 

 

 

Pair Distribution Function g(r)

0 lllllllllllllllllllllllllllllllll

3

 

 

 

 

 

 

 

0 llllllllllIlllllllllllllIllllllIITIIIIIIIIIIIIIIIIIIIIIIII

l 2 3 4 5 6
Distance r (A)

Figure 4.16: Corrected PDFS for the (.'.; H14 molecule at different temperatures. The
inset shows the region betwcen 3 A and 6.5 A on a bigger scale.

I32

4.6.3 Pair Distribution Function from Molecular Vibrations

and Partial PDFS.

In order to understand how PDF of the molecule depends on the molecules
conformation, it useful to consider partial PDFS with respect to the carbon atom
number 1. It. is also informative to calculate partial PDFS from eigenfrequeiicies and
eigenvectors of molecular vibrations in Long and Short conformations and compare
those PDFs between each other and with partial PDFs obtained using results of MD
simulations with our (.'(_)rrection method. Figures 4.17, 4.18, 4.19, 4.20 and 4.21 show

few of such PDFS.

133

 

 

 

 

 

 

 

 

 

 

300 I l I I
t l —
250 — l 298 K ‘
I -
200 — 1‘ —
c I
.g - I. -
o 150 — ‘ —
C II“
E ' l . II ‘
C: 100 '— ( l I") ,' l _
_O _ l I t I I _
5 ‘ l I W.
5 50 _ II | _
r9 l l \I I
E i- l l ’ I ‘\
‘52 1 AL I \‘
Q 0 I I I
.‘.: 0 l 2 5 6
“5 l I I
an 300
t 250 —
a“: _ 298 K ,
a 200 — —
l" (- _
U 150 e _
V _ _
100 — —
50 — _
0 I I I I I
0 1 2 3 4 5 6

Distance r (A)

Figure 4.17: The panel on top shows PDFS with respect to carbon atom number 1
created by the other carbon atoms only at 298K. Blue and green curves were obtained
from molecule in Long and Short conformations correspondingly. The widths of peaks
were calculated in QM approach using eigenfrequencies and eigenvectors of molecular
vibrations. As distance increases there are peaks that correspond to the carbon atoms
number '2, 3, 4. 5 and 6. The red curve is a linear combination of blue and green
curves with corresponding Boltzman weights. The red curve in the bottom panel is
the same as the red curve in the top panel. The grey curve, that also shows PDF with
respect to the carbon atom number 1 due to the other carbon atoms, was obtained
from M D simulations. The black curve was obtained from the grey curve using our
correct ion scheme.

134

 

 

 

 

 

 

300 I I I L I
250 — ll Long
._ ii
.'I

c 200 — I': 298 K a
.9 _ r: _
LL. _ ‘ _
.5 100 — —
E — _
'r: 50 — _
D
.: 0
a“: 0 7
:5 300
g _

9.. 250 -
A .. _
2 200 —

I _ _.
U" 150 — —
V _ .1

100 — —
50 — —

0
0 7

 

Distance r (A)

Figure 4.18: The top panel shows partial PDFS with respect. to the carbon atom
number 1 calculated for the molecule in [long conformation from eigenfrequencies
and eigenvectors of molecular vibrations in the QM approach. The blue curve shows
partial PDF due to the other carbon atoms only (hyrh‘ogen atoms are ignored). The
red curve is clue to all the other atoms (carbons and hydrogens). Thus. the difference
between two curves is due to the hydrogen atoms. The bottom panel shows the same
curves as the top panel. but for the molecule in Short conformation.

135

 

 

 

 

 

 

 

300 I I I I I I
250 — . 298 K _
s: — l‘ —
2 200 j A
*5 ‘ 1 II I ‘
“=3 - l I II -
L“ — I n N 't I’ I ‘
.5 100 - n A II H —
E _ [I' ‘ y [I \\ l E \ 'll.‘ I‘" __
' 50 — i I l ‘ d\ If _
g _ l l l I /\ V , \tl _
Q 0 Jl “V l j \_ I, \ \‘I\
_I_. I I I I I I
<3 0 1 2 3 5 6 7
E 300 I I I I I I
c“ T
g I. a
250 —
i _ 298 K
2‘ 200 ~ _
' .. _
U— 150 ~ —
V I— d
100 — —
50 — —
0 I I I I I I
0 l 2 3 4 5 6 7
Distance r (A)

Figure 4.1.9: The top panel shows partial PDFS with respect. to the carbon atom
number 1 created by all other atoms (carbons and hydrogens). The blue and green
curves were calculated from eigenfrequencies and eigenvectors of molecular vibrations
in Long and Short conformations (~orres].)ondingly. The red curve is the linear combi-
nation of the blue and green curves with corresponding Boltzman weights. The red
curve in the bottom panel is the same as the red curve in the top panel. The grey
curve was extracted from MD simulations in which transitions between conforma-
tions occur. The black curve was obtained from the grey curve using our correction
method.

136

 

I 500 1 1 1 A 1 1

l T T I
l

1000 " 298 K

ITTTI
l

l 1 l l l

5 ()0

T T W‘T I
.--""'”‘
\
L)

if

0 l MA A L

‘ I
0 l 2 3 4
I 500 1 1 1 1

 

-—I
—

I—Ui
-O\
\l

 

l

' 298 K
1000 —

I
1

Pair Distribution Function g(r)

500

 

 

lllJllllll

T I hi I T I T I

0

 

 

 

—I

l l l
O l 2 3 4 5 6 7

Distance r (A)

—
_—

Figure 4.20: The blue and green curves in the top panel were calculated from eigenfre-
quencies and eigenvectors of molecular vibrations for the molecule in Long and Short
conformations correspondingly. They show the total PDFS of the molecule. The red
curve is a linear combination of the blue and green curves with corresponding Boltz-
man weights (mL(T) = 0.79 and IIIS(T) : 0.21). The red curve in the bottom panel
is the same as the red curve in the top panel. The grey curve, that also shows the
total PDF of the molecule, was obtained from MD simulations. The black curve was
obtained from the grey curve using our correction /convolution method.

137

 

1500 1 1 1 1 1 J

l()()() _ 800 K

Soo »— f I —

 

 

O
—"’)—4
A
K)!
O\
\l

 

1500 1 1

00
O
O
7?
I

Pair Distribution Function g(r)

l()()() —- _

500 —

 

 

 

 

 

 

Distance r (A)

Figure 4.21: The blue and green curves in the top panel were calculated from eigenfre-
quencies and eigenvectors of molecular vibrations for the molecule in Long and Short
conformations correspondingly. They show the total PDFS of the molecule. The red
curve is a linear combination of the blue and green curves with corresponding Boltz-
man weights (mL(T) : 0.62 and IIIS(T) : 0.38). The red curve in the bottom panel
is the same as the red curve in the top panel. The grey curve, that also shows the
total PDF of the molecule, was obtained from MD simulations. The black curve was
obtained from the grey curve using our correction /convolution method.

138

 

 

4.6.4 Correction to the Heat Capacity for CGHH molecule

Figure 4.22 shows the dependencies of the total energy (top panel) and the heat
capacity (bottom panel) of the (151114 molecule on temperature.

MD simulations were performed with the time step 1 ferntosccond. The Tinker
MD package provides as an output the values of energy averaged over 100 MD steps
(0.1 picosccond). The total simulation time used to calculate the average energy
values. was 1000 picoseconds for temperatures below 200 K. For ternprn‘atures above
200 K the total simulation time was 100,000 picoscconds. Thus, the average energy
values were found by averaging over 10.000 ('1‘ < 200 K) and 1,000,000 (T > 200 K)
different energy values. that are themselves the average values over 100 MD steps.
The MD simulations were performed at temperatures between 10 K and 1010 K.
Temperature increment. was 10 K for temperatures below 200 K. Simulations also
were performed at, temperatures (190. 195. 200. 205. 210) K, (290, 295. 300. 305, 310) '
K and so on. An additional simulations ware made at temperatures (240. 245. 250.
2:35, 200) K. In the beginning of every simulation run the molecule was in the Long
conformation.

In the top panel of Fig. 4.22 the results of MD simulations are shown as orange
circles (the orange curve connects orange circles and it is a guide for the eye). The
green and blue curves show classical and quantum results obtained in harmonic ap-
proximation from eigenfrequencies and eigenvectors of molecular vibrations. The red
circles (connected by red curve) show corrected MD results.

To obtain the heat capacity from MD simulations (the orange curve in the bottom
panel), the derivatives of the energy curve (4.39) were calculated numerically with
a temperature step 20 K for the temperatures below 200 K. For the temperatures
above 200 K the slopes of the lines that provide best fit (least square fit) to the 5
values of energy (for example at temperatures 190, 195, 200, 205, 210 K) were used

as the values of the heat capacities. Thus. the orange circles in the bottorrr panel

139

 

 

Energy (K/molecule)
A
O
O
O
O
1
I

 

 

 

N b) A U!
O O O O

p.—
O
lllllllllLlllllilLJllll
\
i
A

\
FIIIIIITVIIIII‘erIIIIIIII

Heat Capacitance ( in kb units )

 

 

 

O

I I I I I I

l
800 1000

I I I I I T I I I I

I l
200 400 600
Temperature T (K)

0

Figure 4.22: The dependencies of the energy and the heat capacity on temperature
for the CGH 14 molecule. Blue and green curves represent. the QM and CM results
for the harmonic approximation to the real potential. The orange curves and circles
represent the result. of the MD simulations. The red curves and circles show the
corrected MD results. See more detailed description of the figure in the text.

140

show the heat. capacity obtained from the orange circles in the top panel. In the sarrre
way. red circles in the bottom panel were obtained from the red circles (red energy
curve) in the top panel. The green and blue curves show classical and quantum heat.
capacitances obtained from eigenfreqrrerrcies and eigenvectors of molecular vibrations.

In MD simulations at low temperatures (T < 150 K), the change between the
Long and Short conforrrrations does not occur. At temperatures above 150 K and
below 200 K the molecule spends almost all time in the Long conformation. At
temperatures above 200 K the amount of time that molecule spends in the higher
energy (Short) conformation increases as temperature increases. This leads to the
values of the heat capacity that are on average higher than the values of the heat
capacity obtained in harmonic approximation from eigenfrequericies and eigenvectors
of rrrolecular vibrations. At. temperatures above 700 K when the molecule spends ap—
proximately half of its tirrre in the Long and half of its time in the Short conformations
the heat cari)acity becomes again basically purely vibrational one, i.e. the values of
the classical heat capacity. obtained from MD simulations, coincide with the number
of vibrational modes (there are 33-1 vibrational modes).

In another approach the heat capacity can be calculated from the energy fluctu-

ations using the formula:

 

(<E2>—<E>2)

l
C = —
'1' kfl"

(4.73)

However, this method can not be used with the "TINKER" program because the
temperature is introduced into MD simulations using Berendsen’s algorithm [80] that.
calculates correctly the average value of the energy, but does not calculate correctly
the fluctuations in the energy. If the temperature would be introduces using algorithm
of W'.G. Hoover [81] it would be possible to calculate heat capacity using fluctuation

formula .

141

4.7 Conclusion

In this work in an atterrrpt to develop an accurate technique for the calculation
of PDF for flexible rrrolecules we studied the role of quantum effects that are usually
ignored in traditional Monte Carlo or Molecular Dynamics simulations. We found
that it is very important to take into account the effect of zero-point motion for
the pairs of atoms that bring sharp peaks into the PDF calculated classically. We
(..leveloped a method that allows us to incorporate the effect of zero-point motion into
the PDF, calculated classically from Monte Carlo or Molecular Dynarrrics molecular
trajectories. and thus correct the classical PDF. We found that at large distances,
especially when the molecule becorrres flexible, its motion is almost classical.

In calculations of the total PDF for a conglomerate of molecules that are nec-
essary to make the comparison with the experimental measurements, interactions
between different molecules can be treated classically, since inter-molecular inter-
actions are rnrrclr weaker than intra—rr‘iolecular interactions. Thus, the situation with
inter—molecular interaction is analogous to the intra—molecular at large distances. Our
correction is important if an accurate PDF in the whole range of distances is desired.
It can help in reconstruction of the mutual orientations of different molecules, and
thus to determine possible constrains that can be caused by the inter'molt—écular inter—

actions.

142

 

 

 

Chapter 5

Summary

In the end, as it. is usually done, I want to summarize the main obtained results

and list possible further (.levelopments in the studied projects.

5.1 Charged lattice gas with a neutralizing back-

ground

Our study of the model for the ordering of two types of ions that can occur in
layered double hydroxides has lead to the following results. We found that the model
predicts that only at some particular concentrations of ions would the system have ho-
rrrogeneous distributions of the ions in the planes. At other concentrations the system
would separate into two phases with different homogeneous concentrations of ions in
them. It is interesting that. in spite of the original assumption that was made in the
model about. the homogeneity of the system we found a way to predict and study the
phase separations. Thus, from the point of view of the model itself it would be inter-
esting to consider some other physical situations to which this model can be applied.
On the other hand, it rrray be interesting to compare the results obtained in the frame

of this model with results of a rrrodel that directly allows phase separation. It is also

143

a possible direction for the (.levelopment of the model. It follows from the structure of
the layered double hydroxides that these phase separations can affect the stability of
the compounds or lead to the stacking of layers with different concentrations of ions
in them. We predict that. these phase separations can experimentally be observed
at relatively low temperatures, while at room temperatures these compounds should
have homogeneous concentration of ions in the planes. i.e.; to be stable. Thus while
these corrrpounds. as predicted. should be stable at roorrr temperature they may be-
come unstable as temperature decreases. So far there were no detailed experimental
investigations of the stability and the ordering of ions in these compounds versus ions

concentration and temperature. Thus it would be interesting to have these data.

5.2 Behavior of the PDF at large distances

The main conclusions from our study of the behavior of the pair distribution
function (PDF) at large distances are the following. We found that. the natural
assumption about the origin of decay in the PDF at large distances which, it seems
to us, was shared by many people who work in the area is incorrect. In particularly we
found and proved that an increase in the number of atoms at large distances does not
lead to the decay of in PDF. In other words there is no self-averaging of the PDF. The
PDF at large distances decay in amorphous materials because every atom at large
distances has different atomic environment. In agreement with the last statement, we
found that the PDF for the crystalline materials in which every atom has the same
environrrrent does not. decay at large distances. We understand that non-decaying
behavior of the PDF in amorphous materials originates from relation between the
amplitude of the random fluctuation in the radial density and the peculiarity in the
definition of PDF. Although qualitative behavior of the PDF in crystalline materials is

similar to the behavior of the PDF in amorphous materials, some quantitative features

144

are quite different. The behavior of the PDF in crystalline materials is closely related
to the Gauss circle. problem that CF. Gauss forrmrlated more that 150 years ago.
Both problems require further irrx-cstigations. Our result obtained on the behavior of
the PDF in amorphous case and qualitative similarities in behavior of the PDFS in
crystalline and the arrrorphous materials provide, from our point of view. a somewhat.
different perspective on the Gauss circle problem.

Another possible development. is related to the (.lerivation of a detailed analytical
relation between the excluded volume and the amplitude of the oscillations in the
PDF at large distances in the arrrorphous case.

Since our results can be used to test the amount of dislocations in crystalline
materials and to test the instrumental resolutions in X-ray or neutron diffraction
experiments there are two experiments that are highly desirable.

Thus, it would be valuable to compare the PDFS obtained from the two crystals
made of the same rrraterial. One measurement should be done on a crystal of very
good quality, while another on the crystal with significant amount of imperfections in
the crystal structure. These measurerrrents should be done on the device that provides
good and stable instrumental resolution.

Other two measurements should be done on a single crystal with crystal lattice
of a very high quality, but on two different machines that have different instrumen-
tal resolutions. Although information about the instrumental resolution can be ex-
tracted directly from the X—ray or neutron scattering intensity, I think that these
measurements would be valuable in the modelling of the mathematical forms of the
instrumental resolution functions that. are used to convolute perfect PDFS obtained
from the model structures in order to obtain better agreement between modelled and

experirrrental data.

145

5.3 Quantum corrections for the PDF calculated
classically

In this work, we developed a tecl'mique that allows one to incorporate the quan—
tum effect of zero point motion into the PDF. which can be calculated for molecules
in liquid or gaseous phases, using classical MD simulations. This technique should
lead to an improved agreerrrent between the modelled and experimental PDFS. The
technique can also be useful for calculations of the specific heat.

The next step would be to calculate the PDF for a real liquid, which consists
of relatively complicated molecules. For this purpose, it would be necessary at first
to run MD for a conglomerate of rrrolecules. Then it would be necessary to extract.
from the molecular trajectories the probabilities to find pairs of atoms at a particular
distance. These distributions for some pairs of atoms should be corrected, using the
developed technique. Summation of the PDFS for different pairs should lead to the
total PDF for the conglomerate of molecules that model PDF for the real liquid.

I think that at this stage. it would be valuable to have experimental PDF for the
liquid that consist of the C '6 H 1.; molecules. From these data. that we have it. would be
possible to estimate the role of the quantum corrections by comparing the calculated
PDFS for a. single molecule at different temperatures with experimental PDFS. High
qrrality experimental PDF would also be valuable for the estimation of the level of the
structural information that can be extracted from an experimental PDF on liquids

which consist from flexible molecules.

146

Appendix A

Behavior of pdf at large distances

A.1 Random case in d-dimensions

A. 1 . 1 Ensemble averaging

Now we consider the case of completely random media when points are randomly
distributed in space. In particular, it means that two points can be very close to each
other; the situation that can not occur with atoms in crystals or glasses. we also
assume that 0,-1- = o for every pair of sites i and 3'. In this case it is easy to make
derivations for the general case of d-dimensions. The PDF function is defined in

d-dimensional space as (3.9).

Cd,(r) : rd£1

 

[pd-4r) - pol- (Al)

The origin of the exponent of r will be clear from the following. The area of the surface

of d-dimensional sphere of radius r can be written as Sd( r) : erd‘l, where 0,; is

the total solid angle. Thus 91 = 1, $22 = 27r, {23 = 17c Let us assume that the space
,.d-l

is divided into very sr‘nall boxes. The volume of 2th box is (IV, = d‘lr, = dildo,- (In.

Tire number of particles in the 27th box can fluctuate, but the average (averaging over

147

 

different distributions) number of sites in the ith box is d; = pot/I”, = p,,(l.5'd(r,)dr, =
p(,(1§2d,r:i‘1(lr,. The usual relation (3.14) is applied to the average fluctuation for the
number of sites in the box: [(172, — Err—J3 = [(dn,)'2 — (Tn—J2] :2 (TI—2:. In the following.
while calculating pd,(r). Cd,(r). and Gil-(r) with respect to a particular site we will
drop indexes a’ and i to simplify the notations. Thus for the given distribution of sites

the radial density with respect to particular site can be written as:

r l‘ — 7‘1)2
[)(r) = {exp(—(—20§—]}dn,. (A2)

1 Z 1
Sd(r) i V 27TH;2
Ensemble averaging of (A2) leads to

 

 

I)

w 1 e—nr mamamrwn
[)(T) = ———‘——‘{e.\'}')[— , . l} ..d—l
0 /27T02 20', QatI

Since p(f}) 2 p0 there is no angular dependence and integraticm over the angles cancels
Rd in denominator. For small a the contribution to the integral comes only from the

region r,- E’ 7‘. Thus we conclude that [)(r) 2’ pa. Further:
Gao=rctmi—ar=rtwao—pa (as

In order to calculate p'-’(/‘) it is necessary to calculate drudnj:

 

 

———, 1 (r — r,)2 1 (r — r1)2 (Initial-
r 2 = —— ex — . ex _) — . . A.4
p( ) ;W{ PI 20; llm{ If 20f l} 550,) ( )

In order to calculate (172,-(1rrj we note that:
(17mm!- E [(d-n,)2 — (Efldij + m fin—J,

where the term in square brackets is equal to (1n,- = pOdI/i. If i 7t j then the sunnna-

. . . . . ‘) .
tions over 2 and ] 111 (A4) can be performed separately leading to p5. Thus we can

148

 

 

 

write:

 

, p0 (r — n?) (112,1,7' ‘1— (1r, . K
per: arr) / 2Mae-x1» pr— 0. 1} 9d,, , +p3. (A...)

 

 

 

If the difference between I dand I‘d in the integrand is ignored. integration leads to

 

. 0 1 [)0 ,2
/)(r) _ 2fiSd(I‘) U + [)0

Finally. from (A3) (Sd(~r):f2,1r‘d 1) we conclude that

o —.) U I
SHAH) : “(r)— _

— . A.
p0 Qfigd ( 6)

 

|
S

In the last. expression we introduced index 2'. again to emphasize that 03,0) was

calculated with respect. to a particular site.

A.1.2 Integral approach.

It is useful also to consider another approach for the calculation of G3,(r). In
the following, for simplicity of notations, we drop index (1. Instead of averaging over
different distributions of sites one can calculate, for a particular distribution, the

quantity:

. 1 ”'2 . I
< (If 7‘ >= —_—/ Cf I' (17" A."

In the following we explicitly show that both approaches lead to the same answer

in case of the random distribution of sites. It follows from (3.7,3.9.A.1) that:

122
I = / I'd'l[p,2(r) — 2p,(r)p0 + pild'r. (A8)
R1

149

Further we introduce:

R2
12 :/ -I‘d"12p,-(r)p0dr (A.10)
R1
and
122
132/ rd lp;(1].(A.11)
. R1
It. is obvious that:
R2 a— 1 2 P Qd d d
Ir: 1=—— " R —H, 0 — "<Ng2 >, A.12
.5 [h I po(r 52—; dli lip }_ Sid ( 1

where < N R2 > is an average number of particles in the spherical annulus between
R1 and H2.

For [2 using expression (3.8) for /),-(r), in the assumption that 0,]- : a, we write
(If indices of sunnnation are not shown it is assumed that the summations goes over

suchj that R1 < 73] < H2):

R2
12 = / I‘d—l2p(r)p0dr (A.13)

ngo/Wl—l Z .——.{exm— ————-———( "’frrdr
9d 12, Qde—r . 27ro2 202

2-0 ‘Tij .
: If; Z/II, V2WU2{exp[ (r 20 2 )2 ”(if

#i
a. 21

4'32
<I\RI >.

 

 

 

II?
[\D
‘o
o

 

I?

 

 

 

 

 

It. was assumed during the last. derivation that (95 sign) 0 << R2 — RI so that basically
only sites inside of the annulus contribute to the value of the integral. This assumption
is also used everywhere in the further (.lerivations. We used the sign 2 since for the
particular distribution of sites the number of sites inside the spherical annulus can
deviate from its average value < N]??? > with standard deviation (/< N312 >.

For 11 from (A.9.3.8) we have:

Hi (I‘ — r0)2 1 (7‘ — I‘,I.)2 (1r

[— _ ix — .
1 22/ ,, ”fldmewl 11,.-.

jzfiz' kyéi

 

 

This expression can be split into the two sums 11 = [11 + [12. In the first sum [11:
j = k while in the second [12: j yé k. Thus (If indices of sunnnation are not shown it.

is assumed that the summations goes over such 3' that. R1 < I‘U' < R2):

 

1 1 [R2 {e 1— (""”"")211 d" (A14)
= -——.—- E x , .
ll QfiQZiU R1 \/—0— I) 02 7"1—1
= 2—__\/—,1§220 Z..— d—r_ 1

~——-— ————Q 'r.._ pdr~
2 ,(1—1 d 'li 0 I]
Qfiildo R1 71‘)

_1_ p_._.
= R —R .

It follows from (A.14) that expression for 112 can be written as:

1.. =zzmepr—Mr (A15)

3' kséj

1 ”2 1 1 1 -r,,-+r,-,. .,
-— —— . —+ »— _—_ (I,
523/3 I‘d-1mm“ '4’ 2 ) l ’“

T119 integral in the last expression has non—zero value if r e: (In-j + TM.) / 2 i 0. Since

0‘ ' . . . . . .
18 small (A.15) can be rewrrtterr as (If indices of sunnnation are not shown rt rs

151

assumed that the summations goes over such j that R1 < r1,- < R2):

III?

 

1 (7‘1) _Vr1k12 2 (1—1

(A.16)

Fm‘ther we substitute summation over k by integration assuming mean distribution

of sites (If indices of sunnnation are not shown it. is assumed that the summations

goes over suchj that R1 < I‘,, < Rg):

R2 . 2

2 : ’1’ — 7111‘)

I ”_>_: m ___J______
12 ./Rl V4I:—U—2{(Xp[_ 402 f}

1 2 d— 1 ‘
Qimfgo Qd’ik lPod’Ik
p” R.
g— 1 — — d< I\,
gr, 2
Since 11: [11+ 112:
1 [)0 1 1 R,
1 If? _H + '— 0 < ‘AR' 2 > '
122———\/—52d (CT 1) izdp‘ R1

Finally, for I = 11 — 12 + 13 we write:

1 p0 [’0 I? _2__/)0
I 9: R; —R N 2
77de 0(21)+Q_d < R1 ad
__1_
2 —°.[ —11

YR +p0
< I\Rf> +Q—d

(A17)

(A.18)

< N312 > (A.19)

. . . .R. . . 7 .
Note that In the expression above. the value of the < IN“: > 18 determined w1th

precision (/< N512 > (that is the origin of the 2 sign) and that all terms in which

< N312 > is present directly cancel each other. It is also indirectly present in the

remaining term (see (A.15)). However it does not affect the domination of the main

value of the remaining term.

152

‘3' "a...

 

Thus. for gflR) we have:

gm =<G§.< >>32 1

p0 Qfifld ,

 

(A20)

which is the same expression (Ab) that was obtained from ensemble averaging.

A.2 Amorphous case in d-dimensions: integral ap-

proach

A.2.1 Difference with the random case from the integral ap-

proach: detailed derivation

The derivation for the random case with forbidden volume goes in the same way
as for the completely random case. The only difference comes from the fact. that when
there is forbidden volume around every point expression (A.18) should be modified
into (If indices of summation are not shown it is assumed that the summations goes

over such j that R1 < n, < Hg):

 

132 _ T‘k)2

~ /. gem ——e—]} M
_1: 21 l7rcr'2 402

—— Gd 7:1 1 — 'Ad—l Drink.
91,7, 73-,- +2,,,,ld_1( 7 )p
This change comes from the fact that if. in an amorphous case, there is a point 3' at a
distance 7",,- then point. k: should be placed in such a way that In, —r,k| ~ 0. Otherwise
the exponent in expression (A22) will vanish. But in an amorphous case not the whole
volume Qd'rfj'ldn, is available for placing point 1:. There is forbidden volume around
the site j: V; E M2d’rij, where /\ can be considered as the distance between the

nearest neighbors and 7 is the coefficient that depends on the structure. Thus, see

153

 

 

 

Fig.3.8. effectively eridrm turns into ($2,111, 1 — c,r)\‘1‘1)dr,vk. Further (compare with
(A.18)) (If indices of summation are not shown it. is assumed that the summations

goes mer such j that R1 < 7‘, < R) )z

m afl<Wp mm

. ‘2 1—1 A d—l
— ’z'k) 2( ° ,v‘/\ Podlzk

— Z/Rl 117F02{CQ(N VP )—([ "102 J} ("'0' + "z'lcld—l l

-, 1 1
~_<wu>_fl;&zi_ 4
$2,, R1 9,21 ,rd—l'

ij

Using mean density approximation, we again switch from the sum to the integral:

 

. (1—1 R2
p_o_ ’l’/\ pa 1 ._ . .
[12 ff." Q<Av1ljlz > ——S)—2/ TYQJ": 1,00drl-J- (A.23)
d “d R1 '1‘}
d—l .2
pa [2 “M 1)..
:—-m\2>—————H—R.

Thus, in an arrrorphous case for I = [1 -— 12 + 13 we have:

1 p0 /'/\d—1/)2
— — ———O R1 — H. A24
2[2—l\/_S2d 0' Qd l( 2 1) ( )

So that. for < 0341') > we get:
< (1%) >2 27:77,?[1 — 2f7Ad-1poa]. (A25)
Since p0 = UV:

0
9321“?) =< (1dl(, )> — g

l 0‘
_———1—2 A-. A2"
Po Qfighl fi/Al ( b)

154

 

A.3 1d crystal

From the definition of radial density in 1d (3.8) it follows that:

 

l (r — rm)2 ‘ _
[)(r) : 2,703 "22—; exp[ 202 l, (A 2()

while from the definition of PDF in 1d (3.7):
02m = W) — awe) + p:- (A28)

Due to the iwriodicity of the lattice, in order to find < Car) >, it is enough
to perform integration (3.10) over any interval of length (1.. Because of the structure
of expression (A27) for [)(r) and its physical meaning, the integral of p(-r) over any
interval of length a should be equal to poa = l, i.e. < p(-r) >= p0. Thus, for any R,

we can write:

(7“ — 7m)

R+a 90 1 2
ex —————,—.—— (17" = 1, A29
/R 71ch W{ pl 202 l} ( )

while < 02(7') >:< p2(7‘) > —p3. From (A27) follows that:

 

(r — nu)2 + ('I‘ — ma.)2
202

 

 

U2{exp[— ]} (A.30)

=2:

ll: -CXJ 7712’“;
Since

‘ 71+ TN (1 . n, — 7n 202
(r - ”(1)2 + (1‘ — ma)‘2 : 2(7‘ _ LTLV + (—2)_

expression for p2(7‘) can be rewritten. If i = n — m and j = n + m. then:

°° (A)?

=2: first Pl nzrepi 72—» M

 

155

 

 

 

 

In the terms that arise from the last. expression 27 and j should be both even or odd.

Thus for < p20“) > we get (compare with (A.‘29)):

 

 

 

1 «+0 ') :30 1 [lug
__ " ‘ 1‘ : P" l__r_ . .‘\.32
(I /R I) (I )(I ’22—:30 2fi0a {( \I)l 402 l} ( )
From that it is follows:
1 (7 1 p 00 71'qu
< C2 ' >= ———O l— ‘2 7c— 2——0 .' —+. A33
(7) ‘2 TI'O’[ a]+ ‘2 7r0 ;e\l)[ 402] ( )
Finally:
'2 0° 2 2
(law) . 0 n a ‘ g
. = l— 2 7r— + 2 ex — . A34
9w) [ Wu] 2 pl 402] ‘ l

156

 

 

 

 

 

 

Appendix B

Cut-off for the Morse potential

In order to find (7121mm and < U(T) > for the Morse potential, it. is necessary to
perform integrations (4.14, 4.17, 4.38) with the Morse potential used as U(;r). For the
Morse potential as r —+ 00 the value of potential goes to const, i.e. exp [—U(.r)/ (L‘bT)]
also goes to coast. Thus the integrals under consideration diverge as 7' —+ 00. This
divergence means that in infinite time, a particle would escape from the potential at
any temperature. Due to this (.livergence it is unclear what. should be the upper limit.
of the integration. However, if temperature is low enough, there is a range of upper
cut-offs (sufficiently big, but not too big) that would lead to the same results with
respect to 03,“,(7‘) and < LI’,\1(,~M(T) >.

Figure B.1 shows dependencies of Uilc'adT) and < Mum/(T) > on temperature
for different. values of cut—off. Thus, we see that if T S 10000 K any reasonable value
of cut-off leads to the basically the same results. Thus for the curves on the Fig.4.4,

4.5, 4.6 cut. off 2.5 (A) was used.

 

 

 

 

 

 

 
  
   
    

 

 

 

   

   

 

 

 

 

l 1% l l l l l l l l l i l_l l l l 1
:lSOOO lllllllllllLlllllil C
0.04 1 j l r
: _ — 1.5 (A) /- :
10000~ 2'0 (A) '/ — E
j "i __ 2.5 (A) ” ego :
0.03 — 5 __ * /\\ __
: - 3.0 (A) _ :> :
_ .. _ V _
NA 3 5000— — /:
°< s ‘ * i
V -— _. p— .—
N —« I»—
D 0 02 j ‘ l— T
E] 0 IIIIIITTTIIIIIIIIIT :
I 0 10000 20000 :
3 Temperature T (K) :
0.01—j :—
0 I I I I l I I I I T I I I I l I I I I
0 5000 10000 1 5000 20000

Temperature T (K)

Figure 8.1: Dependencies of 0310““) on temperature (see Fig.4.4) for different
values of the upper cut off in the integrals (4.14, 4.17). Thus, if T S 10000 (K) any
cut-off in the interval (1.5:30) Awould lead to basically the same result. The inset.
shows how < U MC M(T ) > (in temperature units) depends on temperature. Thus for
the potential energy any upper cut off (see (4.38)) in the interval (1.5230) Ais also
suitable.

158

 

 

 

 

 

Bibliography

[1] Y. Xiao, M. F. Thorpe and J. B. Parkinson, Phys. Rev. B 59, 277 (1999).

[2] V.A. Levashov, M. F. Thorpe and B. W. Southern, Phys. Rev. B 67, 224109
(2003).

[3] H.F.VV. Taylor, l\~’Iineral. Mag. 39, 377 (1973).
[4] D.L. Bish and GVV. Brindley, Am. Miner. 62, 458 (1977).
[5] G.VV. Brindley and S. Kikkawa, Am. Miner. 64, 836 (1979)

[6] SA. Solin, D.R. Hines, GT. Seidler and M.M.J. 'Iieacy, J. Phys. Chem. Solids
57,1043(1996)

[7] DR. Hines, GT. Seidler, l\I..\*I.J. Treacy and SA. Solin, Solid State Commun.
101,835(1997)

[8] SA. Solin, D.R. Hines, SK. Yuk, T.J. Pinnavaia and MF. Thorpe, J. Non-
Crystal. Solids 182, 212 (1995).

[9] A. Audemer, A. Delahaye, R. Farhi, N. SacEpee and JM. Tarascon, J. Elec-
trochem. Soc. 144, 2614 (1997).

[10] S.U. Falk and A.J. Salkind, Alkaline Storage Batteries, Wiley, New York, 1969.
[11] R. Barnard and CF. Randell, J. Appl. Electrochem. 13, 97 (1983).
[12] AD. Vidal and M. F iglarz, J. Appl. Electrochem. 17, 589 (1987).

[13] P. Elumalai, H.N. Vasan and N. hr’funicharidraiah, J. Power Sources 93, 201
(2001)

[14] B. Liu, X.Y. Wang, H.T. Yuan, Y.S. Zhang, D.Y. Song and 2.x. Zhou, J. Appl.
Electrochem. 29, 855 (1999)

[15] C. Faure, C. Delmas and P. W illmann, J. Power Sources 36, 497 (1991).
[16] GA. Caravaggio, C. Detellier and Z. Wronski, J. Mater. Chem. 11, 912 (2001).

[17] F.Ducastelle, Order and Phase Stability in Alloys, Elsevier Science Publishers
B.V. Amsterdam, The Netherlands, 1991.

159

 

 

 

[18] A.J. Berlinsky, W.G. Unruh, W.R. McKinnon and RR. Haering, Solid. State.
Commun. 31, 135 (1979).

[19] W. Li, J.N. Reimers and JR. Dahn, Phys. Rev. B 46, 3236 (1992).
[20] AH. Thompson, J. Electrochem. Soc. 126, 608 (1979).

[21] M. Winter, J.O. Besenhard, ME. Spahr and P. Novak, Adv. Mater. 10, 725
(1998)

[22 T. Zheng and JR. Dahn, Phys. Rev. B 56, 3800 (1997).

[23] SA. Safran and DR. Hamann, Phys. Rev. B 22, 606 (1980).

[24] SA. Safran and DR. Hamann, Phys. Rev. Lett. 42, 1410 (1979).
[25] HA. Kramers and CH. VVannier, Phys. Rev. 60, 252 (1941).

[26] M. E. J. Newman and GT. Barkerna, Monte Carlo Methods in Statistical
Physics, Oxford University Press Inc., New York, 1999.

[27] AC. l\-ilaggs and V. Rossetto, Phys. Rev. Lett. 88, 196402 (2002).
[28] ER Wigner, Phys. Rev. 46, 1002 (1934).

[29] G. Meissner, H. Namaizawa and M. Voss Phys. Rev. 13, 4 (1976).
[30] L. Bonsall and AA. Maradudin, Phys. Rev. B 15, 1959 (1977).

[31] This exact result only agrees with the approximate expression in equation (11)
of reference[1] in the limit a: ——> 0

[32] F. Zernieke, J.A. Prins Zeit Physik. 41, 184- (1927).
[33] P.J.W. Debye, H. Menke Physical. Zeit. 31, 797 (1930).
[34] BE. Warren X—ray diffraction (Dover Publications 1990).

[35] J .S. Chung, M.F. Thorpe Local atomic structure of semiconductor alloys using
pair distribution function Phys. Rev. B 59, 1545-1553 (1997).

[36] M.F. Thorpe, J .S.Chung, S.J.L. Billinge and F. Mohiuddin-Jacobs Advances in
Pair Distribution Profile Fitting in Alloys in Local Structure from Diffraction,
157-174, Ed. By S.J.L. Billinge and M.F. Thorpe (Plenum Press, New York,
1998)

[37] M.F. Thorpe, V.A. Levashov, M. Lei and S.J.L. Billinge Notes on the analysis
of data for pair distribution functions, 105-128, Ed. By S.J.L. Billinge and M .F .
Thorpe (Kluwer Academic / Plenum Publishers, New York, 2002)

[38] DA. Dimitrov, H. Reder and AR. Bishop Phys. Rev. B 64, 14303 (2001)

160

 

 

 

 

[39] RH. Toby and T. Egami Acta Cryst. 48. 336 (1992)

[40] For a summary on Lattice Point. lecn‘y, see
http: //mathworld.wolfram.com/ PointLatticehtml

[41] E. Kréitzel Lattice Points, Kluwer Academic Publishers (1988)

[42] P. Erdés, PM. Gruber, and J. Hammer Lattice Points, Longman Scientific and
Technical with John Wiley and Sons, Inc., New York (1989)

[43] C .F. Gauss De nexu inter multitudinem classium, in quas formae binariae secundi
gradus distribuuntur, earumque determinantem. VVerke, vol.2, pp.269-291 (277)

[44] For a summary of what, is known about the Gauss Circle Problem, see
ht t p: / / mathworld.wolfram.c0m / GausssCircleP roblern. html

[15] L.D. Landau, EM. Lifshitz Statistical Physics P1 (1999)

[46] S. Torquato, T.M. Truskett and PG. Debenedetti Is Random Close Packing of
Spheres Well Defined? Phys. Rev. Lett., 84, 2064 (2000)

[‘17] “CM. Visscher and M. Bolsterli Nature (London) 239, 504 (1972) See also
remarks published in Nature (London) 239, 488 (1972)

[48] J. Toboclmik and PM. Chapin J. Chem. Phys. 88, 5824 (1988).
[49] F. W'ooten, K. Winer, D. W'eaire Phys. Rev. Lett. 54, 1392 (1985).
[50] H. He Ph.D. thesis, Michigan State University (1985).

[51] GT. Barkema and N. Mousseau High-quality continuous random networks.
Phys. Rev. B 62, 4985-4990 (2000).

[52] R.L.C. Vink, G.T. Barkerna, M.A. Stijnman and RH. Bisseling Towards device-
size computer models of amorphous silicon. Phys. Rev. B 64, 245214 (2001).

[53] MN. Huxley Exponential Sums and Lattice Points. Proc.London lVIath. Soc.
60, 471-502, 1990

[54] G. H. Hardy Quart. J. Math. 46, 283, 1915

[55] E. Landau fiber die Gitterpunkte in einem Kreise. (Zweite Mit-
teilung.) Nachrichten der K. Gesellschaft der Wissenschaften zu Géttingen.
Mathematisch-physikalische Klasse. 161-171 (1915).

[56] PM. Bleher Duke Mathematical Journal 67, No.3 p.461

[57] PM. Bleher, F.J. Dyson and J .L.Lebowitz Non-Gaussian energy level statistics
for some integrable systems Phys. Rev. Lett. 71, No.19 p.3047-3050 (1993)

161

 

 

 

 

 

[58] J .S. Chung. M.F. Thorpe Local atomic structure of semiconductor alloys using
pair distribution function Phys. Rev. B 59, 1545-1553 (1997).

[59] R. Car and M. Parinello Unified Approach for lVIolecular Dynamics and Density-
Functional Theory, Phys. Rev. Lett. 55, p.2471-2474 (1985)

[60] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R.
Cheeseman, V. G. Zakrzewski, J. A. Montgomery, Jr., R. E. Stratmann, J. C.
Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain,
O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cannni, B. Mennucci, C. Pomelli,
C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P. Y. Ayala, Q. Cui, K.
Morokuma, P. Salvador, J. J. Dannenberg, D. K. Malick, A. D. Rabuck, K.
Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, A. G. Baboul, B. B.
Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gornperts, R. L.
Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M.
Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. W'ong, J. L. Andres,
C. Gonzalez, M. Head-Gordon, E. S. Replogle, and J. A. POple, Gaussian, Inc.,
Pittsburgh PA, 2001.

 

[61] M.J.S. Dewar, E.G. Zoebisch and BF. Healy Gaussian 98, Revision A.11.1, J.
Amer. Chem. Soc. 107, 3902 (1985)

[62] L.D.Landau and EM. Lifshits Quantum Mechanics: Non-Relativistic Theory,
Volume 3, Butterworth-Heinemann; 3rd edition (February 1997)

[63] L.D.Landau and EM. Lifshits Statistical Physics, Volume 5, Butterworth-
Heinernann; 3rd edition, Part 1 (February 1999)

[64] PM. Morse Phys. Rev. 34 57—64 (1929)

[65] N. L. Allinger, Y. H. Yuh and J.-H. Lii, ”lVIolecular Mechanics. The MM3 Force

Field for Hydrocarbons. 1”, J. Am. Chem. Soc., 111, 8551-8566 (1989)
[66] J .-H. Lii and N. L. Allinger, ”Molecular Mechanics. The MM3 Force Field for

Hydrocarbons. 2. Vibrational Frequencies and Thermodynamics”, J. Am. Chem.
Soc., 111, 8566-8575 (1989)

[67] J .-H. Lii and N. L. Allinger, ”Molecular h’lGCllaIllCS. The MM3 Force Field for
Hydrocarbons. 3. The van der Waals’ Potentials and Crystal Data for Aliphatic
and Aromatic Hydrocarbons”, J. Am. Chem. Soc., 111, 8576-8582 (1989)

[68] N. L. Allinger, H. J. Geise, W. Pyckhout, L. A. Paquette and J. C. Gallucci,
”Structures of Norbornane and Dodecahedrane by Molecular Mechanics Calcu-
lations (MM3), X—ray Crystallography, and Electron Diffraction”, J. Am. Chem.
Soc., 111, 1106-1114 (1989) [torsion-stretch]

[69] N. L. Allinger, F. Li and L. Yan, ”Molecular Mechanics. The MM3 Force Field
for Alkenes”, J. Comput. Chem., 11, 848-867 (1990)

162

 

 

 

 

[70] N. L. Allinger, F. Li, L. Yan and J. C. Tai, ”.'\*Iolecular l\‘lechanics (I\-Il\-’l3)
Calculations on Conjugated Hydrocarbons”, J. Comput. Chem., 11, 868-895
(1990)

[71] J .-H. Lii and N. L. Allinger, ”Directional Hydrogen Bonding in the l\/‘Ii\-’13 Force
Field. I”, J. Phys. Org. Chem., 7, 591-609 (1994)

[72] J .-H. Lii and N. L. Allinger, ”Directional Hydrogen Bonding in the l\-‘Il\'I3 Force
Field. II”, J. Comput. Chem., 19, 1001—1016 (1998)

[73] See the Internet web-site of the ” TIN KER” program:
http://dasher.wustledu/tinker/

[74] P. Ren and J. W. Ponder, J. Phys. Chem. B, 107, 5933-5947 (2003)

[75] R. V. Pappu, R. K. Hart and J. W. Ponder, J. Phys. Chem. B, 102, 9725-9742
(1998) '

[76] M. E. Hodsdon, J. W. Ponder and D. P. Cistola, J. Mol. Biol, 264, 585-602
(1996)

[77] M. J. Dudek and J. W. Ponder, J. Comput. Chem., 16, 791-816 (1995)

[78] C. E. Kundrot, J. W. Ponder and F. M. Richards, J. Comput. Chem., 12,
402-409 (1991)

[79] J. W. Ponder and F. M. Richards, J. Comput. Chem., 8, 1016-1024 (1987)

[80] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola and J. R.
Haak J. Comput. Chem., 8, 1016-1024 (1987)

[81] W.G. Hoover Phys. Rev. A 31 1695-1697 (1985)

163

 

]

 

 

 

llllllllllllllll][[l[l][lll[l]l]l]ll

3 1293