Dispersion curves for nearest neighbour FCC lattice in the 100 CUBIC direction

IMPORTANT NOTE: The FCC basis was assumed to be a 4-atom cubic basis & the system considered was a cube with periodic boundaries. The dynamical matrices were only calculated for modes in the 100 direction.

Fig 1: The dispersion curves for all 12 modes.

Fig 2: Dispersion for acoustic mode.

Density of states for 2D slices of a 3D harmonic crystal

Fig 1(top) & 2(bottom): These demonstrate that the density of states scales like the 3rd power of frequency.


Fig 3: The 111 patch in this case crosses the system's periodic boundaries (it doesn't self intersect). The DOS still shows the same scaling as Figs 1 & 2.

Computing Green's function for a 2D slice of a 3D crystal

Crystal: FCC lattice with harmonic nearest neighbour interactions. The domain is assumed to be cubic with periodic boundaries.

FCC Basis: 4 atom Cubic FCC basis. The lattice vectors are simply aligned along the cubic axes.

Calculation Scheme: First the 12x12 Hessian matrices are written down as a function of CUBIC lattice separation. Note: since we have only nearest neighbour interactions, we have at most 27 lattice translations whose corresponding hessian matrices have non-zero entries.

The Dynamical matrices, D(k), are then computing for all CUBIC wavevectors that can fit in the simulation box via taking the DFT of Hessian matrices, that is: D(k) = SUM{{R} H(R)*exp(i k.R)}}, where (R} is the set of all possible translations in the lattice (NOTE: R and -R must be separately counted). This DFT was not done through the FFT since the number of non-zero Hessians is much smaller than the number of lattice translations.

Each dynamical matrix was diagonalized to get the 12 polarization vectors & their corresponding eigenvalues {(P_i, L_i) : 1<= i <=12} and the Green's function in fourier space is computed using: G(k) = SUM{(i=1...12) P_i.outerProduct.P_i/L_i }

Finally the real space Green's function can be obtained for the whole crystal through an inverse Fourier Transform: G(R) = SUM{{k} G(k)*exp(i k.R)} where {k} is the set of all fourier modes that can fit in the lattice.

To obtain the Green's function for the slice, I rotated the tensor obtained and discarded all components that had a 'z' piece. Finally I was left with a 2x2 matrix which gave the 2D correlations in a slice (for a particular separation vector) from a 3D crystal.


Difficulties: Due to the fact, that G(k) is the inverse of D(k), we have problems when computing G(k) for long wavelength fourier modes (since their translational and longitudinal polarizations; eigenvalues are very small, hence these D(k) have high condition number). Specifically we can have G(k) & G(-k) differ significantly enough for G(R) to be not purely real. I solved this problem by explicitly setting: G(k)=G(-k).

Density of States - Boson peak

Fig 1: Convergence of DOS of a 111 plane slice of a 3D FCC crystal (20% disorder width). Note the boson peak at around omega=6.9

Convergence of MSD of soft-sphere particles

In the figure above the x-axis is sampling time while the y-axis is the variance of the MSD of 111 plane computed at a particular timestep. The distribution of MSD of the ordered system must go to a delta function as we sample longer. The variance must go like ~ 1/N (where N=number of independent samples). The measured exponent of the power-law decay in the variance of the ordered system's MSD is ~= 1, in accord with expectation. The disordered system's MSD distribution must have an intrinsic variance which is exhibited by the plateauing of the variance.

FCC lattice: planewave projections of full lattice and of 111,100 plane slices

FCC lattice: Cubic Indexing and mapping to Bravais Index, handling PBCs

For constructing neighbour lists, while accounting for PBCs, we have applied a scheme of using 'cubic indexes' for each atom in our FCC lattice.

NOTE: Our simulation cell is cubic, but it could be cuboidal too & this indexing would still be good.

Cubic Index: Each atom is assigned a unique 4-tuple (m,n,p,s) of integers.

