Phase transitions of two-dimensional potassium chloride

Volume 35, numk

4

CHEhiICAL

PHASE TRANSITIONS RW. HOCKNEY Computer

Science

Received

FWYSICS

OF TWO-DIMENSIONAL

LETTEELS

POTASSILN

1 October

1975

CHLORIDE

and S.P_ COEL * Departmelft,

R&ding

University,

Reading

RG6

2dF,

UK

17 June 1975

The m&i&,

crystallization and boiling of a computed “nmplc” of two-dimensional potassium chloride has been studied of molecular clynamics. A film showing the entire proce~, starting from &e shble solid state with a perfect lattice up to the emcr_eence OF the vapour phase characterized by the dominance of single pairs of potassium and chlorine ions, has been made. Some still pictures from the film are included. Values are obmined for the melting point, boiling point, expansion coefficients, specific and latent heats. These values are shown to apply to a wide range of two-dimensional materials

by the method

and arc compared

with the three-dimcnsionzl

values.

1. htroduction

ionic systems with the slowly-decaying lombic interaction.

Althou& a great deal of work in many fields has been done throu$~ molecular-dynamics (MD) computer simuhtion (see for example, the resource letter recently prepared by the authors Cl] on this topic), the MD method has not been used extensively to study systenis like ionic solids or liquids, in crhich long-range forces, are important. The limited studies done of ionic systems

through

to assembiies

the use of this technique,

of a small

number

of particles,

are confined usually

