The HDLC Optimiser

HDLCOpt is a deprecated module and will be removed in future versions of ChemShell.

Please use the DL-FIND optimiser instead.

The command line arguments for hdlcopt are kept here for reference.

The hdlcopt command

Command line arguments

Argument	Argument type	Mandatory	Default	To specify
theory=	string	no	mopac	Module used to compute energy and forces
coords=	Fragment tag	yes	none	Initial configuration of the system
result=	Fragment tag	no	hdlcopt.result	Optimised coordinates
residues=	Tcl List	no	none	Definition of residues
memory=	integer	no	see note 1	Number of BFGS steps to store for L-BFGS
update_procedure=	Tcl proc	no	none	Procedure to execute after each geometry update, typically used for graphics, see examples section
maxstep=	real	no	1.0	Maximum step size
contyp=	integer	no	0	Connectivity type, 0: primitive internals, 1: total connection
constraints=	Tcl List	no	none	Constraints specification
toler=	real	no	0.00045	Maximum gradient component convergence criterion
icconv=	integer	no	0	Convergence criteria (0: HDLC/Cartesian gradient and displacement, 1: Cartesian gradient)
list_option=	Output keyword	no	medium	how much output to generate
printl=	integer	no	0	Print level (alternative to list_option)
reghdl=	integer	no	0	Interval in optimisation steps to restart the optimiser, regenerate HDLC, and update the pair list.
cfact=	real	no	see note 2	Scaling factor for the Cartesians when constructing HDLC
maxfun=	integer	no	1000000	Maximum number of function and gradient evaluations
active_atoms=	Tcl List	no	all atoms	Specification of the atoms to be optimised (see note 3)
frozen=	Tcl List	no	none	Specification of the atoms to be frozen. Complementary to active_atoms (see note 3)
pdb_file=	filename	no	none	If specified, a pdb file is written with the residue specification as used in the coordinate conversion. Useful for setting up the residue list.
mstart=	integer	no	0	startup procedure: 0 start up from scratch, no checkpoint used at all 1 start up using coordinates from checkpoint file 2 start up using coordinates, resn and ctrl from checkpoint file 3 continue using the complete optimiser status from checkpoint
dump=	integer	no	0	dump the memory (write checkpoint file of the optimiser ) every dump= cycles (0: never)
TS-search specific options
core=	string	no	none	Specify the residue name of the reaction core, in which a TS search is performed. Specification of core or core_atoms will lead to a TS search.
core_atoms=	Tcl List	no	none	Specify the atoms forming the reaction core, in which a TS search is performed. Specification of core or core_atoms will lead to a TS search.
nvar=	integer	no	0	Number of degrees of freedom to be treated by P-RFO. Obsolete: use core or core_atoms now!
mode=	integer	no	1	Hessian mode to search in P-RFO
ctfirst=	integer	no	0	Connectivity type of the first HDLC residue (0: primitive internals, 1: total connection, 2: Cartesians, 3: pure DLC (primitive internals only) 4: total connection scheme (internals only)). This is of use for TS search, as the core is always the first residue.
recalc=	integer	no	0	Interval for recomputing the finite-difference Hessian (Number of Hessian updates before it is recalculated). Use recalc=-1 for no hessian calculation, recalc=0 for calculating it only at the start. Actually, reghdl is the better option for this task.
updmth=	integer	no	see note 5	Hessian update method
lockon=	integer	no	0	Lock the Hessian eigenmode (0: mode switching, 1: mode lock)
delta=	real	no	0.01	Displacement for calculating the finite-difference Hessian. Atomic units (Bohr) if Cartesians or the total connection scheme are used (ctfirst= 1 or 2), otherwise undefined (but close to atomic) units
Options specific to micro-iterative QM/MM search
inner_atoms=	Tcl List	no	none	List of atoms in the inner part (should contain all QM atoms). Specification of inner_atoms or inner_residues will lead to a microiterative QM/MM optimisation.
inner_residues=	Tcl List	no	none	List of residues in the inner part (should contain all QM atoms). Specification of inner_atoms or inner_residues will lead to a microiterative QM/MM optimisation.
include_res=	Boolean	no	false	Only used when inner_atoms is specified: include the whole residue into the inner region if at least one of its atom is included in inner_atoms. This changes the number of atoms in the inner region.
mcore=	Boolean	no	-	If true, microiterations (outer region) are done by representing the QM part by esp charges. This only makes sense with hybrid QM/MM calculations. At the moment it only works with qm_therory being mndo, gamess, or turbomole. In these cases, the default is true, otherwise false.
toler_ratio=	real	no	10	The ratio between the tolerance (convergence criterion) of the macro-iterations and the micro-iterations. Numbers > 1 mean tighter convergence for the micro-iterations (recommended).

