Minimization of internal parameters (min_lapw)

Most of the more complicated structures have free internal structural parameters, which can either be taken from experiment or optimized using the calculated forces on the nuclei. 

 

Starting with WIEN2k_11.1 there are two possibilities to determine the equilibrium position of all individual atoms automatically (obeying the symmetry constraints of a certain space group). One can use either

• the shell script min_lapw, together with the program mini, which will run a scf-cycle, update the positions using the calculated forces and restarts a new scf cycle. This continues until forces drop below a certain value;
• or use the normal scf-scripts run_lapw with a special input in case.inm such that the charge density and the positions are simultaneously optimized during the scf-cycle.

At present we recommend the second (new) option, although there are cases where this scheme can be slower or may even fail to converge.

A typical sequence of commands for an optimization of the internal positions would look like:

• Generate struct file
• init_lapw
• run_lapw -fc 1 [another runXX script or additional options are of course also possible] (this may take some time)
• Inspect the scf file whether you have significant forces (usually at least .gt. 5 mRy/bohr), otherwise you are more or less at the optimal positions (An experienced user may omit the run_lapw step and proceed directly from init_lapw to the next step)

Now you have to decide which method to use:

• min_lapw [options] (this may take some time)
  • it will generate a default case.inM (if not present) by:
    • executing ``x pairhess -copy ; cp case.inM_st case.inM '' (i.e. it sets up the PORT minimization option and calculates an approximate starting Hessian).
    • when -nohess is specified, it will generate case.inM from SRC_templates with the NEW1 option (not recommended).
  • Without -NI switch min_lapw performs an initialization first:
    • removes "histories" (case.broyd*, case.tmpM) if present;
    • copies .min_hess to .minrestart (if present from previous min_lapw or x pairhess).
• or edit case.inm and put MSR1a (or MSEC1a) as ``mixing method''. Then continue with run_lapw -fc 0.5 -ec 0.0001 -cc 0.001 [-it]. It will run x pairhess (unless case.inM is already present) and then run (several hundreds) scf-cycles, simultaneously updating positions and charge densities. Once the forces seem to be smaller than the limit defined in case.inM it will switch to ``mixing method'' MSR1 and finalize the scf-cycle with fixed positions. Because of this, the final forces may not be as small as desired and eventually you have to restart this step using MSR1a again.

When using the second method we recommend you read carefully $WIENROOT/SRC_mixer/Mixer_README_4.1.pdf. Overall the method is very good for semiconductors (or well behaved metals), and allows ``tricks'' like small k-mesh or small RKMax at the beginning of the minimization and using higher accuracy only towards the end.

The following text refers (mainly) to the first method using min_lapw:

When case.scf is not present, an scf-cycle will be performed first, otherwise the corresponding forces are extracted into case.finM and the program mini generates a new case.struct with modified atomic positions. The previous step is saved under case_1/2/3.... Then a new scf-cycle is executed and this loop continues until convergence (default: forces below 2mRy/bohr) is reached.

The last iteration of each geometry step is appended to case.scf_mini, so that this file contains the complete history of the minimization and can be used to monitor the progress (grep :ENE *mini; or :FORxxx ...).

By default (unless switch -noex is specified), min will call the script clmextrapol_lapw after the first geometry step and try to extrapolate the charge density to the new positions. This procedure usually significantly reduces the number of scf-cycles and is thus highly recommended.

mini requires an input file case.inM (see Sec. 8.15) which is created automatically and MUST NOT be changed while min_lapw is running (except the force tolerance, which terminates the optimization).

We recommend the PORT minimization method, a reverse-communication trust-region Quasi-Newton method from the Port library, which seems to be stable, efficient and does not depend too much on the users input (DELTAs, see below with NEWT). The PORT option also uses/produces a file .min_hess, which contains the (approximate) Hessian matrix (lower-triangle Cholesky factor) If you restart a minimization with different k-points, RMT, RKmax, ... or do a similar calculation (eg. for a different volume, ...) it will be copied to .minrestart (unless -nohess is specified), so that you start with a reasonable approximation for the Hessian. The program pairhess, which calculates the first Hessian, also prints out the average Hessian eigenvalue for the symmetric, symmetry-preserving modes in mRyd/au2 as well as the minimum and maximum, and also the vibration frequencies. A list of these is given at the end of case.pairhess. Note that these are not all possible modes, but only the symmetry preserving ones. Therefore if you have prior information about the vibrations of the system you can adjust the rescaling term so the average vibration frequency is about correct. (see the description of pairhess in 9.2). (In addition there is a program eigenhess, which will analyze the Hessian after the minimization has been completed. It also prints vibrational frequencies and may give you hints about dynamical instability of your system. Some more description is given in $WIENROOT/SRC_pairhess/README and at the top of the output file case.outputeig.

When using PORT you may also want to check its progress using

grep :LABEL  case.outputM

where :LABEL is :ENE (should decrease), :GRAD (should also go down, but could sometimes also go up for some time as long as the energy still decreases), :MIN (provides a condensed summary of the progress), :WARN may indicate a problem), :DD (provides information about the step sizes and mode used). Some general explanations are:

