Applying Harmonic Constraints in Jaguar Calculations

Sometimes you don’t want to entirely freeze the value of a coordinate, but allow it to vary within defined limits. You can then optimize the geometry while allowing for small variations of coordinates that you know should remain essentially the same. You can achieve this end by applying harmonic constraints. Harmonic constraints are additional potential energy terms in the form of a harmonic potential with a given force constant k, centered at an atom or on a particular internal coordinate value. The constraint can include a region where the potential is zero inside the constrained area. The form of the potential is as follows:

V = (k/2) |da|2 |d| ≥ a
V = 0 |d| < a

For atomic (Cartesian) harmonic constraints, the value of d is the distance from the atom. For internal coordinate harmonic constraints (bond length, bond angle, or dihedral angle), the value of d is the difference between the coordinate value and its desired value: d = rc, where r is the coordinate and c is the center or target value. The constraint is then more like a dynamic constraint, in which the coordinate is at the target value at the end of the optimization. These extra potential energy terms are listed in the output and used to determine geometry convergence, but they are not included in the final energy.

You can also use an exponent other than 2, to create a cubic or quartic constraint, for example: to do this, set the constraint_exponent keyword in the gen section. Higher exponents give continuous second derivatives at the edge of the flat bottom, and provide a harder wall. The general form of the potential is

V = (k/n!) |da|n.

Harmonic constraints must be set in the coord section of the input file, and can be set on the Cartesian position of an atom or on any bond length, angle, or dihedral angle. To set a harmonic constraint, adding #hc or #HC after the coordinate, followed by the force constant. If you want to specify a region of zero potential, add the half-width of this region, a, after the force constant. If you want to specify a target value c for an internal coordinate, it must follow the radius a. A target value cannot be specified for a Cartesian harmonic constraint. The units of the force constant, the half-width, and the target coordinate values are specified by the gen section keywords iunit and eunit.

The Cartesian position is specified by a single atom label, as in the following example.

&coord
C1 #hc 10.0 
&

The following example specifies a harmonic constraint on a bond length, with a force constant of 10.0 kcal mol−1 Å−2, a width of 0.1 Å, and a target bond length of 1.5 Å:

&coord
C1 C2 #hc 10.0 0.1 1.5 
&

In this constraining potential, the bond length can vary freely between 1.4 Å and 1.6 Å, but the energy rises if it goes outside this region.

Assuming units of kcal/mol for energy and degrees for angles, the following ranges are recommended for harmonic constraints of nonbonded atoms:

  • bond stretch: 10 - 1000
  • angle bend or torsion: 0.0001 to 0.1

Values for bonded atoms can be somewhat larger, but much higher values can begin to freeze the coordinate, while values as large as 10 for a torsional constraint can actually introduce errors into the interconversion of the internal coordinate hessian and the Cartesian hessian, and therefore should be avoided.