Constrained geometry optimizations with common quantum chemistry codes
In molecular modeling, it is sometimes useful to perform a constrained optimization. The simplest case of a constrained optimization involves freezing one or more coordinates (x, y, z) of a chosen atom (or chosen set of atoms). In many cases, however, it is more chemically reasonable to freeze the value of a chosen distance, angle, or dihedral between a set of atoms. Here I describe how to apply these constraints in a few standard quantum chemistry codes. Recall that a distance is defined by two atoms, an angle by three atoms, and a dihedral by four atoms. These quantities are shown below in our test molecule, trans-1,2-dibromoethene.
distance defined by two atoms (Br-C)
angle defined by three atoms (Br-C-C)
torsion defined by four atoms (Br-C-C-Br)1. MOPAC
Begin by obtaining cartesian coordinates of your chemical structure and preparing an input file (.mop) as usual. For trans-1,2-dibromoethene, our coordinates are listed in the .xyz file available here. Recall that in a MOPAC input file, a flag of “1” indicates that the coordinate is allowed to relax while a flag of “0” indicates that it must remain frozen. To apply a constraint, we modify one line in our input file.
constraining Br—C bond length
C 0.17471 1 0.32948 1 0.01531 1
C 0.84673 1 1.14629 1 0.83212 1
Br 0.98631 1 -0.85129 1 -1.16613 1
Br 1.86 0 123.9 1 180 1 2 1 3 H -0.90892 1 10.29585 1 -0.01834 1
H 1.93034 1 1.17992 1 0.86579 1
constraining Br—C—C angle
C 0.17471 1 0.32948 1 0.01531 1
C 0.84673 1 1.14629 1 0.83212 1
Br 0.98631 1 -0.85129 1 -1.16613 1
Br 1.86 1 123.9 0 180 1 2 1 3 H -0.90892 1 10.29585 1 -0.01834 1
H 1.93034 1 1.17992 1 0.86579 1
constraining Br—C—C—Br torsion
C 0.17471 1 0.32948 1 0.01531 1
C 0.84673 1 1.14629 1 0.83212 1
Br 0.98631 1 -0.85129 1 -1.16613 1
Br 1.86 1 123.9 0 180 0 2 1 3 H -0.90892 1 10.29585 1 -0.01834 1
H 1.93034 1 1.17992 1 0.86579 1
On line four of our coordinates above, we specify that the distance between Br#4 and C#2 is 1.86 Å, the angle formed by Br#4—C#2—C#1 is 123.9°, and the dihedral formed by Br#4—C#2—C#1—Br#3 is 180°. A flag of “0” is used to apply a constraint: distance (left), angle (center), and torsion (right).
2. Gaussian
Prepare an input file (.com) with cartesian coordinates as usual and specify the keyword opt=ModRedudant on line two.
constraining Br—C bond length
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
4 2 = 1.86 B
4 2 Fconstraining Br—C—C angle
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
4 2 1 = 123.9 B
4 2 1 Fconstraining Br—C—C—Br torsion
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
4 2 1 3 = 180 B
4 2 1 3 FAfter the last line of the coordinates, we define the relevant constraints with the options “B” and “F” as shown above. “B” builds the coordinate, and “F” freezes it at the specified value.
3. GAMESS
Prepare an input file (.inp) with cartesian coordinates as usual. To apply a constraint, we must introduce the $ZMAT group in our input file.
First, we specify values for the keyword IZMAT. The first number defines the type of constraint IZMAT = 1 (distance), IZMAT = 2 (angle), IZMAT = 3 (torsion). The numbers following the first number (1, 2, or 3) list the atoms involved in the constraint. Two atoms should be listed for IZMAT = 1, three for IZMAT = 2, and four for IZMAT = 3. Next, we specify values for the keyword IFZMAT. The values of this keyword should be identical to those chosen for IZMAT. Finally, we set a value for the keyword FVALUE which specifies the value of our constrained distance, angle, or torsion.
constraining Br—C bond length
$ZMAT IZMAT(1)=1,4,2 IFZMAT(1)=1,4,2 FVALUE(1)=1.86 DLC=.TRUE. AUTO=.TRUE. $ENDC 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
constraining Br—C—C angle
$ZMAT IZMAT(1)=2,4,2,1 IFZMAT(1)=2,4,2,1 FVALUE(1)=123.9 DLC=.TRUE. AUTO=.TRUE. $ENDC 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
constraining Br—C—C—Br torsion
$ZMAT IZMAT(1)=3,4,2,1,3 IFZMAT(1)=3,4,2,1,3 FVALUE(1)=180 DLC=.TRUE. AUTO=.TRUE. $ENDC 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
Note that “1” in parentheses above because we only have one constraint in each case. Additional constraints can be added by incrementing these values, e.g. IZMAT(2), IFZMAT(2), FVALUE(2), and so on.
4. Quantum ESPRESSO
Prepare input file (.in) for geometry optimization with ion_dynamics = ‘damp,’ next, we define our constraints in the CONSTRAINTS group directly following K_POINTS. The first line of our CONSTRAINTS group defines the number of constraints and the tolerance. Each example below has one constraint with a tolerance of 1e-5. The following line in this group defines the type of constraint ‘distance’ (distance), ‘planar_angle’ (angle), or ‘torsion_angle’ (torsion) followed by the indices of the involved atoms and the value of the constraint.
constraining Br—C bond length
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
K_POINTS automatic
6 6 6 0 0 0
CONSTRAINTS
1 0.00001
'distance' 4 2 1.86
constraining Br—C—C angle
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
K_POINTS automatic
6 6 6 0 0 0
CONSTRAINTS
1 0.00001
'planar_angle' 4 2 1 123.9
constraining Br—C—C—Br torsion
C 0.17472 0.32948 0.01531
C 0.84673 1.14629 0.83212
Br 0.98631 -0.85129 -1.16613
Br 0.03517 2.32776 2.01288
H -0.90892 0.29586 -0.01834
H 1.93034 1.17992 0.86579
KPOINTS automatic
6 6 6 0 0 0
CONSTRAINTS
1 0.00001
'torsion_angle' 4 2 1 3 180
