Geometry Optimisation

DL-FIND is a powerful and flexible geometry optimisation library and is the standard 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 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.

Running an optimisation

class Opt

In ChemShell, DL-FIND optimisations can be run using the Opt() driver, which expects a Theory object, e.g.:

my_theory = NWChem(...)
my_job = Opt(theory=my_theory)

The geometry optimisation calculation is run with:

my_job.run()

Then the results will be stored in my_job.result (see: <class Result>).

Methods

Opt.run(dryrun=False)

Parameters:: dryrun (bool, optional) – If the function is called with dryrun=True, the job will only generate input scripts without executing real calculatons

General options

DL-FIND optimisations take the following arguments:

theory: Choose the theory level by specifying a Theory instance

frag2: Initial second configuration of the system. End point of the initial path in NEB calculations, direction of the dimer in dimer calculations. Can also be a list of fragments, in which case it defines points along the initial NEB path (the desired endpoint should be given last).

maxcycle: (default: 100) Maximum number of optimisation cycles

maxene: (default: 10000) Maximum number of energy evaluations

restraints: A list of geometry-dependent Restraint instances that are to be applied to the system.

save_path: (default: False) Instruct ChemShell to save DL-FIND output files regardless of debug level.

Note

Of particular interest is the output file _dl_find/path.xyz, which shows the complete trajectory of the optimisation

Coordinate system options

coordinates

Specifies which coordinate system to use

Allowed values:

'cartesian': (default) Cartesian coordinates
'mass': Mass-weighted Cartesian coordinates
'dlc': Delocalised internal coordinates (DLC). 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 bond - 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

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.

residues

Defines residues (atom groups) for use with HDLC coordinates.

Accepts a dictionary of the form:

{ 'res1' : [atom1, atom2, atom3, ...],
  'res2' : [atom4, atom5, atom6, ...] }

..note ::: For biomolecular Fragments this dictionary can be generated from PDB data using the getHDLCResidues() method on the fragment.

constraints

Define constraints for use with either HDLC or Cartesian coordinates.

For HDLC coordinates, internal constraints can be defined using a list of the form:

[ ('bond', atom1, atom2),
  ('angle', atom1, atom2, atom3),
  ('torsion', atom1, atom2, atom3, atom4),
 ... ]

For Cartesian coordinates, constraints can be defined using a list of the form:

[ ('x', atom1), ('yz', atom2), ... ]

where the first constraint freezes the X coordinate of atom1, the second freezes the Y and Z coordinates of atom2, etc.

active

List of atoms to relax during the optimisation. All other atoms will be held fixed.

e.g. active=[0,1] will relax the first two atoms in the fragment and keep all other atoms fixed for the optimisation.

By default (with no active option specfied), all atoms will be relaxed during the optimisation.

Note

Py-ChemShell counts atom labels from zero, so the first atom in the system is atom 0.

frozen

Inverse of active. List of atoms to hold fixed during the optimisation. All other atoms will be allowed to move during the optimisation.

e.g. frozen=[0,1] will fix the first two atoms in the fragment in place for the optimisation.

Optimisation algorithm options

algorithm

Specifies the optimisation algorithm

Allowed values:

'lbfgs': (default) Limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) 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).

trust_radius

Specifies the trust radius approach

Allowed values:

'constant': (default) 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).

maxstep: (Default: 0.5) Maximum change of one coordinate component in one optimisation step.

tolerance: (Default: 0.00045) Convergence criterion for the maximum gradient component. The criterion for the maximum step component is 4*tolerance, for the RMS step it is 8/3*tolerance, for the RMS gradient it is 2/3*tolerance. All criteria are applied in the coordinate system specified by coordinates.

tolerance_e: (Default: tolerance/450) Convergence criterion for the energy change (in Hartree).

update_method

Update algorithm for the Hessian (e.g. in P-RFO TS search or quasi-Newton optimisation).

Allowed values:

'bofill': (default) Bofill update of the Hessian.
'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.

initial_hessian

Type of Hessian to calculate at the start or when resetting.

Allowed values:

'external': (default) External Hessian.
'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.

Nudged elastic band options

neb

Specifies if and how the nudged elastic band (NEB) method is to be used

Allowed values:

'no': (default) No NEB calculation will be performed. neb is automatically set to 'no' if no frag2 is specifed.
'free': The endpoints of the NEB path are minimised
'frozen': The endpoints of the NEB path are frozen
'perpendicular': The endpoints of the NEB path are minimised perpendicular to the NEB path, but their distance to the nearest other point is kept fixed.

nimages: (Default: 1) Specifies the number of images in a NEB calculation.

nebk: (Default: 0.01) Force constant for NEB calculations.

Dimer method options

dimer

Specifies if the dimer method for transition state (TS) search is to be used.

Allowed values:

False: (default) Do not use the dimer method
True: Use the dimer method for TS search

delta: (Default: 0.01) Dimer separation (in internal coordinates)

tsrelative: (Default: False) If True, the coordinates given in frag2 are given as displacements relative to the coordinates of frag (i.e. the frag2 coords will equal frag.coords + frag2.coords).

Microiterative options

Microiterative optimisation methods improve the efficiency of QM/MM optimisation by separating the active atoms into an inner region (which should contain the QM region) and an outer environment. After each step taken in the inner region, the environment is relaxed completely. The idea is to minimise expensive inner region (macroiterative) evaluations at the cost of increasing the number of environment (microiterative) cycles.

Microiterative methods only save overall calculation time if QM region calculations are not performed during the microiterations. If micro_esp_fit is set to true, the electrostatic influence of the QM region will be approximated by fitting point charges at the QM atom sites to an electrostatic potential generated by the QM code.

Microiterative optimisation has been implemented for minimisation with L-BFGS, transition state optimisation with P-RFO and the dimer method, and reaction path optimisation with NEB. In all cases the outer environment region is relaxed using L-BFGS. For the transition state and reaction path methods, this is equivalent to specifying spectator degrees of freedom (setting weights to zero) in a standard optimisation. This is useful for eliminating complications that can be caused by irrelevant degrees of freedom. For P-RFO this also means that the Hessian is only calculated over the inner degrees of freedom, which can dramatically reduce the cost of the Hessian calculation.

microiterative: (Default False) Specifies if microiterative optimisation is to be used.

inner_atoms: List of atoms in the inner region (which should contain all QM atoms). These atoms will have their positions relaxed in each macroiteration.

maxmicrocycle: (Default: 100) Maximum number of microiterative cycles permitted before switching back to macroiterations.

micro_esp_fit: (Default True) For microiterations, approximate the QM region using charges fitted to the electrostatic potential.

Warning

ESP charge fitting is only supported for MNDO, DFTB+ and GAMESS-UK QM theories.