DL-FIND

1. Introduction

DL-FIND is a powerful and flexible geometry optimisation library and is the recommended optimiser module in ChemShell. The code is under active development but the features documented here are those that are available in the version distributed with ChemShell. DL-FIND can also be downloaded as a standalone library from the CCPForge repository.

DL-FIND supports a wide variety of optimisation methods. Standard local minimisation techniques are supported (steepest descent, conjugate gradient, Newton-Raphson, BFGS) as well as damped dynamics and limited memory BFGS (L-BFGS). Transition state methods include nudged elastic band (NEB) for finding reaction paths and P-RFO and the dimer method for refining transition states. Conical intersection optimisations are supported using penalty function, gradient projection and Lagrange-Newton algorithms. Parallel optimisation methods include a genetic algorithm, stochastic search and parallel NEB.

DL-FIND accepts input structures in Cartesian coordinates. These may be converted into a number of internal coordinate systems for the optimisation, including delocalised internal coordinates (DLCs), total connection, and hybrid delocalised internal coordinates (HDLCs) for large systems. Constrained coordinates are supported as well as a general mechanism for imposing restraints.

A detailed description of DL-FIND's features can be found in:

Johannes Kästner, Joanne M. Carr, Thomas W. Keal, Walter Thiel, Adrian Wander and Paul Sherwood, J. Phys. Chem. A, 2009, 113 (43), 11856-11865. DOI: 10.1021/jp9028968

Please cite this paper if you use DL-FIND in published work.

2. Command line arguments

2.1 General options

2.2 Coordinate system options

• cartesian: Cartesian coordinates
• mass: Mass-weighted Cartesian coordinates
• dlc: Delocalised internal coordinates. A redundant set of bonds, angles, torsions, and sometimes inversions is created. A non-redundant combination of them is found by diagonalising the spectroscopic G-matrix. The optimisation is performed in this non-redundant set.
• tc: Same as dlc, but the redundant set consists only of bonds—all atoms in the system are connected (total connection scheme).
• hdlc: Hybrid delocalised internal coordinates. The system is partitioned into fragments. Delocalised coordinates (as in DLC) are used within each fragment. The fragments are coupled via Cartesian coordinates. This version is recommended for large systems, as the coordinate transformation scales linearly with the number of atoms.
• hdlc-tc: Hybrid delocalised internal coordinates based on the total connection scheme.

Note that the dlc option is not equivalent to specifying hdlc with a single residue. The dlc case gives pure internal delocalised coordinates (6 external degrees of freedom are dropped), whereas with hdlc Cartesian coordinates are added to the primitive list even if there is only one residue, and the external degrees of freedom are retained.

2.3 Optimisation algorithm options

• lbfgs: Limited-memory Broyden–Fletcher–Goldfarb–Shanno optimisation. Recommended for minimum search, NEB, and the dimer method. The time and memory requirements spent for determining the search direction and step length scale linearly with the system size. The number of steps kept in memory is specified by lbfgs_mem.
• prfo: The partitioned rational function optimisation method. Recommended for transition state search (unless done by the dimer or NEB methods). It requires the calculation of the Hessian. The keyword update_method can be used to specify how the Hessian should be updated.
• cg: Conjugate gradient method following Polak–Ribiere.
• sd: Steepest descent.
• dyn : Damped molecular dynamics. Four additional parameters can be specified. timestep is the time step in atomic units. fric0 is the start friction, fricfac is the factor to reduce the friction each time the energy decreases, and fricp is the friction to apply if the energy increases. The frictions are defined so that 0 corresponds to free (undamped) dynamics, and 1 corresponds to steepest descent. This variable friction facilitates convergence to an energy minimum.
• nr: Basic Newton-Raphson/quasi-Newton optimiser.
• ln: Lagrange-Newton constrained optimiser (for conical intersection searches only).

The trust radius approach (or line search) is controlled via the keyword trustradius. Possible choices are:

