For each local structure, the "flatness" is define as the minimal mean squared distance of the neighbor atoms to the optimized plane fixed on the central atom. If the position of the central atoms is (px,py,pz), the normal vector of the plane is (A,B,C), then the function of the plan is:

A(x-px) + B(y-py) + C(z-pz) = 0

If the positions of the N neighbor atoms are {(x[i],y[i],z[i]) | i=1...N}, the flatness is

F = sqrt( sum( (A(x[i]-px) + B(y[i]-py) + C(z[i]-pz) )² ) / (A²+B²+C²) / N )

Represent the normal vector by the spherical coordinates (r, theta, phi):

A = r * sin(theta) * cos(phi)
B = r * sin(theta) * sin(phi)
C = r * cos(theta)

The flatness can be rewritten as

F= sqrt( sum( ( sin(theta) * cos(phi) * (x[i]-px) + sin(theta) * sin(phi) * (y[i]-py) + cos(theta) * (z[i]-pz) )² ) / N )

To find the global minimum of F is a two-dimensional optimization problem, so it can be solved by search all of the (theta,phi) space. The space has been devided to (100*100) cells and the minimal value of F and associated (theta,phi) have been found. Downhill Simplex method (Numerical Recipes in C, second edition, W.H.Press et al, Page 411) has been employed to further minimize (F, (theta,phi)) to the local minimum in terms of the precision of F² is better than 10^-15.

The variable space has also been devided to (1000*1000) cells but the value of minimized F is lowered only about 10^-4 angstrome. So (100*100) cells are considered enough to locate the position of the global minimum.

The histogram of F for different temperature levels are drawn below. More data are needed to reduce the statistical fluctuations. (One configuration only has about 280 surface atoms, thus the same amount of surface local structures.)

1) 100K

2) 200K

3) 300K

4) 400K

5) 500K

6) 600K

7) 700K

8) 800K

9) 900K

10) 1000K