Homogenous Dislocation Nucleation

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

2010-10-01T16:50:00.000-07:00

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

2010-10-01T16:41:00.000-07:00

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

2010-10-01T15:04:00.000-07:00

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

2010-09-17T09:14:00.001-07:00

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

2010-08-31T09:17:00.000-07:00

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

2010-08-26T15:22:00.000-07:00

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

2010-08-18T07:00:00.000-07:00

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

2010-08-17T15:52:00.001-07:00

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

2010-08-17T15:51:00.000-07:00

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

2010-08-17T15:35:00.001-07:00

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

2010-08-17T15:32:00.000-07:00

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

2010-08-17T15:16:00.000-07:00

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

2010-07-15T16:42:00.001-07:00

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

2010-07-01T11:01:00.000-07:00

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)

2010-05-21T12:46:00.001-07:00

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

2010-05-21T09:15:00.000-07:00

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

2010-05-20T23:37:00.000-07:00

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

2010-05-20T23:27:00.000-07:00

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

2010-05-20T08:20:00.000-07:00

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

2010-05-19T07:16:00.000-07:00

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

2010-05-18T13:35:00.000-07:00

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

2010-02-09T11:18:00.000-08:00

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

2010-01-19T12:19:00.000-08:00

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

Update: smooth indenter & analytical calculation of velocity & hessian

2009-12-20T20:11:00.000-08:00

Smooth Indenter: Wrote a new class in smoothIndenter.py which encapsulates all information about the smooth indenter. This class produces the correct LAMMPS script needed on "moving" the indenter. The file minimizeWithIndenter.py was suitably modified for this type of indenter.
On running the indenter I found that, x-quadruple gives the expected nucleation pattern BUT x-triple nucleates 2 dislocations even after moving the indenter (horizontally) some 5-6 times.

Analytical Velocity: This has the advantage that the indenter needn't be given a small perturbation to measure the velocity field. The equation used is in another post. In order to make solution of the linear system seamless and not involve any MATLAB use, I have decided to adopt the scipy sparse.dok_matrix format for all matrices and vectors. The newTheoHessian.py file has all the changes but is yet to be tested. More results on this will be posted here when available.
Fig1: 'Crystal response forces' are the analytically calculated forces in the crystal due to an infinitesimal downward motion of the indenter. The displacement field is computed from 2 configurations 1e-6 apart in indenter depth. Numerical velocity is computed by moving the indenter by 1e-7 and calculating the resulting crystal displacement.

Analytical Hessian: is finally working for LJ smooth (& LJ cut) potential. In fig 2. (below) are the lowest mode of x-triple. Also quoted below are the 4 lowest eigenvalues (calculate from full hessians using Lanczos). There is very close agreement between the analytical and numerical results.

4 lowest Eigenvalues of full hessian of x-triple :-
Numerical : -0.0058792, -0.24143, -0.59427, -0.77723
Analytical: 0.0058764, 0.24143, 0.59427, 0.77723

Fig 2: Lowest eigenmodes of x-triple.

Analytical velocity field and hessian for 2-body potentials

2009-12-16T15:12:00.000-08:00

Fig 1: 'Velocity' is defined as the rate of change of positions of crystal atoms with respect to downward motion of the indenter. The above expression assumes that the system is in mechanical equilibrium so that the response is essentially linear.

Fig 2: Derivation of elements of the Hessian matrix for a 2-body potential in a Cartesian coordinate system.

Evolution of lowest eigenvalues near dislocation nucleation

2009-12-08T07:54:00.000-08:00

Fig 1: Plot of 4 lowest eigenvalues for 4 system geometries. The 4 'lowest' eigenvalues of the full system are plotted against indenter depth.

Fig 2: Fitted functions of the form y=-A*sqrt(B-x) - the cyan curves -to the tails of the eigenvalue ~ Depth plots (where they approach zero); A & B are constants. Except for x-triple, all systems show the existence of 1 critical mode whose eigenvalue seems to go to zero as the square root of (Dc - D) & for all systems, the parameter B was the critical depth. (in accord with expectations from a 'saddle node bifurcation' based nucleation mechanism).

Fig 3: The 4 lowest modes from the system x-triple corresponding to the last configuration before nucleation (roughly del_D ~ 1e-6). This image shows that more than 1 normal mode is active during nucleation partially accounting for different character of xtriple's eigenvalue plots. Mode 1 corresponds to the lowest eigenvalue (the one closest to zero).
Fig 4: The 4 lowest modes for the system x-quadruple corresponding a configuration with (Dc-D) ~ 1e-6. Notice the strong contrast with the modes of x-triple. Mode 1 is again the lowest.