• constant: Constant trust radius (equal to maxstep)
• energy: Trust radius based on the energy change
• gradient: Trust radius based on the projection of the gradient on the current step. Can be used when no energy expression is available (e.g. NEB).

For optimisers requiring a Hessian, such as the Newton-Raphson or P-RFO algorithms, the type of Hessian to be calculated is specified by the keyword initial_hessian. Possible choices are:

• external: Calculate with the ChemShell routine $theory.hessian (default)
• onepoint: Build by one point finite difference of the gradient
• twopoint: Build by two point finite difference of the gradient
• diagonal: Build a diagonal Hessian by a single one point finite difference
• identity: Set the Hessian to be an identity matrix

In the case of initial_hessian=external, two point finite difference will be used if an external Hessian is not available.

The Hessian will be updated depending on the keyword update_method. Possible choices for update_method are:

• bofill: Bofill update of the Hessian (default)
• powell: Powell update of the Hessian
• bfgs: BFGS update of the Hessian. (Note: this is a full BFGS Hessian update, not to be confused with the limited-memory version).
• none: No update, recalculate the Hessian in each step

2.4 Nudged elastic band options

The NEB method may be run in parallel using the task-farming framework.

2.5 Dimer method options

2.6 Conical intersection options

The states to be optimised are chosen using the options state_i and state_j. Specific algorithm options can be set using pf_c1 and pf_c2 for the penalty function method, gp_c3 and gp_c4 for the gradient projection method, and ln_t1 and ln_t2 for the Lagrange-Newton method. The default values are those recommended in the MNDO implementation of these algorithms.

By default the gradient projection and Lagrange-Newton methods require the calling program to provide the interstate coupling gradient. However, if the interstate coupling gradient is zero by definition (for example, for states of different multiplicities), the algorithms can be used with the coupled_states option set to false, in which case only the two state gradients are required. The penalty function method never requires the interstate coupling gradient, but is usually less efficient than the other two methods.

Note that the pf and gp conical intersection optimisation algorithms are independent of the choice of optimiser in principle (in practice nr is recommended with bfgs updating). The Lagrange-Newton method can only be used with the Lagrange-Newton optimiser.

2.7 Population optimisation options

The genetic algorithm and stochastic search methods may be run in parallel using the task-farming framework.

3. Microiterative optimisation

4. Hessian calculations

The one-point or two-point finite-difference Hessian may be calculated in parallel using the task-farming framework.

5. Rate Calculations

5.1 Rates Calculated by Harmonic Transition State Theory

The Hessians at both the reactant state (RS) minimum and the transition state (TS) saddle point geometries have to be calculated by setting thermal=true. This writes the files qts_reactant.txt and qts_hessian_rs.txt in case of the reactant, and qts_ts.txt and qts_hessian_ts.txt in case of the transition state. In case of a TS, the crossover temperature T_c is calculated:

• Classical rates
• Rates with quantised vibrations (which includes the zero-point vibrational energy)
• Quantised vibrations plus tunnelling correction via the simplified Wigner correction to second order:

  

  The rate with quantum vibrations is multiplied by the expression above.
• Exact analytical quantum rates for a symmetric Eckart barrier fitted to the particular system (barrier hight and ). All degrees of freedom perpendicular to the reaction coordinate are approximated as quantum mechanical harmonic oscillators.
• For temperatures above the crossover temperature T_c, the full Wigner-corrected rates is also given:

  

KIEs can be calculated directly by first running DL-FIND with rate=true using the Hessians for the light isotopologue. The file arrhenius of this run can be copied to the file rate_H in the directory where the rate with heavier isotopes is to be calculated. There, the same class.in as in the light case (at least the same temperature parameters) should be used. The rates obtained with the light nuclids are read and the KIE is directly calculated and written to a file called kie.

5.2 Rates by Instanton Theory