The 3-tuple (m,n,p) indexes the cubic Bravais cell. Each cubic cell has 4 atoms which are distinguished by their 's' index according to the following rule: (atom_position) --- ('s' index)
v0: (0,0,0) -- s=0
v1: (1/2,1/2,0) -- s=1
v2: (0,1/2,1/2) -- s=2
v3: (1/2,0,1/2) -- s=3
(v1,v2,v3 are the FCC basis vectors)

Each time a lattice is constructed with cubic indexing, a pickled dict: 'cubicIndexMap.pickle' is specifying a map from (m,n,p,s)-->ID is also written.

A map from the integer Bravais indices (a,b,c), which are integers specified along the primitive FCC basis vectors, can be specified to the cubic indices using only integer math.
F: (a,b,c)-->(m,n,p,s)
(To construct F we write the position of the atom first in Bravais space and then expand the basis vectors into 3D cartesian space, as specified above, noting that EITHER 2 vector components OR none can at most be half-integers - depending upon the value of s).

The map F: (a,b,c)-->(m,n,p,s) is,
m,v[1] = divmod(a+c,2)
n,v[2] = divmod(b+a,2)
p,v[3] = divmod(b+c,2)

s is assigned by checking if 'v' is v0,v1,v2 or v3.
NOTE: divmod () is a built-in python function which returns the quotient & the remainder of a division.

Applying PBCs:
When finding Bravais neighbours, one can quickly map the cubic index and then to the ID (using the pickled map). In case periodic boundaries need to be handled, we merely need to 'mod' cubic indices (m,n,p) with the number of cubic cells along a cartesian axis. NOTE: 's' should be unchanged.
See: the function getCubicIndex() in
/home/ahasan/usr/local/simulations/fcc-3d/disordered/lattice/pickFCCplane.py

Covariance Matrix of 111 plane of FCC, NVE solid -- Ordered

For a 32x32 parallelogram shaped region (1024 atoms) of a 111 FCC plane, I constructed the covariance matrix sampling 4000 snapshots separated by 10tau (the decorrelation time).

MATLAB eig() was used to diagonalize this 2048x2048 matrix for the eigenvectors and the eigenvalues. It took less than 5 mins.

Fig 1: is a histogram of 1/sqrt(covariance eigenvalues) for the ordered system.

Green's function for 111 plane

Using the cubic basis I computed lattice and time averaged Green's functions for atoms separated along the x-cartesian axis. For the ordered and the disordered lattice the Green tensor components looked qualitatively similar to Fig. 1. In Fig 2 & Fig 3, I show the decay exponent of the Green's function with increasing separation.

NOTE:
(i)The averaging was done over 4000 independant frames and NOT 4000 tau as indicated below.
(ii) I also computed the VARIANCE of each Green's function tensor component (NOT PLOTTED).


Fig 1: The green's function (showing all the 6 tensor components)



Fig 2: The green's function on a log scale showing its decay with increasing separation



Fig 3: For the DISORDERED lattice

Measuring decorrelation times for FCC 3D, soft spehere, NVE, T=0.005

I computed the auto-correlation of the projection of displacement fields on the longest wavelength plane wave for the 111 plane of my FCC, NVE ordered and disordered system.

NOTE: The autocorrelation was measured for the COMPLEX fourier coefficients and NOT the square amplitudes of the projection onto plane waves.


Based on the plots below, 10tau seems to be a reasonable time interval to sample statistically independent configurations.

Sorting .gz dumpfiles

Am now sorting .gz files and writing them to disk in unsorted format. All .gz files except for the first one will be deleted. I should also gzip the sorted dumpfile for additional space efficiency.

the file sorting program is sitting in:
/home/ahasan/usr/local/simulations/fcc-3d/disordered/sortDumps.py

New kinds of simulations: bonded, disordered, ordered ensemble

1) I ran a purely harmonic system with the neighbouring atoms bonded to each other. The configuration files were unfortunately deleted :(. However the setup files are on navier at: /home/ahasan/usr/local/simulations/fcc-3d/bondSimul/
-- The makeLattice.py file doesn't have cubic indexing yet :(
-- I wrote a module which computes the 3D neighbour list in O(n) time. It's sitting in this folder.

2) I ran a simulation of system where there were 10 different kinds of species with their individual interaction strengths evenly spaced. Each particle was assigned to some specie with equal probability. This simulation is sitting in:
/home/ahasan/usr/local/simulations/fcc-3d/disordered/

