Update: smooth indenter & analytical calculation of velocity & hessian

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

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

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

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...

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)

Explanation of Specific Points:

1. 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).
2. 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).
3. 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.
4. 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

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

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

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

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

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

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()

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