Skip to content

DL-Find Based Optimization

Select minimization as your objective in the input file with the 

run minimize

key-value pair. The default will be to use the DL-Find library1 to minimize the energy of the input structure. DL-Find can also be used for transition state searches. In that case, the "run" keyword should be set to "ts" (see below).

Choosing coordinates

The efficiency of optimization can often be improved by using suitable coordinates. There are three classes of coordinates supported by DL-Find:

  • Cartesian: These are the usual xyz coordinates. Optimization in Cartesian coordinates usually takes more iterations than internal coordinates, but it may also be more robust.
  • DLC: These are nonredundant internal coordinates. The set of coordinates is computed once at the beginning of the optimization and is not updated. Thus, when/if the coordinates become poor, the optimization might stall or fail. In that case, restarting the optimization with the lowest energy geometry from a previous iteration may succeed.
  • HDLC: These are similar to DLC coordinates, but allow for the possibility of multiple fragments in the molecular geometry. Internal coordinates are used for each fragment, and Cartesian coordinates to describe the position of each fragment are added to the coordinate list.

The DLC and HDLC coordinates can also be used in a "total-connection" scheme, where the primitive internal coordinates are the entire set of atom-atom distances in the molecular system. Our experience is that these total-connection DLC and HDLC coordinates rarely perform better than the usual DLC and HDLC coordinates, but it may be worth trying them if you are having trouble. In summary, there are five different options for the coordinate system used in DL-Find: Cartesian, DLC, HDLC, DLC-TC, and HDLC-TC. We usually start with DLC and switch to Cartesian if there is any trouble.

Constraining degrees of freedom

It is often desired to optimize a molecular geometry while constraining certain degrees of freedom, for example in "dihedral scans" that are often performed for fitting empirical force fields. In such cases, some coordinate is specified as fixed and all other coordinates are relaxed to minimize the energy. DL-Find only allows coordinates that are in the original coordinate set to be constrained. For example, if one uses Cartesian coordinates, the only constraint which is possible is one (or more) of the Cartesian coordinates of one (or more) atoms. Hence, it is often the case that one uses DLC or HDLC coordinates for constrained optimization.

Input for constraints in DL-Find is formatted as follows:

$constraints
atom hydrogens
atom heavy
atom 1 5
atom 7
$end

where we show some example constraints. If Cartesian coordinates are being used, only "atom" constraints are valid. The lines in the example above will first constrain all hydrogen atoms, then constrain all heavy (non-hydrogen) atoms, and then constrain atoms 1, 2, 3, 4, and 5 (the first atom in the input geometry is considered as atom 1, i.e. we count from 1 here). Finally, we constrain atom 7 to be fixed. Note that this example is silly - we have constrained all atoms in the molecule! But it serves to show the possibilities.

The general format for a line in the $constraints section is

constraint_type atom1 atom2 atom3 atom4

where constraint_type can be atom, bond, angle, or dihedral. Possibilities for the atom constraint type were given in the example above. For bond, angle, and dihedral, TeraChem will expect 2, 3, or 4 atom indices (respectively). The atoms are expected to be listed in the conventional order. For water, the angle would list H1, O, and H2 atoms in that order. For ethylene, the dihedral about the CC bond would be listed as H1a, C1, C2, H2a, where H1a is bonded to C1 and H2a is bonded to C2. If Cartesian coordinates are selected, only the atom constraint type is valid. If any of the other four coordinate types are selected (DLC, HDLC, DLC-TC, HDLC-TC), only bond, angle, and dihedral constraint types are valid.

Output files for minimizations

The usual text output file (stdout) from TeraChem will provide information about the progress of the optimization. The geometries along the optimization path (with energies) are in the optim.xyz file found in the scratch directory. A list of the SCF, HOMO, and LUMO energies is found in the optlog.xls file in the scratch directory.

Worked examples

The minimal input for a DL-Find optimization is just run minimize:

opt_spiropyran.in
basis           6-31g
method          rhf
coordinates     sp.xyz
charge          0
run             minimize
timings         yes
end

If sp.xyz contains additional structures after the first, they are ignored during minimization (but used by NEB — see below).

Nudged elastic band (NEB) calculations

The NEB method is used to find a minimal energy path between two structures. Unlike the minimization methods, this is selected by

run ts

For NEB calculations, the input geometry must contain at least two structures (the starting and ending structures). Note that the path will try to connect atom 1 in the reactant with atom 1 in the product (and so on for each of the atoms). This means that it is very important to pay attention to how the atoms are ordered in the starting and ending structures. One can have more than two structures in the input geometry if one wants to specify an initial guess for the minimum energy path. Coordinates can be specified as described above. Beware that internal coordinates can be highly problematic for NEB calculations (the chosen internal coordinates must be well-chosen for both the initial and final structures - this is hard to accomplish in general). We thus usually use Cartesian coordinates for NEB calculations.

A minimal NEB input that searches for the path between two endpoints stored in the same coordinate file:

neb_spiropyran.in
basis           6-31g
method          rhf
coordinates     sp.xyz
charge          0
run             ts
nstep           200
end

Number of images and spring constants

The number of images along the NEB path defaults to 20 and can be changed with min_image. The default spring constant connecting adjacent images is 0.01 a.u. and is set with min_nebk:

min_image  16        # 16 images instead of the default 20
min_nebk   0.02

Variable spring constants help resolve the path near the barrier — higher spring constants are used for images near the highest-energy point and lower ones for images far from the barrier. To enable variable springs, set both min_nebk and max_nebk and ensure max_nebk > min_nebk (otherwise TeraChem reverts to the constant-spring scheme and prints a warning):

min_nebk   0.01
max_nebk   0.05      # higher constant near the top of the barrier

The number of initial guess vectors per image (used by the Davidson-like solver internally) is guess_per_image (default 2).

Force-modified NEB (FMNEB)

If a steering file is present, or the input file contains a steering keyword pointing at a steering file, the NEB calculation runs as force-modified NEB: the same constant pulling forces used in steered MD are added to every image of the NEB path. This is useful for studying mechanochemistry, where the barrier of interest is modulated by an external pulling force. FMNEB currently supports up to two attachment points; the details of the steering-file format are documented in the Steered MD page.

Output files for NEB calculations

  1. J. Kastner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander and P. Sherwood, DL-Find: An Open-Source Geometry Optimizer for Atomistic Simulations, J. Phys. Chem. A 113 11856-11865 (2009) ↩