Plotting Info: Used my numerical hessian routine to compute negative-definite hessian matrices. An ARPACK based Lanczos implementation built into matlab was used to diagonalize the hessian matrices. Only configurations with max_force <= 1e-8 were for the computation of hessians. And the minimizations were done with the pre-MR minimizer, ie. backtrack_quadratic.

Structure of paper for unsccm'09 procs

2009-11-22T12:53:00.000-08:00

1x - triple, quadruple; 2x- triple, quadruple

delta(s) curves; delta_dot fields.

Standard notation summarized elsewhere, to mirror the proposal as closely as possible to avoid ambiguities.

Particular system: 2xtriple

Questions: Where shall we describe our indentation
Defect Structure(fig refs from proposal)
1) Saddle node scaling - deltaDot(depth)
--Fig 7 left
2) Rescaling of deltaDot wrt sqrt(delD)
--Fig 4 (left & right)
Story for 1 &2)
-> Usual saddle node square root criticality arising from (atleast) a 3rd order polynomial representing energy as it can form inflection/saddle points. Locus of minima is a diverging parabola.
-> Why we work at 0K: want to track the minima as closely as possible
-> Can demonstrate: evidence of the square root scaling, onset of scaling from as far as 1e-2 : should try to relate this to energy change in the crystal as this shall be a measure of barrier height?
-> Geometrical Instability: role of lowest mode(s) leads smoothly to point 4 (so switch 3 & 4?).

3) Small values of delta(s) @ criticality
--Fig 6 (all 3) + Picture of the velocity vector field without circles {in a grid of 2-by-2}
**Talk about the Pierels framework -- look harder at the proposal.

4) Lowest eigs wrt to deltaD (have this only for MD, do we need for energy minimization???)
--Standard 2 pictures (eigVals evolution + zoomed view) from the Chaos submission.
** Usual story from the Chaos submission. May have to modify if we decide to have energy minimization modes.

Mesoscale Resolver:
1) Min_eig (x,r,R;delD)
Have picture for 2xtriple.
--Fig 8 (left), maybe fig 8(centre) -- not sure as it doesn't present a clear picture.
--Festering Wound Picture from the blog (15 Oct, '09 post).
Have log scale plot but nothing to show dependence on delD.
-- Might be interesting to demonstrate pictorially the observation: 'as long as the mesoscale probing region intersects the defect core, we obtain the characteristic dislocating lowest mode'. What sort of picture shall one make here ? A 2-by-4 grid might do a good job ?? Just found the pictures, they look somewhat strange ...(in Mac in ~/Desktop/mesoModes/)

2) delta_dot(s) (x,r,R;delD)
--Will have to make picture, not entirely clear what we want. What does it mean to look at meso-scale velocity/displacement field?

___________________________________________

More to do on indentation:
-- Use previous displacement as an initial guess with decimation and see what happens.
-- Writing class smoothIndenter to abstract fix_indent. Main problem: efficiently (& cleanly) changing a single line in a text file.

Update: running the force-zeroing minimizer - bactracking, problems...

2009-11-17T10:27:00.000-08:00

The energy minimizer (adapted from MR) tries to find force zeros. However when the minimizer tries to take a large step, we need to take a smaller step (backtrack). Two strategies are currently applied for backtracking: (i) Decimation - new step is 1/10 of previous step, (ii) AH - linearizes force and solves for the zero of the straight line.

Decimation is significantly slower than AH because it uses no information about the magnitude of the forces & hence tries a large number of moves whereas AH zeroes onto the correct move very quicky. Since decimation samples a large number of points of the search space, most of which are far from the starting point, the chance of it finding a minimum far away from the staring point is very high.

The latter case has been observed in indentation runs where decimation frequently shows a load drop before AH, especially when stepping by 1e-6 & sometimes when we step by 1e-5 too.

Problems/Observations with indentation runs:
(i) The biggest problem is that 4x-quadruple nucleates along 2 modes. Despite shifting the indenter horizontally, the problem persists. The best I was able to do was with AH, where upto steps of 1e-4, I was able to obtain nucleation along a single mode. However I still got a featureless velocity field.
ASIDE: An MD run, with fix viscous on, showed only one nucleation mode. The indeter depth was ~1 particle spacing lower than energy minimization.