1) The algorithm takes steps along what it considers are good directions (using some internal logic), provided that these steps are smaller than what is called the trust-region radius. After a good step (e.g. large energy decrease) it expands the trust-region; after a bad one it reduces it. Sometimes it will try too large a step then have to reduce it, so the energy does not always go down. You can see this by using ":DD'' and ``:MIN" .

2) A grep on :MIN gives a condensed progress output, in which the most significant terms are E (energy in some rescaled units), RELDF (last energy reduction), PRELDF (what the algorithm predicted for the step), RELDX (RMS change in positions in Angstroms) and NPRELDF (predicted change in next cycle). Near the solution RELDF and RELDX should both become small. However, sometimes you can have soft modes in your structure in which case RELDX will take a long time before it becomes small.

3) A warning that the step was reduced due to overlapping spheres if it happens only once (or twice) is not important; the algorithm tested too large a step. However, if it occurs many times it may indicate that the RMT's are too big.
4) A warning "CURVATURE CONDITION FAILED" indicates that you are still some distance from the minimum, and the Hessian is changing a lot. If you see many of these, it may be that the forces and energy are not consistent.

Sometimes PORT gets "stuck" (often because of inconsistencies of energy and forces due to insufficient scf convergence or a very non-harmonic potential energy surface). A good alternative is NEW1, which is a "sophisticated" steepest-descent method with optimized step size. It can be very efficient in certain cases, but can also be rather slow when the potential energy surface is rather flat in one, but steep in another direction (eg. a weakly bound molecule on a surface, but constraining the sensitive parameters, like the bond distance of the molecule, may help).

Another alternative is NEWT, where one must set proper "DELTAs" and a "FRICTION" for each atom. Unfortunately, these DELTAs determine crucially how the minimization performs. Too small values lead to many (unnecessary) "geometry steps", while too large DELTAs can even lead to divergence (and finally to a crash).
Thus you MUST control how the minimization performs. We recommend the following sequence after 2-3 geometry steps:

grep :ENE *mini
:ENE  : ********** TOTAL ENERGY IN Ry =        -2994.809124                    
:ENE  : ********** TOTAL ENERGY IN Ry =        -2994.813852                    
:ENE  : ********** TOTAL ENERGY IN Ry =        -2994.818538

Good, since the total energy is decreasing.

grep :FGL001 *mini
:FGL001:  1.ATOM                         0.000          0.000         18.219    
:FGL001:  1.ATOM                         0.000          0.000         12.375    
:FGL001:  1.ATOM                         0.000          0.000          7.876

Good, since the force (only a force along z is present here) is decreasing reasonably fast towards zero. You must check this for every atom in your structure.

When you detect oszillations or too small changes of the forces during geometry optimization, you will have to decrease/increase the DELTAs in case.inM and rm case.tmpM. (NOTE: You must not continue with modified DELTAs but keeping case.tmpM.) Alternatively, stop the minimization (touch .minstop and wait until the last step has finished), change case.inM and restart.

You can get help on its usage with:

min -h or min_lapw -h

PROGRAM:        min

USAGE:          min [OPTIONS]

OPTIONS:
-j JOB ->       job-file JOB (default: run_lapw -I -fc 1. -i 40 )
-noex ->        does not extrapolate the density for next geometry step
-p ->           adds -p (parallel) switch to run_lapw 
-it ->          adds -it (iterative diag.) switch to run_lapw 
-it1 ->         adds -it1 (it.diag. with recreating H_inv) switch to $job
-it2 ->         adds -it2 (it.diag. with reinitialization) switch to $job
-noHinv ->      adds -it -noHinv (it.diag. without H_inv) switch to $job
-sp ->          uses runsp_lapw instead of run_lapw
-nohess ->      removes .minrestart (initial Hessian) from previous minimization
-m ->           extract force-input and execute mini (without JOB) and exit
-mo ->          like -m but without copying of case.tmpM1 to case.tmpM
-h/-H ->        help
-NI ->          without initialization of minimization (eg. continue after a crash)
-i NUMBER ->    max. NUMBER (50) of structure changes
-s NUMBER ->    save_lapw after NUMBER of structure changes

CONTROL FILES:
.minstop        stop after next structure change

For instance for a spin-polarized case, which converges more difficultly, you would use:

min -j ``runsp_lapw -I -fc 1.0 -i 60''