Optimisation of the first instanton starting from the classical transition state structure (qts=true): The file qts_hessian_ts.txt has to be renamed to qts_hessian.txt. All geometrical data are read in from qts_hessian.txt. However, coords and coords2 still have to be provided (for historic reasons, number of atoms, ...), but are ignored (as in all instanton optimisations and rate calculations). A finite value of distort specifies how far the images will be spread along the unstable mode of the classical TS, see [1]. Newton–Raphson optimisation (optimiser=NR) is recommended [1, 2]. The number of images of the instanton path is specified by nimage. Instanton search is activated by qts=true. The NR optimiser was modified to avoid convergence to higher-order saddle points [2]. This avoids the collapse of the instanton path to the classical TS.

Instanton searches are performed in mass-weighted coordinates with masses consistent with atomic units (electron mass, m_e). That is, the mass of a hydrogen atom (¹H) is 1837.15 m_e. This scales all distances up by a factor of 42.695 (=(atomic mass unit/m_e)^1/2) compared to mass-scaled coordinates. Thus, the tolerance criterion (tolerance) has to be smaller by the same factor to achieve equivalent convergence. A tolerance of 10⁻⁷ (input as 1.E-7 in ChemShell) is usually sufficient, but convergence even to a tolerance of 10⁻⁸ can often still be achieved, which results in more reliable rates. Since NR converges quadratically, the more stringent tolerance generally does not increase the number of steps dramatically.

If NR (or P-RFO) is used, the updated Hessian will be used to calculate a preliminary estimate of the rate (if qts_reactant.txt is available). In that case, qts_hessian_upd.txt will be written, which contains only the updated Hessian. The file qts_coords.txt will be written. It acts as input for subsequent recalculation of the Hessian and a rate calculation.

Restarting instanton searches: Proper restart information (check files) is not written for the time being. Using NR, a restart is possible, though, by renaming qts_hessian_intermediate.txt (which is written after each step) to qts_hessian.txt and starting the simulation again. It will start from the Hessian and the geometry after the last full set of energies has been obtained.

Instanton rate calculation (qtsrate=true): qts_coords.txt from a previous instanton optimisation is read (coords and coords2 are ignored). The temperature is also read from qts_coords.txt. The rate calculation is selected by qtsrate=true. Hessians at all images and the rate are calculated as described in [2]. qts_hessian.txt is written, which acts as input for subsequent instanton optimisations.

Restarting of rate calculations is also only possible by using the Hessian information written for each image in qts_hessian_imageX.txt. For these files to be read, set inithessian=6.

Next instanton optimisation in sequential cooling: Starting from qts_hessian.txt at a previous (in general higher) temperature, another instanton is calculated. Distort should be zero, all other parameters are the same as above. The number of images may be increased. For optimal interpolation, the number of new images P_n should be related to the number of old images P_o by P_n = k P_o − k + 1 with k>1 being an integer. This ensures k−1 new images between each pair of old images.

Instanton KIEs can be calculated by starting out from a Hessian (qts_hessian.txt) for a different isotopologue and changing the masses in the input. The Hessian will be re-weighted accordingly. The instanton geometry has to be re-optimised. The file qts_reactant.txt obtained with changed masses can not be used. Instead, a file qts_hessian_rs.txt can be provided. The Hessian of the reactant obtained from that file will also be re-weighted.

In an approximation called fixed-path approximation (FPA) one can keep the instanton geometry fixed and just change the masses [3]. This is done by calculating an instanton rate with inithessian=5 (read the Hessian from file rather than recalculating it) and changing the masses.

For lower temperature (compared to T_c) the number of images necessary can be kept at bay by adapting the integration grid to the potential energy along the instanton path [2]. This only makes sense if the instanton path has reached the reactant minimum. The adaptive grid can be achieved by setting nebk=1 (This is not interpreted as the NEB force constant, but as a parameter which can vary from zero to one. One corresponds to a fully adaptive grid).

An instanton rate can be calculated from an existing Hessian (i.e. from the file qts_hessian.txt) without any additional energy calculations by setting inithessian=5.

Hessians of the individual images can be read in (all or just a part) from files qts_hessian_imageX.txt) by setting inithessian=6.

