Fortran Analysis Structure
The cluster program is self contained. It reads in a configuration of atomic positions
and produces a circular linked list for each cluster identified within the configuration.
The best value of the critical separation r_cl depends on the particular physical system
being considered. To illustrate, we have provided a file cluster.inp consisting of
N=256 atoms at low overall density, generated by a short quench from a disordered system
at high temperature to a low temperature inhomogeneous state (not yet at equilibrium).
This file should be copied into the working directory before running the program.
With the default value r_cl=1.5, six well-separated clusters should be identified.
The results in this case are moderately insensitive to the value of r_cl, but increasing
it above 3 includes all atoms in a single cluster, while reducing it below 1.15 will start to
separate isolated atoms into clusters of their own.
Clustering algorithms are part of the standard toolkit of data analysis, and in practical
applications it may be more efficient and convenient to use a packaged implementation of
an algorithm such as dbscan
- M Ester, H-P Kriegel, J Sander, X Xu. (1996). Proc. Second Int. Conf. on Knowledge Discovery and Data Mining (KDD-96) p 226 (Eds: E Simoudis, J Han, UM Fayyad; AAAI Press, 1996).
Fortran implementations of dbscan are available from various sources.
For systems in periodic boundaries, rather than supplying the atomic positions, the user should
compute a distance matrix using the minimum image convention, and supply that to the routine,
as suggested by Turci
in the context of a Python version.
The program pair_distribution reads in a set of configurations and calculates
the pair correlation function g(r).
We limit the number of configurations to a maximum of 1000 (numbered from 000 to 999)
simply so as to use a fixed naming scheme for the input configurations;
in a practical application, a trajectory file would be used instead.
We have tested it on a set of 500 configurations
of N=256 Lennard-Jones atoms,
cut (but not shifted) at Rc=2.5σ,
at the usual state point ρ=0.75, T=1.0.
The interval between configurations was 100 MC sweeps.
This data set is provided in the
file pair_distribution_data.zip in the Data repository.
Using the default resolution of 0.02σ,
the results shown below were obtained for g(r).
The program grint.f90 reads in a set of configurations and calculates
the pair correlation function for a system that is inhomogeneous in the
z direction. It is assumed that the configurations consist of a liquid slab,
surrounded by gas, with the interfaces lying in the xy-plane.
The two interfaces are located by fitting the instantaneous density profile
to a difference of two tanh functions. Then the single-particle density function,
relative to each of the interface locations, is calculated.
To make the process as robust as possible, an initial guess at the mid-point
of the liquid slab should be provided, and this is updated automatically as
successive configurations are read in, so as to shift the liquid slab into
the middle of the periodic box, before fitting.
Also, the results of one fit are passed on as the starting point of the next one.
The program handles cubic boxes only;
the modifications necessary to handle non-cubic boxes are fairly easy to make,
but we will not do that here.
No attempt is made to correct for capillary-wave fluctuations affecting the
width of the interface.
Having located, and combined, the two interfaces, the histograms for calculating the one-body and two-body densities are accumulated, in a coordinate system which has its origin at the interface. The one-body density is then fitted by a single tanh curve: this could easily be changed by the user if desired. For simplicity, in the normalization of the two-body density to give the pair correlation function, we use the fitted tanh form of the single particle density. This is obviously an approximation, and could be replaced by a better fit, or an interpolation scheme, to evaluate the function at arbitrary z. Finally, the program writes out the average, overall, density profile, the single-particle density (with fit) and slices at selected values of z and c through the pair correlation function.
In testing this program, it is important to use a large enough system so that all values of z1 and z2 of interest (measured relative to the interface position) lie far from the other interface position.
The program was tested on a system of N=10000 atoms, interacting through
the Lennard-Jones potential cut (but not shifted) at Rc=2.5σ,
in a cubic box of side 30σ, at a temperature T=0.90.
For this system, ρG ≈ 0.024, ρL ≈ 0.713
(see A Trokhymchuk, J Alejandre, J Chem Phys, 111, 8510 (1999)). A set of 100 configurations from this run, together
with the output of grint.f90 with default parameters, are provided in the
file grint_data.zip in the Data repository.