(ii) I just observed that one gets a featureless velocity field in 2xquad if the state had been reached with decim but we try to define velocity with AH (or vice versa, can't recall). This is very peculiar and should be looked into.

(iii) All of the other systems gave defect-incipient velocity vector field that allowed computeCrossStrain.py to do its job 'fromScratch', although most velocity fields weren't 'clean' - especiallly those from 'decim' (probably because of its tendency to jump far).

A new line search algorithm (adapted from Miller & Rodney)

2009-10-27T08:00:00.000-07:00

Explanation of Specific Points:

Choosing initial alpha: We want to avoid moving more than dmax in the entire linesearch
and we want to make sure that atleast 1 move is made. The initial alpha ensures two things: the first move (alpha_init*hmaxall) will not be greater than dmax & the first move is not too small (since alpha_init is either 1.0 or 0.5*dmax/hmaxall).
Allowing an energy increase smaller than 1e-12: First of all, a move that increases the energy (by 1e-12) is allowed iff the gradient condition is immediately satisfied, else that move is undone. Now a small energy increase is permitted because changes in energy below some small tolerance ( I chose 1e-12) is not significant because we are working in finite precision (the number of particles is ~1e4 and we have about 14 significant figures for the energy of individual particles). The choice of 1e-12 is therefore some bound on measurable energy change of the whole system (the specific choice is a LAMMPS default value).
Rubberbanding: Notice that the term F_(i)/{F_(i-1) - F_(i)} = 1/(relative change in directional derivative). When linesearch routine reaches this step it means: that we just made a move which decreased the energy and the directional derivative but the directional derivative wasn't sufficiently decreased, so we want to continue moving in the search direction. Rubberbanding says: linearize the directional derivative and solve the resulting linear equation for its zero, this gives the 'z'. The bound of 4.0 is put on the boost factor('z') to avoid making alpha too large, as that'd run the risk of attempting to step outside the allowed search space.
Backtracking: Here we do 2 things:- undo the move made in the current iteration and compute a new value of alpha (which must be smaller than the previous alpha). Decreasing alpha by a factor of 10 was a conservative choice made by me. The MR expressions for decreasing alpha were too complicated. We could also try to linearize force at this point or fit it to a parabola of the form: f = (alpha - a)^2 + b (a,b are the parameters) and choose the new step by solving for the zero of the force equation.

More on Using mesoScale probes

2009-10-15T07:23:00.001-07:00

Figure 1: Atoms coloured according to the lowest vibrational eigenvalues of mesoscale regions of various radii (given by the 'size') centred at themselves. Note the sharp transition in pictures of size 6 & greater.

How I generated these plots:
For system 2X-triple I took a region of radius = 20.0, at each atom belonging to this region a mesoscale region (of radius indicated next to each picture) was centered.
For each mesoscale region I plotted the lowest vibrational eigenvalue on the atom which was its centre.
NOTE: the hessian is negative-definite so more negative eigenvalue means more stable.
Figure 2: Lowest eigenmodes of mesoscale regions of radius 8.0 with centers displaced parallel to the critical plane with respect to the center of the defect core. The sense of arrows is arbitrary.

Probing with mesoscale regions

2009-10-12T11:55:00.000-07:00

1) Procedure: For system 2X-triple I took a region of radius = 8.0 and moved it perpendicular and parallel to the slip from the hottest atom. For each mesoscale region I plotted the modes and delta fields and computed the lowest eigenvalue. No mesoscale region hit the surface of the crystal.

2) The plot at the top shows the meso scale regions' lowest eigenvalues. The transition at around 5-6 in the trend of eigenvalues was most interesting. A transition in the qualitative nature of the corresponding eigenmode of the mesoscale happens at the same time, which is reflected in the delta field plots.

3) I made plots of the delta curves on the slip plane for each of these mesoscale regions too. They are consistent with the information from the full delta fields.

4) Issue of bulk eigenvalue -?

News At This Point: Transliteratin Fortran code

2009-10-06T07:12:00.000-07:00

Spent the last week working on figures for Craig's proposal. The proposal is a good summary of work done to date, except it doesn't say anything about work on the minimizer.

Now I am starting work on transliterating Ron Miller's F77 program's linesearch to C++ such that it can be seamlessly used with LAMMPS.
The relevant file is: mod_solve.f and the subroutine is cgstep().

The central idea in Ron Millers linesearch:
1)In each linesearch try to achieve:
(i) abs(final gradient)/abs(init gradient) <= 0.1, OR
(ii) detect the point where the gradient just changes sign
NOTE: gradient here means directional derivative.

2) The energy is only used to avoid increasing the function value.

3) If this direction doesn't do well for us then request a new direction, many times reset CG by using the bare gradient as the search direction.

*** The minimizer section is not yet complete ***

A bad situation for Armijo backtracking & Next Steps

2009-09-12T15:56:00.000-07:00

In the configuration on the right, we start from a position in the energy landscape where the energy is very steep.

---
Next Steps:
Try running molecular dynamics with zero viscosity and try both indenter types. Use the 6000 atom system.
Analyze new LAMMPS line_minimizer (7Sep version).

Results of running plain Armijo bactracking

2009-09-10T21:23:00.000-07:00

