Geometry questions

Q. How can I check whether my geometry optimization really converged?

A. The first step is to read the <seed>.castep file and look for any errors or warnings. Has CASTEP complained about mixing symmetry and constraints? Or warned that there is nothing to optimize? Or warned that it has run out of search directions? There will be some explanatory text with each of these, to tell you what to do about it. Or has the SCF convergence failed for a particular structure? In which case, there is no ground state wavefunction or forces and so the geometry optimizer has to stop.

If there are none of these issues, then you can look to see if it has worked and what (locally optimal) structure it has found. The initial input configuration corresponds to 'step 0' and at the end of calculating the energy and forces (and stress if a variable cell optimization is being performed) you will see a block of text like this in your <seed>.castep file:

 LBFGS: finished iteration     0 with enthalpy= -1.30119012E+003 eV

 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS
 | Parameter |      value      |    tolerance    |    units   | OK? | <-- LBFGS
 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS
 |  dE/ion   |   0.000000E+000 |   2.000000E-005 |         eV | No  | <-- LBFGS
 |  |F|max   |   9.757539E-001 |   5.000000E-002 |       eV/A | No  | <-- LBFGS
 |  |dR|max  |   0.000000E+000 |   1.000000E-003 |          A | No  | <-- LBFGS
 |   Smax    |   4.056629E+000 |   1.000000E-001 |        GPa | No  | <-- LBFGS
 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS

This gives you a concise summary of the main convergence parameters at the end of step zero. The geometry optimizer will stop if all the criteria are met (all Yes in 4th column) OR if it has exceeded the maximum number of search steps (default = 30 but can be changed by the geom_max_iter parameter in the <seed>.param file). The convergence parameters have default values but can also be user-set in the <seed>.param file.

As the calculation progresses, the different geometry optimization methods have different algorithms for computing the next structure. Some also do a 'line search' which might mean calculating several closely related structures in a given 'step' to find the optimal 'distance' to move for that step. If this happens, you may see a block of text like this in your <seed>.castep file:

 +------------+-------------+-------------+-----------------+ <-- min LBFGS
 |    Step    |   lambda    |'  |    enthalpy     | <-- min LBFGS
 +------------+-------------+-------------+-----------------+ <-- min LBFGS
 |  previous  |    0.000000 |    0.038230 |    -1301.295179 | <-- min LBFGS
 | trial step |    1.000000 |    0.037008 |    -1301.295746 | <-- min LBFGS
 |  line step |   31.286152 |    0.015359 |    -1301.294832 | <-- min LBFGS
 |  quad step |   52.772393 |    0.004920 |    -1301.284019 | <-- min LBFGS
 +------------+-------------+-------------+-----------------+ <-- min LBFGS

which shows the minimizer (in this example geom_method=LBFGS) trying different step lengths 'lambda' in order to find the optimal step. This can be seen by the value of '' approaching zero. For a fixed-cell optimization this will also correspond to the lowest value of enthalpy in this step.

Finally, after a number of steps, you should see something like this in your <seed>.castep file:

 +------------+-------------+-------------+-----------------+ <-- min LBFGS
 |    Step    |   lambda    |'  |    enthalpy     | <-- min LBFGS
 +------------+-------------+-------------+-----------------+ <-- min LBFGS
 |  previous  |    0.000000 |    0.001152 |    -1301.956728 | <-- min LBFGS
 | trial step |    1.000000 |    0.000115 |    -1301.956779 | <-- min LBFGS
 +------------+-------------+-------------+-----------------+ <-- min LBFGS

 LBFGS: finished iteration    12 with enthalpy= -1.30195678E+003 eV

 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS
 | Parameter |      value      |    tolerance    |    units   | OK? | <-- LBFGS
 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS
 |  dE/ion   |   3.146456E-006 |   2.000000E-005 |         eV | Yes | <-- LBFGS
 |  |F|max   |   1.075030E-002 |   5.000000E-002 |       eV/A | Yes | <-- LBFGS
 |  |dR|max  |   7.645730E-004 |   1.000000E-003 |          A | Yes | <-- LBFGS
 |   Smax    |   4.612923E-004 |   1.000000E-001 |        GPa | Yes | <-- LBFGS
 +-----------+-----------------+-----------------+------------+-----+ <-- LBFGS

 LBFGS: Geometry optimization completed successfully.

which tells you that CASTEP was able to successfully geometry optimize this structure.

It may be that your structure was optimizing nicely, but CASTEP stopped due to exceeding the number of allowed steps. If that is the case, you will see something like this in your <seed>.castep file:

LBFGS: WARNING - Geometry optimization failed to converge after 30 steps

If, having checked the convergence summaries (e.g. use a grep command) then you might want to finish the calculation by allowing some more steps. You can do this by adding the following lines yo your <seed>.param file:

geom_max_iter : 50
continuation : default

which will restart the calculation from the last geometry step, and do (at most) another 50 steps.

Q. My BFGS geometry optimization did not converge. What can I do?

A. The usual reason is that either you ran out of steps and a few more will get it to work (see above for how to add more steps and/or change the stopping criteria) or that the forces are not accurate enough to find a good 'downhill' direction.

You can reduce the noise / improve the precision of your forces by a tighter SCF convergence. You can do this by making either

elec_energy_tol and/or

smaller in your <seed>.param file.

You can improve the quality of your forces by improving the quality of your wavefunction, by increasing the cut-off energy, the fine grid size, the number of k-points, etc.

Q. Which geometry optimization method should I use?

A. CASTEP supports a number of different optimization algorithms, which can be chosen via the geom_method keyword in your <seed>.param file. The allowed values are:

  • BFGS
    This was the default value until CASTEP v18. It uses the fractional atomic positions and cell strains as the variables to optimize, and looks for a state of minimum enthalpy (fixed cell) or zero force and stress (variable cell). It allows both variable and fixed cell calculations, with/without linear constraints. BFGS is reasonably fast and robust, but can require a lot of memory for large systems.
    This is the default value from CASTEP v18. It is the low-memory version of BFGS.
    This uses delocalized coordinates in a modified BFGS algorithm. It can use non-linear constraints but only works for a fixed cell calculation. It can be more efficient than normal BFGS/LBFGS for a large system where the initial configuration is close to optimum (e.g. derived from good experimental data). It searches for a state of minimum enthalpy.
    This is a form of optimally damped molecular dynamics. It is low-memory and reasonably robust, and can work when BFGS-type methods fail, but is generally less efficient than BFGS-type methods.
  • FIRE
    This is a molecular dynamics based optimizer, that has an equation of motion that includes local and global force and velocity components, but only works for a fixed cell calculation. It is low-memory and reasonably robust, and can be competitive with BFGS-type methods.
  • TPSD
    This is a two-point steepest descent optimizer. Unlike BFGS-type methods it does not require a line search. It is very robust and low-memory. In some systems (e.g. those that do not require a line search) it is almost as efficient as BFGS-type methods. It is the best choice when doing constrained-cell calculations.

Q. Should I trust the "Estimated frequency and bulk modulus" printed by a geometry optimisation?

A. Maybe.
These numbers are output at the end of a BFGS or LBFGS calculation, and are calculated by analysis of the Hessian matrix that has been accumulated during the optimization. If this Hessian is reasonably complete, then the numbers may be quite accurate. If the Hessian is incomplete (e.g. CASTEP converged in only a few steps) or has been updated lots (e.g. CASTEP required very many steps and traversed a number of different potential energy valleys) then the numbers can be inaccurate. The intention is that these numbers can be used to more efficiently initialize the Hessian for future geometry optimizations on this or similar systems, by setting the value of


in your <seed>.param file.