Geometry Optimizations

Relaxing a structure to the nearest minimum with the internal ANCOPT optimizer.

Table of contents

1. The --opt runtype
2. Optimization levels
3. Optimization with frequencies
4. Optimizing a whole ensemble

The --opt runtype

--opt relaxes the input structure to the nearest local minimum. Since CREST 3.0 the optimizer is a standalone implementation and no longer a call to xtb: ANCOPT, a rational-function optimizer with a BFGS update, working in approximate normal coordinates rather than in Cartesians. Working in ANCs is what makes the optimization converge in a handful of cycles instead of dozens.

The effect is easiest to see on a large, floppy molecule, where a Cartesian-space optimizer spends most of its effort untangling coupled soft modes:

Taxol is a good stress test for an optimizer: a fused taxane core with several ester side chains and a flexible benzamide tail, so the soft torsional modes are strongly coupled to each other. ANCOPT reaches the minimum in roughly 240 steps where L-BFGS needs about 425, and it gets the energy down faster throughout, not just at the end. The RMS force is the noisier quantity in both cases, which is why the energy change alone is not a sufficient convergence criterion.

crest struc.xyz --opt --gfn2 -T 4

# Geometry optimization
input   = "struc.xyz"
runtype = "optimize"            # alias: "ancopt"
threads = 4

[calculation]
optlev = "normal"               # optimization threshold preset

[[calculation.level]]
method = "gfn2"

            │           GEOMETRY OPTIMIZATION SETUP           │
          ┝━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┥
          │   algorithm                     ANCOPT          │
          │   optimization level        normal (0)          │
          │   max. optcycles                   200          │
          │   ANC micro-cycles                  20          │
          │   degrees of freedom                30          │
          ├─────────────────────────────────────────────────┤
          │   RF solver                   davidson          │
          │   Hessian update                  bfgs          │
          │   write crestopt.log.xyz          True          │
          │   linear?                        False          │
          │   energy convergence        0.5000000E-05 Eh    │
          │   grad. convergence         0.1000000E-02 Eh/α  │
          │   maximium RF displ.         1.0000000          │
          │   Hlow (freq-cutoff)        0.1000000E-01       │
          │   Hmax (freq-cutoff)         5.0000000          │
          │   S6 in model hess.         20.0000000          │
          └─────────────────────────────────────────────────┘

                            [ ... ]

│                                 CYCLE    5                                 │
└────────────────────────────────────────────────────────────────────────────┘
 * total energy  :   -14.5585939 Eh     change ΔE    -0.3534498E-05 Eh
   gradient norm :     0.0004548 Eh/α   predicted    -0.3009033E-05 ( -14.87%)
   displ. norm   :     0.0080889 α      lambda       -0.7691226E-06
   maximum displ.:     0.0077877 α      in ANC's #1, #2, #3, ...
   converged δE/grad :  True / True

   *** GEOMETRY OPTIMIZATION CONVERGED AFTER 5 ITERATIONS ***

                            [ ... ]

 TOTAL ENERGY         -14.5585939330 Eh
 GRADIENT NORM          0.0004548395 Eh/a0

  Important files written by CREST during this run
 ------------------------------------------------------------------------
 crestopt.xyz                │ final optimized geometry
 crestopt.log.xyz            │ step-by-step optimization trajectory
 ------------------------------------------------------------------------

The setup block printed at the start lists everything the optimizer was configured with, including the two convergence criteria that have to be met simultaneously: the energy change between cycles and the gradient norm. Each cycle then reports both, together with the displacement, and flags them as converged δE/grad : True / True once satisfied.

For the 1-propanol example this takes five cycles and gains 0.69 kcal/mol over the input geometry.

Two files are written:

Optimization levels

How tightly the structure is relaxed is controlled by an optimization level, given directly after --opt on the command line or as optlev in the [calculation] block. The level scales both convergence thresholds at once:

Run on the same input, the levels differ in cost and in how far they push the energy down:

The spread between crude and vtight is 0.019 kcal/mol here, well below the accuracy of the underlying method. This is why the conformational search workflows optimize at loose levels for most of their intermediate steps and only tighten up at the end.

Optimization with frequencies

--ohess runs the optimization and follows it directly with a numerical Hessian, giving vibrational frequencies and a full thermochemical analysis at 298.15 K:

crest struc.xyz --ohess --gfn2 -T 4

          ::       THERMODYNAMICS at      298.15 K        ::
          ::::::::::::::::::::::::::::::::::::::::::::::::::
          :: TOTAL FREE ENERGY        -14.479960446075 Eh ::
          ::----------------------------------------------::
          :: total energy             -14.558594533089 Eh ::
          :: ZPVE                       0.106452687996 Eh ::
          :: G(RRHO) w/o ZPVE          -0.027818600982 Eh ::
          :: G(RRHO) total              0.078634087014 Eh ::
          ::::::::::::::::::::::::::::::::::::::::::::::::::

Low-frequency modes are treated with the modified (quasi-)RRHO approximation rather than the harmonic one, with the cutoff and entropy model reported above the table. In addition to the optimization output, this writes:

Optimizing a whole ensemble

To relax every structure of an existing ensemble instead of a single geometry, use the ensemble optimization runtype, which parallelizes over the structures:

crest --mdopt crest_conformers.xyz --gfn2 -T 8

This is documented in detail on the Ensemble Optimization page, which also covers --screen, the variant that additionally sorts and filters the result.