.64 or 2 16 [2,3 J _This arises because the number of computer operations that is required to simulate a system of N particles, u&g conventional methods of molacular dynamics, is proportional to N2. The computer simulation of systems of a large number of particles is therefore very expensive. One may contrast this situation with that of a Lennard-Jones (LJ) sy-stem such as argon, in which the force of interaction between atoms decays rapidly with increased separation. In this case the interaction of atoms, separated by more than a certain cut-off radius, may be ignored without significantly affecting the accuracy of the results. Using this approximation, systems of seveial thousand argon atoms may be simulated [4]. However this approximation is not valid.in the study of

Recently

Resent

address:

Kurukrhe&,

Physics India.

Depar&nent.

Kurukshetrn

have

been

made

cou-

in the techniques

available for the simulation of ionic systems based on the use of a mesh for th.e calculation of the long range forces. The particle-particle/particle-mesh (PPPM) scheme devised by Hackney et al. [5-g], which is used in this investigation, enables ionic systems of lo4 or more particles to be simulated. Flowever, in order to permit a rapid survey of all the phase transitions in a reasonable computer time, we have used a relative small sample of 1024 ions for most of the results of this paper. The survey includes the solid-liquid, liquidsolid, and liquid--vapour transitions and is available as a 16 mm fdm from the authors. In order to assess the accuracy of the results the solid-liquid transition has also been studied under more careful conditions with a larger sample of 6400 ions. The study of the properties of the two-dimensional material, which can only be done by computer simulation, is of considerable interest in its own right. In particular we direct our attention to the differences be!ween the properties of the two-dimensional and three-dknensional systems.

2. The PPPM

l

advances

long-range

scheme

University, The

PPPM scheme has been fully described

elsewhere

Volume 35, number 4

CHEMICAL

PHYSICS

[.5-71. Its main features are its freedom from the li_ mitations of the conventional MD method and its capacity to study systems of up to 20000 particles. It has already been used to simulate the condensation of 10000 potassium chloride ions [S] The scheme in brief is as follows. The force of interaction between particles is divided two parts. The long-range electrostatic part, which is determined by the solution of Poisson’s equation on a fixed fmite-difference mesh covering the region of calculation; and a short-range repulsive part which is found, as in conventional MD, by summing contributions from particles within :I certain radius. In other respects the calculation proceeds step-wise in time in the same way as conventional molecular dynamics. The hybrid method thus combines the mesh techniques, long used in plasma physics for the Coulomb interaction [9], with the conventional MD method and gets the best of both worlds. It is however only practical after some rather subtle tampering with the Fourier transforms during the solution of Poisson’s equation into

KEITERS

1 October

usual rules for the multiplication of signs in the product (vii>are understood to apply. The first term represents the coulombic interaction and the second ten-n the strong short-range repulsion which dominates when

the electronic shells of ions overlap. This simple analytic form of the force !aw is adopted to allow the derivation of and the comparison with theoretical results, whilst retaining the essential qualitative features. This law will be referred to as a (l-10) force law, where the digits refer to the powers of r in the coulombic and repulsive terms respectively. The equations of motion of an ion are: ,prid’&dt2

i = +,-,

=

excepting the ion for which the equation is written. We have taken the fo!lowing values for the constants: ‘*

=

f;-

ri’=

In this work, we assume that the region of interest is covered by a square mesh consisting of(256 x 256) cells. The 512 potassium (K+) and 5 12 chlorine (Cl-) ions are placed, alternately in 32 rows, in 3 32 x 32 square lattice in the middle of the region under consideration (see fig. la). The square lattice is chosen as it represents the minimum potential energy configuration. Periodic boundary conditions are imposed on the twodimensional region of space in which the ions are kept. An ion going out of t!ze region in one direction is assumed to enter it from the opposite side. The periodic repeat unit is the square box drawn in figs. 1 and 2. The crystal occupies only the central = iO% of the computed area and the external forces from the image systems on the test sample are negligible. For this reason the experiment can be considered as being conducted at a constant external pressure of zero. The inter-ionic force is assumed to have the following form: Fij(T) = Q_&r)

(1 f (ii)(!-$,9},

(1)

where i and j = +.for R’ icns, = - lor Cl- ions and zhe

(2)

where -s! is the coordinate of the ion in question and the summations are taken over all ions of the indicated type,

=

-+ = 3.14 x 1W8 cm,

rQ

[71. 3. The model

1975

2.66

x

10e8 cm,

ri-

(3) = 3.62 x 10m8 cm; ’

where we use the notation r 0, without superscripts, to denote the equilibrium separation for a pair of ions of

opposite polarity; m + = 39-l/>,

amu/cm,

)~z- = 35.457/2r, arnulcm,

(4) rl+ = 4.8 x lo- “/2ro

esu/cm,

q- = -qf

esulcm,

where amu = 1.66

x lo-*” g. The above constants were chosen to represent KC1 as close as possible in a two-dimensional simulation T!te amount of mass and charge per unit length is largely arbitrary. Distributing the atomic quantities over a distance (2ru): as we have done here, ensures that tire vibration frequency for an isolated pair in the two-di-

mensional model is the same as that for the equivalent three-dimensional model with a force :aw: F(r) = (&i/G)

{l + (zj~(r~lr>g}_

(5)

The time step chosen was lo-” s and the space mesh distance H =r0/2 = 1 Si A. The initial crystal size is 120 x 120 AZ and the region of calculation 400 x 400 A’. Eqs. (1) and (2) may be put in dimensionless form by expressing distances in terms ofro, and times in units of

501 ..

CHEMICAL PHYSICS LETTERS

.-Volume 35; number 4

of four parts - the four corners - in an expanded but nearly ordered form, interacting rather weakly with one another. This behaviour of the lattice seems to be similar to that of a two-dimensional square assembly of connected springs which can withstand only limited tension. The springs along the midclle lines nearly give way when

., One then obtains:

p

F

tt = t].r.

r&,

Since p+ = /3’- = 1 ther e are only three constants describing the problem namely or-, Q’” and fl--. Hence, the resuhs quoted here apply exactly to any other twodimensional (l-10) alkali halide withfl- = 0.907, 0%f = 0.846 and P_- = 1.15. They apply for any pair equilibrium distance r,, or ionic mass mf by suitable

reinterpretation of the time scale. It is interesting to ‘. note that the unit of velocity

and the unit of kinetic energy and temperature am?a(q f ) ‘, is independent.of r. (a resuit that is not true x three dimensions). This means that the interpretation of temperature of the simulation depends only on the charge per ion, in particular it is ind,-.pendent of the equilibrium ion separation and mass. Hence we can conclude that the characteristic temperatures, that is ‘to say melting-point, boiling-point and critical tempera

ttires, are the same for all sir&-ionised (1-l 0) +kzli halides.

two-dimensional

Other properties that do depend on p and pare not expected to vary much with p’” and /3- since these repulsive terms are quite small (I= 0.2%) for ions separated by e So as they are in the solid or liquid state.

4. The initial

lattice

It was found in a series of trial runs that a regular square lattice, representing a two-dimensional solid, does not remain stable (even when quite cold and no heat is supphed) if two adjacent ions are separated by the @r equilibrium distance r,,. The !at;rice has a strong initial tendency to expand and the momentum generated is sulfftcient to burst the crystal along the X and Y directions passing through the mid point. Fig. la shows a c,rystaJ’at 400 K with lattice spa&g r,-, and fig. 1b the system 6 x 1O-‘2 s later. The lattice seems to consist 502

‘,

..

1 October 1975

‘,

the stresses created are too strong to sustain. We show this result for two reasons, firstly because it depicts an interesting result of relevance to elastic behaviour of the system, and secondly because it compelled us to choose the inter-ionic distance to be 1 Jr,-, to obtain a stable initial lattice. This is shown in fig. 1 c.

5. The 1024 ion survey experiment The 1024 ions constituting the above stable lattice are initially given randomly oriented velocities according to a maxwellian velocity distribution such that the solid attains a temperature of = 400 K. The lattice was allowed to reach thermal equilibrium by letting the ions interact with one another without any input of energy from an external source for a period of 6 x 10S1” s. This time corresponds roughly to 60 vibration periods of a lone pair of unlike ions. The two-dimensional sample was then heated by multiplying all particle velocities by a factor, o (actually 1 .OOOS or 1.OOl), at every time step. This corresponds to an input of energy of (a2 - I)kT per atcm per step, or, if the heat input is slow enough to be reversible, to a constant increment of entropy ((u2 - 1)k per atom per step. The instantaneous temperature and the specific area were measured at every step. The instantaneous temperature, which does not necessarily imply thermal equilibrium, is defined as the kinetic energy per atom divided by Boltzmann’s constant, k. The specific area is obtained by associating an area of (3 x 3) mesh cells with each particle (= 2.25 $)_ A cell of the (256 x 256) mesh covering the whole region of calculation, is.then considered to be occupied by the material if any part of it is occupied by an atom.

6. Melting Melting was detected by observing a flat average temperature in the heating curve (temperature versus time) over a period of 110 time steps (lO_‘” s) preceded and

Vdume

1 October

CHEMICAL PHYSICS LETIiRS

35, number 4

(al

L

L

197.5

(b)

r

Cd)

(c)

(e)

(f)

L

Fig. 1. Solid nnd liquid phzs

fracture,

(c) stable c~std,

in twc-dimensiond

XII.

K* ions SSGWIJ8s plus, CL- ions show-

(d) melting, (5) hot liquid. ad

zs dot. (2) Unrtablc

crystal,

@J

(f) boilinp.

SC3

Volume 55, nun&r 4

’

1 I)ctobeer 197.5

CL-LEMICAL PHYSICSLETTERS

followed by a constant rise in terrperature. -In order to observe this, it is-necessary to average the temperature over ten time steps, which corresponds to the dominant vibration period of the system. The measured values at tZ-& phase transition are given below with the values for three-dimensional sylvin in brackets: melting point = 945 5 5 K (1043 I<>;

latent heat, L, 4 104 f 10 k erg/atom (1600 k); area expansion on melting = A4//i = 1 .O -t 0.2% (x20%). We note that the melting point is lower in two dimenslons due to the weaker nature of the two-dimensional forces. Most significant however is the small value of the latent heat and area expansion in two dimensions. This is a measure of the geometrical similarity of the liquid and solid phases in two dimensions as has recently been pointed out by Finneyf [lo]. In fact Finney predicts that there is no discontinuity of area at the melting point. Our measurements are best fitted by a small change which is verified by the more precise measurements on the 6400 ion sample (see section 10). Fig. Idshows the sample at the end of the tiat section -of the heating curve. At this stage, one can see the beginning of the break-up of long-range order mainly as the presence o.f numerous isolated ,defects. However ,+here is very little visual distinction between the solid and liquid at the melting point.

7. Boiling As heating proceeds (now with c’ = 1.00 1) there is rise in the sample and a gradu.al loss of long-range order. Fig. 1e shows the : liquid at .1280 K aUd.fig. 1fat 1400 K. At this stage disorder dominates in. the s&mple and ;I second flat is obseryed in the-heating curve which we associate with a boiling point: a steady expansion and temperature

.boiling point + 14OO.s 10 K (1650 K). Further heating causes complete vapourisation of .the sample, first to an assembly of long chain molecules, ~&. 2c and finally to ah assembly m:ainly consisting of KC! molecules, fig. 2d, fil!ing the entire area of the ex.- periment. The nominal instantqeons temperatures in ..the’se twdfigures are 4280 K and16000 K but the simu,. ~ .’

lation was not run for long enough, without the input of heat, in either case to ensure thermal equilibrium. It is clear however that long chain and loop structures are very stable and characteristic of the cold twa-dimensional vapour. It would be interesting to know whether there was any evidence for such structures in actual threedimensional materials. Further evidence for a boiling transition at 1400 K was obtained when the sample was heated past this value to 1500 K and then allowed to remain unheated for 6 x AO-” s. As thermal equilibrium was reached the temperature dropped to = 1420 K and expansion took place. Many internal ionic bonds were broken and holes appear throughout the sample. This transition is seen between figs. If and fig. 2a. Thus the latter figure may be considered as partially vapourised material.

8. Recrystallisation Fig. 2b shows the effect of cooling the partially vapourised liquid, fig. 2a, for 1.2 x lo-l1 s with a = O-999. The foal temperature recorded is 670 K in the solid region and several clear crystal domains are visible.

9. Properties of solid and liquid During the heating of the solid and liquid, mcasurements were made of the specific heats and expansion coefficients: SOLID (600-945

K).

Area expansivity at m.p., u=&A/AAT=

1.1 iO.l

x10-‘K-l;

specific area, A = 1.41 [l + 1.3 x lO~l.(K)]r~’

cm2/atom;

density: p = (0.707/ri)[l

- 1.3.x lOA qK>] atom/cm’;

linear separation between atoms, rr= 1.19[1 +0.65~lO~“~(K)]r,, specific hkat,

cm;

CHEMICAL

Valume 35, number 4

PHYSicS

1 October 1975

LETTERS

------I

.‘..I

..::

J

(Ill

CrvstIIs

. :.

_,

..

-.-~

--___

--.,:-.

K)

0.75 + 0.1 x :0-4 K;

p = (0.675/r;)

[i - 0.81

a = 1.22 [I + 0.40 x 10-4

“?

x

cm?‘/atom;

LOm4 T(K)] qK)]rO

‘_ ‘.

plus,

shown a~ dot.

Cl-ions

p(solid) = I..99 (1 - 1.05 x lo-$ TJ dcm3 ; P(liquid) = 1.97 (I - 2.94 x lo* T) g/cm3,

cp = 1.7 -C0.3 k e&/atom K.

A = 1.48 [l + 0.81 Y.1O-4 T(K)]r;

.,

(d) hot “apour

Fig. 2. Crystallization and vapourization in two-dimensional KCI. Kt ions showq as

01 =

‘_

-.. _.

(cl cold vapour

LIQUID (945-1400

_._, .._ ..-7

atom/cm

2. ,

but we notice a measurably smaller expansion coeffcient in the two-dimensional liquid compared with the solid. The specific heat is in nice agreement with the two-dimension21 form of the Dulong and Petit law which

predicts

a mmimum

of c

P

= 2k erg/atom

K.

cm;

= 1.7 t 0.3 k erg/atom K.

The’hensity a:., ,xpansion properties are of approximately the same magnitude as the three-dimensional values;

IO. The 6400 ion experiment

The survey experiment on the~lO24 ion sample, described above, was performed under non-equilibrium

Yglume 35, number 4

CHEMICAL

1 October

PHYSICS LETTERS

1975

‘I 850

950

900

1000

T’K

?jg

llrr!llrlllllll 850

950

900

moo

T°K Fig. 3. (a) Variation of internal energy, U, with temperature, T, at the fust order “solid” to “liquid” trcnsition in the 6400 ion system. (b) Variation of specific area, A, with temperature at the first order ‘~soiid” to “liquid” trznsition in the 6400 ion system.

conditions of constant heating. The errors quoted are a true estimate of the uncertainty of the measurements made from the data so obtained. In order to survey all the phase transitions in a reasonable computer time, a heating rate equivalent to about a IIS increase of the kinetic energy per shortest lattice vWation period was used. The adv&tage of this method is that continuous measurements are available of the thermodynamic variables. However we did not know to what extent the’resuits would be changed if the heating rate were much slower or if the systeln were allowed to reach equili-. brium. In addition some 12% of the ions in the 1024 sample are surface ions which obviously depart from .,bulk behaviour, and we donot know how sensitive the quoted resultS are 10 surface effects.

5

Fig. 4. (a) Structure of the “solid” phase at 906 K (point A of fig. 3). (b) Structure of the “Liquid” phase at 977 K (point B of fig. 3).

To answer both ,these &estions we have studied the solid--liquid transition with an Xl x 80 srrmple of 6400 ions, for which the equiiibrium separation at 800 K was found to be 1.25 ‘n. The surhcd atoms have now been

Volume 35, number 4

CHEMICAL

PHYSICS

reduced to 5% of the sample and the experiment was conducted under as near equilibrium as possible, as follows. The thermodynamic variables of a sequence of equilibrium states was obtained at intervals of 100 time steps. For the first 25 steps of the intct-Ial a heat pulse was injected by constant heating with 01= 1.00025. For the remaining 75 steps the sample was allowed to equilibrate and the values of the internal energy and temperature of the thermodynamic state were obtained by averaging over the last 50 steps. This corresponds to an average rate of heating one-ci&th of that of the survey experiment. Fig. 3a shows the variation of internal enera with temperature during the solid-liquid transition. The open circles denote measurements at each equilibrium state and a clear first order phase transition is seen at 912 K with a latent heat of 140 k per ion. The measured specific heat adjacent to the transition is 2.4 k/ion K on both the solid and the liquid side. Fig. 3b shows the variation of specific area versus temperature. A fist-order change is observed at 940 K with a relative area change of 2.7%. The results are in excellent agreement with the value of 945 f 5 K obtained for the transition temperature in the survey experiment. The qualitative result that both the latent heat and area expansion at the transition are about 10% of the three-dimensional values remains unchanged but it is clear that the survey experiment underestimated the actual values. Figs. 4a and 4b show the structure of the “solid” and “liquid” phases at the points A and B respectively ot fig. 3. The so-called liquid differs from the solid only in the presence of a greater density of defects and va-

LETTERS

1 October

197.5

cancies. There exists a clear long-range crystal structure throughout the “liquid” and its radial distribution func-

tion shows many more peaks than are characteristic of three-dimensional liquid. It is an open question whether

a

the

two

further

phases will respond differently to shear and experiments are planned to find out.

Acknowledgement The authors wish to thank the staff of the Rutherford High Energy Laboratory and the Atlas Computer L.aboratory where the calculations were performed, and the SRC and UKAEA for their support.

References S.P. Gocl and R.W. Hackney, Rev&n Brosileiva dc Fisim 4 (1974) 121. 121 A. Rahman, R.H. Fowler and A.H. Nartcn, I. Chcm. Fhyr

ill

57 (1572)

3010.

Chcm Phys. Letters 10 (1971) 257. I31 L.V. Woodcock, Cdmputer Phys. Commun. S (1973) 17. 141 P. Schofield, S.P. Soeland J.W. Ezstwood, Chcm. Phys. ISI R.W. Hackney,

Letters 21 (1973) 589. S.P. Coel and J.W. Eastwood, 161 R.W. Hockncy, [‘I

tional Phy-s 14 (1974) 148. J.W. Eastwood and R.W. Hackney,

J. Computa-

J. Computational

Phys

16 (1974) 342. J.W. Eastwood, Optimal Particle-Mesh Algoritilms, J. Computational Phys. (1975), to bc published. [9] R.N. Hockncy, Phys. Fluids9 (1966) 1826. [lo] J.L. Finney, Nature 242 (D73) 398.

i81