5.4 Force Field Calculations
5.4.1 Purpose
uff preoptimizes a structure and calculates an analytical Hessian which can be used as a start Hessian in a geometry optimization. This will accelerate the convergence of an optimizations. For optimizations in Cartesian space this will be faster by a factor of two for any molecule.
5.4.2 How to Perform a Uff Calculation
You have to generate Cartesian coordinates (file coord), nothing else. You can start an single-point calculation calculation by typing
uff
To start an uff geometry optimization, one has to change the number of cycles (parameter maxcycle) in the block $uff in the file control. The ouput is the optimized structure (file coord), the analytical gradient (file uffgradient) and the analytical cartesian hessian (file uffhessian0-0). Furthermore the control file will be modified:
$forceinit on
carthess
$uffhessian file=uffhesian0-0These commands have the effect to initialize the force constant matrix for a geometry optimization with the hessian one.
In some cases uff cannot recognize the connectivity, then one can specify the connectivity in the file ufftopology. The program will calculate the bond, angle, torsion, inversion and non-bonded terms (force field terms) based on the connectivity specified in the topology file.
5.4.3 The Uff implementation
The uff implementation follows the paper by Rappé [11]. The energy expression in uff is as follows:
\[\begin{align}
\mathrm{E}\,{}_{UFF} =& \sum^{N_B} \frac{1}{2}\cdot K_{IJ} \cdot
\left( r-r_{IJ} \right)^2 \\
+& \sum^{N_A} \left\{
\begin{array}{rl} \frac{K_{IJK}}{4}
\left( 1 - \cos (2\theta) \right) & \text{linear case}\\
\frac{K_{IJK}}{9} \left( 1 - \cos (3 \theta) \right) &
\text{trigonal planar case}\\
\frac{K_{IJK}}{16} \left( 1 - \cos (4 \theta) \right) &
\text{quadratic planar case}\\
\frac{K_{IJK}}{16} \left( 1 - \cos (4 \theta) \right) &
\text{octahedral case}\\
K_{IJK} \cdot \left( C^A_0 + C^A_1 \cos \theta + C^A_2
\cos(2\theta) \right) &
\text{general case}\\
\end{array} \right. \\
+& \sum^{N_T}\frac{1}{2}\cdot V_{\phi} \cdot
\left(1-\cos \left(n\phi_0\right) \cos(n\phi) \right) \\
+& \sum^{N_I} V_\omega \cdot
\left(C^I_0 + C^I_1\cos \omega + C^I_2 \cos 2 \omega \right) \\
+& \sum^{N_{nb}} D_{IJ} \cdot \left( -2 \left(\frac{x_{IJ}}{x}\right)^6
+ \left(\frac{x_{IJ}}{x}\right)^{12} \right) \\
+& \sum^{N_{nb}} \frac{q_I \cdot q_J}{\epsilon \cdot x}
\end{align}\](5.1)
The Fourier coefficients \(C^A_0,C^A_1,C^A_2\) of the general angle terms are evaluated as a function of the natural angle \(\theta_0\): \[\begin{align} C^A_2 &= \frac{1}{4 \sin^2 {\theta_0}} \\ C^A_1 &= -4 \cdot C^A_2 \cos{\theta_0} \\ C^A_0 &= C^A_2 \left( 2 \cos^2{\theta_0}+1 \right) \end{align}\](5.2–5.4) The expressions in the energy term are:
- \(N_B, N_A, N_T, N_I, N_{nb}\)
-
the numbers of the bond-, angle-, torsion-, inversion- and the non bonded-terms.
- \(K_{IJ}, K_{IJK}\)
-
force constants of the bond- and angle-terms.
- \(r,r_{IJ}\)
-
bond distance and natural bond distance of the two atoms \(I\) and \(J\).
- \(\theta,\theta_0\)
-
angle and natural angle for three atoms \(I-J-K\).
- \(C^A_0,C^A_1,C^A_2\)
-
Fourier coefficients of the general angle terms.
- \(\phi,\phi_0\)
-
torsion angle and natural torsion angle of the atoms \(I-J-K-L\).
- \(V_\phi\)
-
height of the torsion barrier.
- \(n\)
-
periodicity of the torsion potential.
- \(\omega\)
-
inversion- or out-of-plane-angle at atom \(I\).
- \(V_\omega\)
-
height of the inversion barrier.
- \(C^I_0,C^I_1,C^I_2\)
-
Fourier coefficients of the inversions terms.
- \(x,x_{IJ}\)
-
distance and natural distance of two non bonded atoms \(I\) and \(J\).
- \(D_{IJ}\)
-
depth of the Lennard–Jones potential.
- \(q_I,\epsilon\)
-
partial charge of atoms \(I\) and dielectric constant.
One major difference in this implementation concerns the atom types. The atom types in Rappé’s paper have an underscore "_". In the present implementation an sp\(^3\) C atom has the name "C 3" instead of "C_3". Particularly the bond terms are described with the harmonic potential and the non-bonded van der Waals terms with the Lennard–Jones potential. The partial charges needed for electrostatic non-bonded terms are calculated with the Charge Equilibration Modell (QEq) from Rappé [65]. There is no cutoff for the non-bonded terms.
The relaxation procedure distinguishes between molecules wih more than 90 atoms and molecules with less atoms. For small molecules it consists of a Newton step followed by a linesearch step. For big molecules a quasi-Newton relaxation is done. The BFGS update of the force-constant matric is done [66, 67, 60, 68]. Pulay’s DIIS procedure is implemented for big molecule to accelerate the optimization [69, 59].
The coordinates for any single atom can be fixed by placing an ’f’ in the third to eighth column of the chemical symbol/flag group. As an example, the following coordinates specify acetone with a fixed carbonyl group:
$coord
2.02693271108611 2.03672551266230 0.00000000000000 c
1.08247228252865 -0.68857387733323 0.00000000000000 c f
2.53154870318830 -2.48171472134488 0.00000000000000 o f
-1.78063790034738 -1.04586399389434 0.00000000000000 c
-2.64348282517094 -0.13141435997713 1.68855816889786 h
-2.23779643042546 -3.09026673535431 0.00000000000000 h
-2.64348282517094 -0.13141435997713 -1.68855816889786 h
1.31008893646566 3.07002878668872 1.68840815751978 h
1.31008893646566 3.07002878668872 -1.68840815751978 h
4.12184425921830 2.06288409251899 0.00000000000000 h
$end