I disabled the quadratic branch in the linesearch, which forces backtracking by cutting alpha by a factor of 2. This new scheme performs much worse as I get failed linesearches in systems where the indenter has barely interacted with the lattice. Besides, there are failed linesearches on each step of the indenter :(

On inspecting the output from the linemin, it became clear that the failed linesearch arose due to |de_ideal| being less than IDEAL_TOL. I think that this condition is caused because we are very close to the minimum to start with. In this specific instance f.s didn't even change sign (it remained negative).
Note: All the above observations were made with IDEAL_TOL = 1e-10

Next I tried to minimize the first indentation configuration with IDEAL_TOL = 1e-18. With this new IDEAL_TOL the backtracker performed much better but I still got a failed linesearch. Despite f.s changing sign, the backtracker failed to identify the energy minimum because the Armijo sufficient decrease condition wasn't satisfied. Eventually, because of cutting alpha by a factor of 2, a failed linesearch arose when: |de_ideal| < IDEAL_TOL.

Conclusion: backtracker needs more muscle !

Min::linemin_quadratic() flowchart

2009-09-08T16:02:00.000-07:00

Concerns:
(1) We may be abandoning plain Armijo backtracking when it may be having some remaining mileage . The danger in this is that the quadratic mode may satisfy Armijo by jumping to the initial point.

(2) The quadratic mode assumes success when:
(de <= de_ideal || de_ideal >= -IDEAL_TOL) == True
At the initial point this is automatically true (even when forces aren't zero).
Note that: de_ideal = -0.4(alpha)(initial_f.s)
and in quadratic mode: alpha = alpha0

(3) Jumping to the initial point can cause the minimizer to go into an infinite loop causing it to exit with either MAX_FORCE_EVALUATIONS or MAX_ITERATIONS (ie. max cg loops).

(4) Example of linemin_quadratic() when it exit with success by jumping to the initial point:

----
count fdotdirall energy alpha start
0 2.436356006811224e-03 -2.719796543157099e+00 3.460098962684194e-01 initial
1 -1.801321567525278e+01 -1.684889377507936e+00 3.460098962684194e-01
2 -1.017649600012320e+00 -2.658368603360779e+00 1.730049481342097e-01
3 -2.847044105111986e-01 -2.708236347904628e+00 8.650247406710485e-02
4 -1.597636733426530e-01 -2.715590284316164e+00 5.290154968720939e-02 insideQuad
5 -1.297466257039866e-01 -2.716985938194178e+00 4.325123703355242e-02
6 -6.561635732195255e-02 -2.719095785009185e+00 2.162561851677621e-02
7 -3.280817866097375e-02 -2.719629373969908e+00 1.081280925838811e-02
8 2.436356006811224e-03 -2.719796543157099e+00 1.665334536937735e-15 insideQuad
9 2.436356006811224e-03 -2.719796543157099e+00 1.665334536937735e-15
9 2.436356006811223e-03 -2.719796543157099e+00 1.660130366509804e-15 0.000000000000000e+00 -1.617867436214352e-18 exitingLinemin
---
The 3rd last entry in line 9 (last line) is the decrease in energy compared to the initial point & the 2nd last entry is de_ideal.

This caused an infinite loop and at the end the minimizer quit due to MAX_FORCE_EVALS, leaving the final state with large forces.

----

Lammps: Min::linemin_quadratic()

2009-08-31T23:40:00.000-07:00

Some observations on linemin_quadratic():

1) 'dir' -- passed as argument by the CG routine isn't a unit vector. In the 1-D instance it's the force. There is a check for: f.s > 0 ('f' is force and 's' search direction), the linemin() immediately exits with error if this isn't satisfied.

2) normflag = 0, hence 'fdotdirall /= atom->natoms' doesn't have any effect. Further, nextra seems to be 0 too as fdotdirall += fextra[i]*hextra[i] doesn't have any effect. The only non zero component of the 9-component global force is vector is f[3] (this is the x-component of atom 2's force).

3) initAlpha = min(1.0,(dmax/hmax) ; dmax is a supplied parameter & is = 0.3 in my case ; hmax = biggest component of search direction in absolute terms.
lammps checks -- if (hmax==0): linemin() exits with error
So for the 1-D case, initAlpha is always 1.0

4) At the start of each linemin(), the calling CG routine passes the current config, current energy & forces. Assume that these initial quantities' consistency isn't suspect.

5) All steps in linemin() are made ONLY with respect to the initial configuration, ie. x_new = x_init + alpha_currrent*dir;

6) In my file (linemin) I'm dumping: count,f.s,f[3],energy,x[3],x[6],alpha
f[3],x[3] - are force and x-position of particle 2, x[6] is x-position of particle 3
In some cases a last string field is there to indicate where in the code is the dump coming from.

Energy minimization

2009-08-28T11:01:00.000-07:00