Scalar Properties, Static Correlations and Order Parameters

A & T pgs. 46-58, 182-185, 199-204


Scalar Properties:


Today we discuss the basic properties we can calculate with simulation and also can be compared to experiment. What we can calculate with MD are microcanonical averages but to within factors of order (1/N) these should equal canonical averages. Later we will calculate the same things with Monte Carlo. The main scalar properties are:

It would also be possible to estimate the specific heat as the difference in energy with runs at two different temperatures. The disadvantage of this approach is that the statistical errors are amplified as the temperatures get close together and that separate runs are required.


Homework 3.1

Write the formula for the statistical error of the specific heat versus difference in the temperature dT. See detailed formulas.

    1. How does the systematic error depend on dT?
    2. Given a fixed amount of computer time, chose dT to make the ''statistical error'' = ''systematic error''.
      (Equating errors is legitimate, and easier, because we do not want one type of error to dominate).
    3. Finally how does this error (now actually dT) depend on computer time (or sample data points)?

It is not easy to calculate the free energy directly. We will discuss special methods to calculate it later.  Scalar quantities themselves are not usually very enlightening on the state of the system (e.g. liquid or solid). It is better to look at what is known as "order parameters"; this is one of the main advantages of computer simulation over experiment. We will discuss correlation function in real and k-space next time.
 

Static Correlations and Order Parameters

The most important information that comes from simulations are the static correlation functions. You should develop a habit of looking at these quantities routinely since they tell you basic information about what the particles are doing, namely where they are and how they arrange themselves.  The homework and exams require you to understand basic properties and relations of the static correlation functions. By static we mean properties that depend on properties computed with single "snapshot" or time step.  This is contrasted with "dynamical properties" they depend on how properties change in time.  We will talk about that next time. Most of the information from this talk is on the formula page .

The types of static correlations can generally be divided into

The density

The density tells you where the particles are (or have been). In a translationally invariant system it will be constant except for noise if you average long enough.

You find it by histogramming, that is subdividing the total volume into a set of subvolumes. How to histogram:

How do we set the histogramming size? By balancing the statistical error with the systematic error. The statistical error goes inversely as the square root of the volume size times the length of the run. It blows up for radial functions on a uniform mesh since the volume size goes at r2. The systematic error depends on the physical system, if the system is spread out, the density is pretty smooth, only a few bins will be needed.


Homework 3.2

Assume that you are computing the density in an inhomogenous system with N particles and M independent samples and the average density f(r) depends only on the radius r from the origin. Set up a radial grid with spacing dr. (This is a 3D system!)

    1. What is the estimator for f(r)?
    2. Assuming the particle postitions are uncorrelated, what is the statistical error of the estimate of f(r)?
      • Hint: Hits form a bitstream for a particular bin with NM uncorrelated bits.
      • Hint: Probability for any single particle ini a particular bin is estimator for p = f(r) dV/N, which is clearly correct for f(r)= constant = N/V.
      • Hint: Probability for any single trial hitting a bin is proportional to the size of the bin, i.e. p = dV/V. For the NM trials, the estimator for p = (# of hits)/NM.
    3. How does the systematic error depend on dr and f(r)?
      • Hint: How much does f(r) vary within 1 bin?
      • Hint: Assume that dr is small, i.e. dr << 1.
    4. What is the value of dr which minimizes the total error? How does it scale with M and r?
      • Again, rather than minimize total error, equate the systemmatic statistical error so that one does not dominate the error.

 

Three dimensional histograms are a problem because, first the number of hits/bin is low so noise level is high,  second they are hard to visualize, third they take alot of memory. You might want to use symmetry or projections toreduce a 3d function to 1 or 2 dimensional densities.

It is also good to examine the density in Fourier space. In periodic boundary conditions  the allowed K-vectors are on a grid. The Fourier density can be used to either: 1) smooth out the real space part by setting to zero the FT density larger than some value and fourier transforming back. (Disadvantage: the resulting density might go negative.) 2) look for periodic structures.

Pair correlations

The pair correlation function (also called the radial distribution function or g(r)) is the density of other particles around a given particle. Its first peak occurs around density -1/3. Normalization is usually defined such that at large r, g(r)=1. The theory of liquids is based on pair correlations because you can calculate the energy, pressure etc if you know g(r) (as long as the only potential is a pair potential.)  In a liquid or gas, g(r) depends only on the distance between two particles, not on the angle. Actually PBC break isotropy, g(r) could be different between the box axis and the body diagonal. Usually one ignores this to obtain a spherically symmetric g(r) for r less than half the box side. The pair correlation function is needed to apply 'tail' corrections to the potential energy and pressure coming about because of the truncation at the edge of the box.

The Fourier transform of g(r) is called the static structure function Sk. It is important to calculate because it can be compared directly with experiment, either neutron or X-ray scattering can measure it. There are very expensive instruments at Argonne, Grenoble, Brookhaven that are dedicated to measuring Sk and its time dependent generalization. If the agreement is good it gives a lot of confidence concerning the potential model. But don't be carried away; agreement to 5% is typical if the potential is at all reasonable.

The first peak of Sk is around 6 density 1/3. Normalization is such that at large k: Sk=1. It can be computed in two different ways, either as a FT of g(r) or the direct way. The direct way is slower but better at low k-values or a reciprocal lattice vectors of a possible solid structure because there are not the systematic errors of ignoring anisotropic correlations and histogram effects. Also g(r) must be extended to large r to do the FT.

One should use the peak of Sk to decide is a system is solid. For a liquid Sk is smooth and order 1 everywhere. A solid shows a delta functions with peak high a fraction of the number of particles at the reciprocal lattice of the solid. Solid structure shows up in g(r) in a more subtle way: instead of a nice damped oscillations, there are bumps at the nearest neighbor, next-nearest neighbor, etc. (Where they are depends on the type of lattice). See the figures for illustrations of g(r) in liquid helium.

The structure factor also signals separation in a two phase region. An example is the formation of droplets such as occurs when the number of particles corresponds to neither a stable liquid of gas, but a mixture of liquid and gas. In that case the structure function will diverge at low wave vector. The value at zero (S0) equals N by definition (in the canonical ensemble) but this can be different than limit {k Þ 0} Sk which is proportional to the compressibility of the system.

Order parameters

A key concept from the theory of phase transitions is the idea of an "order parameter" which characterizes a phase. For the moment we will just explain how this works for the two most common phase transitions: liquid-gas and liquid-solid.

The liquid-gas transition: the order parameter is the density. a liquid and gas are identical except a liquid has higher density. At the critical point the difference disappears. There is a first order phase transition between liquid and solid so it is quite possible to do a simulation in the two phase region. What will happen? Liquid droplets will form assuming 1) the surface tension allows it an 2) the simulation is long enough. When this happens Sk will have a characteristic behavior.

The liquid-solid transition: A solid is characterized by a periodic density. For example the Fourier transform of the density will be very large for certain wavevectors (the reciprocal lattice vectors) CLAMPS allows you to input a particular set of vectors to monitor melting and freezing. (keyword READ_K file )
 


Calendar

© 2000 D. M. Ceperley

© 1999,2001 D.D. Johnson