-- The lattice in this simulation has cubic-indexing.

3) Ran an ensemble of 5 ordered simulations starting the same initial FCC lattice and at the same temperature.This simulation is sitting in:
/home/ahasan/usr/local/simulations/fcc-3d/ordered/

-- The lattices in this simulation have cubic-indexing.

MSD of fcc, soft-sphere, gamma=0 runs

The LAMMPS runs whose MSDs are plotted below were discussed earlier in the following post:
http://dislocationnucleation.blogspot.com/search?updated-
max=2010-05-19T07%3A16%3A00-07%3A00&max-results=7

Both these runs have the same number of particles (~350,000), are NVE & run at different temperatures.

On the y-axis I have plotted the average value of the square of the displacement computed at a given timestep (x-axis).

T=0.01 shows very large movements from the initial position, while for T=0.005 the perparticle MSD plateaus.

Bond Length History in 2D dics

Fig 1: Bond length fluctuation from equilibrium separation plot for gamma=1.0 (colourbar on top)

NVT sample (Bond fluctuation history)

Temp:               0.0025
Disorder widths:    1.0, 0.5, 0.25, 0.125
timstep:            0.1
# obs:              2 million (2e6)
sampling frequency: 100

The data is located on navier at:
/home/ahasan/usr/local/simulations/disordered/half/bondHistory

Velocity participation factors in lowest 4 modes (2x triple)

Calculation of shear modulus for 2D, dense, bi-disperse, soft sphere granular mixture

On making an affine shear transformation of a square cell of granular, bi-disperse particles at equilibrium we get non-affine displacements to preserve equilibrium.

I divided the domain into many pieces. Below is an example of the domain divided into 4 pieces and their individual non-affine corrections plotted together for convenience.



In the following picture, the domain is divided into 400 pieces and the non-affine correction to the shear modulus is computed for each sub-domain and plotted.


NOTES: (i) these are computed using the formulas using the formulas in Craig's PRE of 2006. (ii) there were 2 critical bugs in the coding the potential & I had confused the diameter as the radius making the system super dense. (iii) the code that does this stuff is sitting on research-1 at:
/second_drive/asad/granularShearData/

Lattice averaged Green's function from the covariance matrix

Using the covariance matrix I computed the average value of the following components of the green's function for each reciprocal lattice vector (k): Grr, Gnn, Grn.
NOTE: (i) 'r' is along 'k' while 'n' is normal to 'r'.
(ii) name of the module: computeGreensFunction.py (in ..../disordered/2d-22/ on navier).

Plots of Grr, Gnn, Grn (respectively) follow:




NOTE: this program hasn't been thoroughly debugged.

2D hex, soft sphere, disordered system density of states

For my most disordered system (gamma=1.0), I assembled the covariance matrix for 4 million steps and diagonalized it for its eigenvalues. This system has 4200 particles. The plot above shows the convergence of 1/eigenvalues for various sampling intervals.

The plot below shows the convergence of histograms of the eigenvalues plotted above.



In the plot above, I have shown the DOS of various window sizes (shown in the legend). All eigenvalues were obtained by diagonalizing the covariance matrix obtained after sampling 4 million steps. For small frequencies, we expect linear DOS. Note the bimodal nature of the DOS. For comparison, below is a plot of the ordered (gamma=0.0) system's DOS.

100 plane, FCC 3D, mean square amplitude of displacement field FT onto plane waves

The plane waves are identical to those shown in the previous post.

The SLOPE of the linear part in the above plots was measured to be nearly -1.2 which was close to the expected 2*(2/3).

Fourier transfrom of 2D, disordered, soft-sphere, NVT displacement field

For my most disordered 2D system (gamma=1.0) I computed the average projection of displacement fields onto transverse plane wave traveling along the x direction. Essentially i computed the quantity:
U' = sum(n: 0-->N-1) {exp(2*pi*i*(n/N)*x)*Uy}