Note 1: The lower value of either 3*natoms and 40 is taken as default

Note 2: Default: 1/(number of DF in the residue), cfact>0: cfact is used, cfact<0: -cfact/(number of DF in the residue) is used

Note 3: If residues are specified, they are intersected with active_atoms, and the complement of active_atoms is added as residue "fixed". If no residues are specified, a residues list containing the complement of active_atoms as residue "fixed" is created

Note 5: The default depends on the parameter mode: if mode=0, updmth=2, else updmth=1

Update_procedure (example)

proc hdlcopt_update { args } {
    parsearg update {coords} $args
    write_xyz coords=$coords file=update.xyz
    exec cat update.xyz >> optimise.xyz
    end_module
}

If you want the coordinate and result files to be updated each step, include the statement flush_all_objects in the update_procedure.

Practical hints

• Storage problems:
  • Lower the number of remembered steps (memory=...)
  • Move atoms from the reaction core to the environment (nvar=... and residues=...)
  • Use Cartesian coordinates for the bulk of a large system, especially for the atoms far from the QM region in a QM/MM calculation (residues=...)
• Convergence problems:
  • HDLCopt reacts very sensitively to an inconsistent energy / gradient pair. The result is no or very slow convergence
  • Make sure you provided the option keywords=precise to MOPAC
  • MNDO provides generally a very accurate gradient
  • Force fields: make sure that the MM code never updates the nonbonded pair list on its own (see below)
• Definition of residues:
  • Use the argument residues=... and the commands for the manipulation of residue membership lists
  • Recommended size of a residue is 3 to 50 atoms
  • The atoms of the residue named "fixed" are kept fixed
  • The argument active_atoms=... can be used if most of the atoms are to be kept fixed
  • If no core= or core_atoms= are given, the first residue represents the reaction core in a TS search
  • Try to cut single bonds only
  • Try to avoid cutting rings and other crosslinked subsystems
  • Do not put more than one molecule into a residue
  • Use the PDB fragmentation for biomolecules
• Number of remembered L-BFGS steps (memory=...):
  • Never exceed the number of steps between two restarts (reghdl=..., default: never restart)
  • For force fields, a value between 5 and 20 is sufficient
  • For a potential from electronic structure calculations, values up to the number of degrees of freedom make sense, values up to 200 are recommended
• Reaction core (nvar=...):
  • Including the atoms up to the second or third neighbours of the transferred atoms (approx. 6 to 20 atoms) into the reaction core is generally sufficient
  • Use the total connection scheme for the reaction core (ctfirst=1) and primitive internal coordinates for the environment (contyp=0)
  • Don't forget: the first residue defined in the residue membership list (residues=...) is the reaction core. nvar=... must be three times the number of atoms in the first residue