5.3 Bimolecular Rates

5.4 Tunnelling Splittings

5.5 References Covering the Implementation of Tunnelling Rate Methods

• [1] Judith B. Rommel, T. P. M. Goumans, and Johannes Kästner, J. Chem. Theory Comput., 2011, 7, 690. DOI: 10.1021/ct100658y
• [2] Judith B. Rommel, Johannes Kästner J. Chem. Phys., 2011, 134, 184107. DOI: 10.1063/1.3587240
• [3] Jan Meisner, Judith B. Rommel and Johannes Kästner J. Comput. Chem. 2011, 32, 3456. DOI: 10.1002/jcc.21930

6. Examples

• Minimisation in internal coordinates (DLC):

  dl-find coords=c \
      theory= mndo: $mndo_args \
      coordinates=dlc \
      result=minimum.c
  

• Transition-state search using the P-RFO algorithm (equivalient to "Baker" in newopt):

  dl-find coords=c \
      theory= mndo: $mndo_args \
      optimiser=prfo trustradius=const \
      delta=0.01 update_method=bofill maxupdate=50 \
      result=ts_prfo.c result2=tsmode_prfo.c tsrelative=true 
  

• Minimisation with internal constraints. Note that internal constraints (bonds, angles, ...) are only available in DLC and HDLC coordinates. Atoms can be fixed (frozen) in all coordinate systems using the active_atoms or frozen options.

  dl-find coords=c \
      theory= mndo: $mndo_args \
      coordinates=dlc \
      constraints= { { bond 1 2 } { torsion 1 2 3 4 } } 
  

• Nudged-elastic band: It is recommended to use a larger tolerance in NEB calculations than in normal TS search. The TS can subsequently be refined with the dimer method.

   
  dl-find coords=rs.c coords2= { ts.c ps.c } \
      theory= mndo: $mndo_args \
      coordinates=cartesian \
      optimiser=lbfgs trustradius=const \
      neb=frozen nimage=10 nebk=0.01 \
      tolerance=0.0045 \
      maxcycle=400 maxene=3000 \
      result=ts_neb.c restult2=tsmode_neb.c \
      maxstep=0.9
  

• Dimer method for subsequent refinement of the NEB transition state:

  dl-find coords=ts_neb.c coords2=tsmode_neb.c \
      theory= mndo: $mndo_args \
      coordinates=cartesian \
      optimiser=lbfgs tolerance=0.00045 trustradius=const \
      dimer=true delta=0.01 \
      maxcycle=400 maxene=300 \
      result=ts_dimer.c result2=tsmode_dimer.c \
      maxstep=0.5
  

  The transition mode can be visualised by writing a sequence in xyz format:

  tsmode_xyz coords=ts_dimer.c coords2=tsmode_dimer.c \
      n=20 delta=1. file=tsmode_dimer.xyz
  

• To minimise the system starting from a transition state, use the keyword distort. Note that any coordinate system can be used (e.g. coordinates=dlc) even though the transition mode has been calculated in cartesians.

  dl-find coords=ts_dimer.c coords2=tsmode_dimer.c distort=0.1 \
      theory= mndo: $mndo_args \
      coordinates=dlc \
      optimiser=lbfgs trustradius=energy \
      result=minimum.c
  

• Conical intersection optimisation using the Lagrange-Newton algorithm:

  dl-find coords=c \
      theory= mndo: $mndo_args \
      coordinates=cartesian \
      conint=ln state_i=1 state_j=2 \
      optimiser=ln tolerance=0.00045 trustradius=const \
      initial_hessian=diagonal update_method=bfgs maxupdate=1000 \
      maxcycle=400 maxene=300 \
      result=conint.c \
      maxstep=0.1 minstep=0.0
  

  Note that minstep should be set to zero as the Hessian must always be updated to ensure that the Lagrange-Newton constraints are imposed. The ln optimiser must be used with BFGS updates recommended.

• A potential-energy-surface scan is described in the tutorial.