Here is an example of the real part of the plane wave with n=5 & N=60:The above is the 'cos' part, the 'sin' is merely shifted by half a fringe-width.

Now the following is a plot of , averaged over 4 million timesteps at T=0.0025.

The different curves are for different window sizes, the number in the legend denotes the number of particles along the horizontal axis (radius=1.0). On the vertical axis, the square amplitude is measured in units of particle size, while on the hor axis the wavevector is in unit on 1/particleSize. For comparison, above is the corresponding plot for the ordered system, gamma=0.0.

Melting in FCC, 3d, soft sphere crystals

I found that using LAMMPS' Nose-Hoover NVT setup the melting temperature was very low ( below 1E-5). To see this I simply tracked a particle for some time after starting the simulation. Why this is the case is unexplained especially since in 2D i have been getting reasonable statistics at T=0.0025.

Switching to NVE, after some short sampling at various T, it seemed that melting seemed occurred between T=0.05-0.1. I ran 2 simulations at T=0.01 and T=0.005. The following plots are the Y versus X coordinate of a SINGLE particle sampled at 10,000 timesteps from 0-4 million steps.
The MSD of the particle in the upper plot was HUGE ~27.5 !! The main contribution was ~= 26.2. It seems that T=0.005 is suitable for analysis.

NOTE: the good config (lower picture) is in the: fcc-3d/dumpfiles3/ folder and its 100 plane is extracted to the folder: fcc-3d/extractDump3/planeDump3/.

Work done from 20th Jan '10 till today

1) Projected velocity vectors for 2x-triple on lowest 4 eigenmodes. To compute the velocity vectors I found that I couldn't use the linearSover, MATLAB's cgs() since it was giving me unexpected results, however using gaussian elimination solved the problem.

2) The projections of the velocity vector on modes was in accord with expectations: the most active mode showed the highest participation in the velocity field.

3) I then started working on mesoScale analysis for 3x-triple & I performed the following calculations:
a) Made a lowest eig as a function of (x,y,R) picture for R={4,6,8,12,16} while choosing (x,y) from a circle of radius 20.
b) Centered the mesoscale region on the hottest guy and made pictures of the lowest mode. Noticed that the character of the lowest mode changed radically between 24

Ineffiencies in computing the numerical hessian

1) Maybe due to starting LAMMPS 2N times.
2) Maybe due to use of popen.
3) Could be due to heavy file I/O, sice that 's the only way to communicate with LAMMPS ...

Below are some timing results:

Legend:
totalTime - total time take for computing the hessian
timeNeib - time taken to compute the neighbour lists
timePerturb - total time spent inside perturb() per a FULL hessian calculation
timeLAMMPS - total time spent inside the popen() call which runs LAMMPS
fileTime - total time spent doing file I/O for atlking to LAMMPS
timeHess - time spent manipulating the hessian dict structure

totalTime timeNeib timePerturb timeLAMMPS fileTime timeHess
49.7423419952 0.434784173965 48.5797655582 36.8780109882 11.70175457 0.727792263031
63.6533300877 0.439907073975 62.4323446751 50.4112865925 12.0210580826 0.781078338623
70.9052629471 0.468246936798 69.3777544498 55.5753176212 13.8024368286 1.05926156044
72.0013990402 0.472586154938 70.4885210991 55.9716053009 14.5169157982 1.04029178619
75.1916220188 0.476332902908 73.6517179012 59.05403018 14.5976877213 1.06357121468
75.3300080299 0.479758024216 73.775904417 59.0307564735 14.7451479435 1.07434558868
75.5174300671 0.566131830215 73.8336679935 59.0396382809 14.7940297127 1.1176302433
75.4272179604 0.481420993805 73.8312358856 59.0277297497 14.8035061359 1.11456108093
75.4155731201 0.489336013794 73.8398208618 59.0288999081 14.8109209538 1.08641624451
75.5248250961 0.570980072021 73.8176853657 59.0196013451 14.7980840206 1.13615965843