• HDLC breakdowns (the G matrix is only diagonalised at start and restarts. Otherwise, no updated Hessian could be used):
  • The default value of the scaling factor of the Cartesian coordinates when delocalising cfact=... is inappropriate for large systems, since the external coordinates play an important role for a fragment surrounded by other fragments. Try using values between 0.25 and 1. Both efficiency and likeliness of coordinate breakdowns increase with decreasing values of cfact=...
  • If the HDLC of the same residue persistently break down, the optimiser automatically removes this residue from the list. For subsequent optimisations, e.g. a transition state search, it can be removed manually from the list (residues=...) using the command inlist)
  • For large systems, consider doing the first 200 to 1000 optimisation steps (maxfun=200) in Cartesian coordinates (no residues=... argument) or try using frequent restarts (reghdl=10) if starting very far from a minimum
• Pair list update (reghdl=...):
  • At each restart of hdlcopt, the pair list for nonbonded interactions is updated automatically
  • Force regular pair list updates using (reghdl=...). Recommended intervals are between every 20 and every 200 steps
  • Under no circumstance may the MM code update its pairlist on its own. This is especially important when using an external force field code such as CHARMM (set inbfrq to zero). The GROMOS interface automatically prevents pair list updates unless forced by reghdl=...

Related reading

S. R. Billeter, A. J. Turner and W. Thiel, Linear scaling geometry optimization and transition state search in hybrid delocalised internal coordinates, Phys. Chem. Chem. Phys., 2 (2000), 2177-2186.

• HDLCopt:
  • S. R. Billeter, A. J. Turner and W. Thiel, Linear scaling geometry optimization and transition state search in hybrid delocalised internal coordinates, Phys. Chem. Chem. Phys., 2 (2000), 2177-2186.
• L-BFGS:
  • J. Nocedal, Updating quasi-Newton matrices with limited storage, Mathematics of Computation, 35 (1980), 773-782.
  • D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Math. Prog., 45 (1989), 503-528.
• P-RFO:
  • A. Banerjee, N. Adams, J. Simons and Ron Shepard, Search for Stationary Points on Surfaces, J. Phys. Chem., 89 (1985), 52-57.
  • J. Baker, An Algorithm for the Location of Transition States, J. Comput. Chem. 7 (1986), 385-395.
• Microiterative optimisation:
  • A. J. Turner, V. Moliner and I. H. Williams, Transition-state structural refinement with GRACE and CHARMM: Flexible QM/MM modelling for lactate dehydrogenase, Phys. Chem. Chem. Phys. 1 (1999), 1323-1331.
  • S. R. Billeter, A. J. Turner and W. Thiel, Linear scaling geometry optimization and transition state search in hybrid delocalised internal coordinates, Phys. Chem. Chem. Phys., 2 (2000), 2177-2186.
• Internal coordinates:
  • H. B. Schlegel, Adv. Chem. Phys. 67 (1987), 249.
  • P. Pulay and G. Fogarasi, Geometry optimization in redundant internal coordinates, J. Chem. Phys. 96 (1992), 2856-2867.
  • J. Baker, Techniques for Geometry Optimization: A Comparison of Cartesian and Natural Internal Coordinates, J. Comp. Chem. 14 (1993), 1085-1100.
  • C. Peng, P. Y. Ayala, H. B. Schlegel and M. J. Frisch, Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States, J. Comput. Chem. 17 (1996), 49-56.
  • J. Baker, A. Kessi and B. Delley, The generation and use of delocalized internal coordinates in geometry optimization, J. Chem. Phys. 105 (1996), 192-212.
  • O. Farkas and H. B. Schlegel, Methods for geometry optimization of large molecules. I. An O(N-2) algorithm for solving systems of linear equations for the transformation of coordinates and forces, J. Chem. Phys. 109 (1998), 7100-7104.
  • M. v. Arnim and R. Ahlrichs, Geometry optimization in generalized natural internal coordinates, J. Chem. Phys. 111 (1999), 9183-9190.
  • J. Baker, D. Kinghorn and P. Pulay, Geometry optimization in delocalized internal coordinates: An efficient quadratically scaling algorithm for large molecules, J. Chem. Phys. 110 (1999), 4986-4991.