Numerical Hessian for a plannar molecule with n atoms

n atoms:(x1, y1, 0), (x2, y2, 0), (x3, y3, 0), ...

A total of (1 + 9*n^2) energy calculations need to be performed:

1: at equilibrium geometry.
n*1: for 1z-1z 2z-2z,3z-3z,....
[n*(n-1)/2]*2: for 1z-2z, 1z-3z, 1z-4z, 2z-3z, 2z-4z, 3z-4z.
[(2n)^2]*4/2: for the second derivatives with respect to x(s) and y(s)

1 + n + n*(n-1) + 8*n^2 = 1 + 9*n^2

For example, Numerical Hessian for a plannar molecule with 3 atoms

3 atoms:(x1, y1, 0),(x2, y2, 0),(x3, y3, 0)

82 energy calculations need to be performed:

1: at equilibrium geometry.
3: for 1z-1z 2z-2z,3z-3z
3*2: for 1z-2z, 1z-3z, 2z-3z.
2*6^2: for the second derivatives with respect to x(s) and y(s)

1 + 3 + 6 + 72 = 82 (verified by an AM1 method=fullnum Hessian calculation on HOCl)

For a planar molecule such as trans-HNOH with 4 atoms, 1+9*4^2=145 energy calculations need to be carried out to obtain the numerical Hessian.

 $force temp(1) = 1100.00, 1200.00, 1300.00, 1400.00, 1500.00,
                  1600.00, 1700.00, 1800.00, 1900.00, 2000.00
  method=fullnum projct=.t.
 $end
 $contrl maxit=200 mult=2 scftyp=uhf runtyp=hessian coord=unique $end
 $basis gbasis=am1 $end
 $system timlim=1000000 mwords=80 $end
 $data
 HNOH
 Cs

 O           8.0   0.0094634591  -0.0934516366   0.0000000000
 N           7.0   1.2729078083   0.1027207551   0.0000000000
 H           1.0  -0.4112060327   0.7952578868   0.0000000000
 H           1.0   1.7288347653  -0.8045270053   0.0000000000
 $end

For a non-planar  molecule with 4 atoms, see the post on H2NO.

Numeric Hessian for a non-planar H2NO molecule with a Cs Point Group

Both H atoms are bonded to N; the point group is Cs; the symmetry plane contains N-O bond.

$force temp(1) = 10.00, 20.00, 30.00, 40.00, 50.00,
60.00, 70.00, 80.00, 90.00, 100.00
method=fullnum projct=.t.
$end
$contrl maxit=200 mult=2 scftyp=rohf runtyp=hessian coord=unique $end
$contrl ispher=+1 $end
$contrl cctyp=cr-ccl $end
$basis gbasis=cct $end
$system timlim=1000000 mwords=400 memddi=0 $end
$scf soscf=.t. dirscf=.t. $end
$statpt nstep=100 opttol=0.00001 $end
$data
H2NO
Cs

N 7.0 -0.1357158331 0.0890286127 0.0000000000
H 1.0 -0.3822861363 -0.3680309994 0.8687014208
O 8.0 0.0165139267 1.3587083668 0.0000000000
$end

The following 2nd derivatives are evaluated:

Between atoms 1 and 1 9 (2x-x 4x-y 2y-y 1z-z)
Between atoms 1 and 2 (not calculated because 1-2 = 1-3)
Between atoms 1 and 3 (4*3*3=36)
Between atoms 2 and 2 (not calculated because 2-2 = 3-3)
Between atoms 2 and 3 21
Between atoms 3 and 3 18
Between atoms 4 and 1 18
Between atoms 4 and 2 (not calculated because 4-2 = 4-3)
Between atoms 4 and 3 (4*3*3=36)
Between atoms 4 and 4 9

Total = 147 + 1 = 148

Note that due to the Cs symmetry, E(z+dz) = E(z-dz) at z=0


E"(z=0) with respect to z = [E(z+dz) + E(z-dz) - 2E(z)]/dz^2


E"(z=0) = 2[E(z=dz) - E(z=0)]/dz^2


Also, at z=0, E(x+dx,z+dz) = E(x+dx,z-dz);
E(x-dx,z-dz) = E(x-dx,z+dz) due to the Cs symmetry.

Therefore,
E"(z=0) with respect to x and z
= [E(x+dx,z+dz) + E(x-dx,z-dz) - E(x-dx, z+dz) - E(x+dx, z-dz)]/4dxdz = 0

Numerical Hessian

The total number of energy calculations for a N-Cartesian-coordinate system is 2N^2+1.

If N=3n assuming each particle has 3 Cartesian coordinates, then the number of energy calculations = 2*(3n)^2+1 = 18*n^2+1


The single point energy of the geometry is evaluated first.


Then the coordinates are disturbed in the following ways:


For each diagonal element, disturb only 1 coordinate (x) at a time. Evaluate E(x+dx) and E(x-dx). The 2nd derivative can be derived in the following ways:


E'(x+dx/2)=[E(x+dx)-E(x)]/dx


E'(x-dx/2)=[E(x)-E(x-dx)]/dx


E" = [E'(x+dx/2)-E'(x-dx/2)]/dx


E" = [E(x+dx) + E(x-dx) - 2E(x)]/dx^2


For each off-diagonal elements, disturb 2 coordinates (x and y) at the same time. Evaluate E(x-dx,y-dy), E(x-dx,y+dy), E(x+dx,y-dy), and E(x+dx,y+dy).


The 2nd derivative can be derived in the following way:


E'(x+dx,y) with respect to y = [E(x+dx,y+dy)-E(x+dx,E(y-dy))]/2dy


E'(x-dx,y) with respect to y = [E(x-dx,y+dy)-E(x-dx,E(y-dy))]/2dy


E"(x,y) = {[E(x+dx,y+dy)-E(x+dx,E(y-dy))]/2dy - [E(x-dx,y+dy)-E(x-dx,E(y-dy))]/2dy}/2dx


E"(x,y) = {[E(x+dx,y+dy)-E(x+dx,E(y-dy))] - [E(x-dx,y+dy)-E(x-dx,E(y-dy))]}/4dxdy


E" = [E(x+dx,y+dy) + E(x-dx,y-dy) - E(x-dx, y+dy) - E(x+dx, y-dy)]/4dxdy


Note that E(x,y) is not used in the above off-diagonal element calculation.


In summary, 2N+1 energy calculations need to be done to obtain the N diagonal elements of the energy second derivative matrix (Hessian); 4[N(N-1)/2] energy calculations need to be done to obtain the N(N-1)/2 off-diagonal elements of the Hessian matrix.


2N+1+4[N(N-1)/2]=2N^2+1


For example, a 4-particle system has 12 Cartesian coordinates. 2(12^2)+1=289 energy calculations need to be carried out to obtain the Hessian matrix. A 5-particle system has 15 Cartesian coordinates. 2(15^2)+1=451 energy calculations need to be carried out to obtain the Hessian matrix.


The use of the internal coordinates reduces the number of coordinates by 5 or 6. A linear system has 3n-5 internal coordinates while a non-linear system has 3n-6 internal coordinates, where n is the number of particles.


The number of coordinates can be further reduced if the studied system has a symmetry element.