Geometry optimization and vibrational analysis¶
In this example, we show how to optimize the minimum energy geometry and calculate the harmonic frequencies and normal modes of a polyatomic molecule.
Geometry optimization¶
We first need a potential energy surface (PES) to work with.
We will use one for H2O published by
Polyansky et al (J. Chem. Phys. 105, 6490, 1996), which is
available as a built-in DFun object.
>>> import nitrogen as n2
>>> from nitrogen.pes.library.h2o_pjt1996 import PES
>>> PES.nx
9
The PES object returns the potential energy surface
(in cm-1) as a
function of the 9 Cartesian coordinates (in Angstroms)
in H, O, H order.
We wish to perform a geometry optimization with respect to curvilinear internal coordinates. To do so, we must first define a coordinate system transforming between internal and Cartesian coordinates, and then construct a PES with internal coordinate arguments.
>>> zmat = """
... H
... O 1 r1
... H 2 r2 1 theta
... """
>>> cs = n2.coordsys.ZMAT(zmat)
>>> V = cs ** PES
This code snippet defines the internal Z-matrix coordinate system
in terms of the two OH bond lengths and HOH bond angle and then
defines V as a new DFun that accepts
internal coordinate arguments, transforms these to Cartesians, and then
evaluates the original PES function.
NITROGEN has various optimizer routines. A good general purpose
optimizer is nitrogen.pes.opt_bfgs(), which is a standard
Broyden–Fltetcher–Goldfarb–Shanno (BFGS) algorithm. Let’s call it
with an initial geometry of r(OH) = 1.0 Å and \(\theta\)(HOH)
= 100\(^\circ\) and extra information printed.
>>> qmin,Vmin = n2.pes.opt_bfgs(V, [1.0, 1.0, 100.0], disp = True)
Step Value |grad|
-----------------------------
1 7.0178e+02 2.0302e+04 ... [ 1. 1. 100.]
2 1.9115e+01 4.0473e+03 ... [ 0.95115476 0.95115476 104.70196102]
3 4.2499e-01 5.9103e+02 ... [ 0.95893928 0.95893928 104.46238834]
4 2.2943e-04 1.3751e+01 ... [ 0.95794415 0.95794415 104.49873555]
5 2.8546e-09 4.8817e-02 ... [ 0.95792042 0.95792042 104.49964887]
6 3.6116e-15 3.6321e-06 ... [ 0.9579205 0.9579205 104.49964697]
7 1.3489e-18 8.6905e-08 ... [ 0.9579205 0.9579205 104.499647 ]
Convergence reached, |g| = 8.690e-08
7 gradient(s) and 1 Hessian(s) were calculated.
>>> qmin # Minimum energy geometry (rOH, rOH, aHOH)
array([ 0.9579205, 0.9579205, 104.499647 ])
>>> Vmin # Minimum energy value (cm^-1)
1.3489154554107497e-18
nitrogen.pes.opt_bfgs() has various options for changing the convergence
tolerance, constraining one or more arguments to fixed values, supplying a pre-computed
Hessian, etc.
Curvilinear vibrational analysis¶
Now that we have optimized the equilibrium geometry, we can perform a harmonic vibrational normal-mode analysis. We will first do this explicity using our Z-matrix coordinate system using the standard curvilinear GF approach (see Wilson, Decius, and Cross 1955). Before doing so, we need to define the atomic masses.
>>> masses = n2.constants.mass(['H','O','H'])
>>> masses
[1.00782503224, 15.9949146196, 1.00782503224]
>>> omega, nctrans = n2.pes.curvVib(qmin, V, cs, masses)
nitrogen.pes.curvVib() returns harmonic frequencies and a
linear transformation object containing the normal-mode displacement vectors.
>>> omega # harmonic frequencies (* hc, in cm-1)
array([1649.58906249, 3830.38088976, 3940.96386738])
>>> nctrans.T # columns are the normal-coordinate displacement vectors
array([[ 7.11357961e-03, 6.74635431e-02, 6.76662154e-02],
[ 7.11357961e-03, 6.74635431e-02, -6.76662154e-02],
[-1.25106206e+01, 9.57536309e-02, -8.98152710e-16]])
The displacement vectors are scaled to equal reduced dimensionless normal coordinates, \(q\), i.e., the coordinates in which the harmonic potential is \(V = \frac{1}{2} \omega q^2\), where \(\omega\) is the harmonic frequency (in energy units).
A new coordinate system can be constructed using these curvilinear normal coordinates. Let’s build this and verifying that the equilibrium geometry is \(q = 0\) and that the Hessian in this coordinate system is diagonal with elements equal to the harmonic frequencies.
>>> cs2 = nctrans ** cs # The normal-mode coordinate system
>>> V2 = nctrans ** V # The PES w.r.t normal-mode coordinates (q)
>>> qmin,Vmin = n2.pes.opt_bfgs(V2, [0.1, 0.2, 0.3], disp = True)
Step Value |grad|
-----------------------------
1 2.4367e+02 1.2706e+03 ... [0.1 0.2 0.3]
2 4.3457e+00 1.8472e+02 ... [-0.00926409 -0.0339949 -0.03199741]
3 6.3725e-02 2.2157e+01 ... [0.00097352 0.00393639 0.00411375]
4 1.5009e-05 3.4152e-01 ... [9.69510200e-06 5.54551059e-05 6.77409462e-05]
5 2.1711e-10 1.2347e-03 ... [-1.98193817e-07 -3.00724345e-07 7.64542669e-08]
6 2.4936e-13 4.2101e-05 ... [-6.66707056e-09 -7.97745234e-09 6.80267296e-09]
7 2.0227e-18 1.2603e-07 ... [-4.31359815e-11 1.06617424e-11 3.09024816e-11]
Convergence reached, |g| = 1.260e-07
7 gradient(s) and 1 Hessian(s) were calculated.
>>> hes = V2.hes(qmin)[0]
>>> hes
array([[ 1.64958906e+03, 3.39518280e-09, -7.26838984e-09],
[ 3.39518280e-09, 3.83038089e+03, -5.39473322e-08],
[-7.26838984e-09, -5.39473322e-08, 3.94096387e+03]])
>>> np.allclose(np.diag(hes), omega) # diagonal elements equal omega?
True
Everything checks out.
Rectilinear vibrational analysis¶
The standard Watson Hamiltonian is based on rectilinear normal
coordinates. These can be calculated using
nitrogen.vpt.calc_rectilinear_modes(), which first requires
calculating the Hessian with respect to Cartesian displacements. To do that,
we first evaluate the equilibrium Cartesian position using our curvilinear
equilibrium geometry from above and rotate it to the principal axis system.
>>> Xe = cs2.Q2X(qmin)[0] # Cartesian equilibrium geometry
>>> Xe,_,_ = n2.angmom.X2PAS(Xe, masses) # Rotate to PAS
>>> hes = PES.hes(Xe)[0] # The Cartesian Hessian at Xe
>>> omega_rect, T = n2.vpt.calc_rectilinear_modes(hes, masses)
>>> omega_rect
array([5.00232981e-03, 1.36059436e-05, 7.38904032e-06, 2.88647058e-05,
2.93979778e-03, 6.45802187e-03, 1.64958906e+03, 3.83038089e+03,
3.94096387e+03])
The harmonic frequencies in omega_rect include the 3 translational and
3 rotational modes, which equal zero. The vibrational frequencies equal
the those calculated with the curvilinear GF method above.
The normal-mode Cartesian displacement vectors are returned as the columns of
T. By default, calc_rectilinear_modes() normalizes
the vibrational vectors to the same reduced dimensionless normal coordinates,
\(q\), as above.
To request the displacement vectors with respect to mass-weighted Cartesians,
use the norm keyword. These are normalized to unity modulus.
>>> omega_rect, L = n2.vpt.calc_rectilinear_modes(hes, masses, norm = 'mass-weighted')
>>> np.allclose(L.T @ L, np.eye(len(L))) # L is orthonormal
True
The mass-weighted displacement vectors can be used to calculate Coriolis coupling constants
>>> Lvib = L[:,6:] # The vibrational vectors
>>> zeta = n2.vpt.calc_coriolis_zetas(Lvib)
>>> zeta # (mode i, mode j, axis k)
array([[[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 7.21324074e-17, -2.51543557e-17, -2.54098520e-10],
[-7.78383835e-17, 1.35562504e-16, 9.99931994e-01]],
[[-7.21324074e-17, 2.51543557e-17, 2.54098520e-10],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[-9.08401587e-17, -1.21849559e-16, 1.16621858e-02]],
[[ 7.78383835e-17, -1.35562504e-16, -9.99931994e-01],
[ 9.08401587e-17, 1.21849559e-16, -1.16621858e-02],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]])
Harmonic centrifugal distortion constants can also be calculated
>>> B0,CD = n2.vpt.analyzeCD(Xe, omega_rect[6:], Lvib, masses, printing = True)
==========================================
Harmonic centrifugal distortion analysis
==========================================
cm^-1 MHz
-------------- --------------
Ae 2.73812E+01 820869.10098
Be 1.45785E+01 437052.06382
Ce 9.51334E+00 285202.71759
A' 2.73832E+01 820927.74239
B' 1.45804E+01 437110.70523
C' 9.51040E+00 285114.75547
A' - Ae 1.95607E-03 58.64141
B' - Be 1.95607E-03 58.64141
C' - Ce -2.93410E-03 -87.96212
------------- (Ir) -------------
A(A) 2.73816E+01 820878.74053
B(A) 1.45862E+01 437283.44557
C(A) 9.50628E+00 284991.01699
A(A)-Ae 3.21541E-04 9.63955
B(A)-Be 7.71806E-03 231.38175
C(A)-Ce -7.06157E-03 -211.70060
A(S) 2.73817E+01 820883.30410
B(S) 1.45854E+01 437259.53118
C(S) 9.50695E+00 285011.28053
A(S)-Ae 4.73765E-04 14.20312
B(S)-Be 6.92037E-03 207.46736
C(S)-Ce -6.38565E-03 -191.43706
------------ (IIIr) ------------
A(A) 2.73717E+01 820582.11772
B(A) 1.46211E+01 438329.09419
C(A) 9.48129E+00 284241.99118
A(A)-Ae -9.57273E-03 -286.98326
B(A)-Be 4.25971E-02 1277.03036
C(A)-Ce -3.20464E-02 -960.72640
A(S) 2.73974E+01 821351.94126
B(S) 1.45755E+01 436963.44257
C(S) 9.50614E+00 284986.77628
A(S)-Ae 1.61058E-02 482.84029
B(S)-Be -2.95609E-03 -88.62125
C(S)-Ce -7.20303E-03 -215.94131
sigma ............ 6.050359
---------------------------------------
Kivelson-Wilson parameters
---------------------------------------
DJ 9.59431E-04 28.76303
DK 2.47000E-02 740.48861
DJK -3.75105E-03 -112.45370
R5 1.05198E-03 31.53761
R6 -1.02158E-04 -3.06262
---------------------------------------
A-reduced (Ir) parameters
---------------------------------------
DeltaJ 1.16375E-03 34.88826
DeltaK 2.57216E-02 771.11477
DeltaJK -4.97695E-03 -149.20510
deltaJ 4.64437E-04 13.92346
deltaK 3.68404E-04 11.04448
---------------------------------------
S-reduced (Ir) parameters
---------------------------------------
DJ 1.13330E-03 33.97555
DK 2.55694E-02 766.55120
DJK -4.79428E-03 -143.72882
d1 -4.64437E-04 -13.92346
d2 -1.52224E-05 -0.45636
---------------------------------------
A-reduced (IIIr) parameters
---------------------------------------
DeltaJ 1.20005E-02 359.76656
DeltaK 2.57216E-02 771.11477
DeltaJK -3.74873E-02 -1123.84000
deltaJ 4.95395E-03 148.51569
deltaK -1.78079E-02 -533.86879
---------------------------------------
S-reduced (IIIr) parameters
---------------------------------------
DJ 7.03185E-03 210.80954
DK 8.78263E-04 26.32968
DJK -7.67524E-03 -230.09788
d1 -4.95395E-03 -148.51569
d2 -2.48434E-03 -74.47851
==========================================