nitrogen.pes

Potential energy surface utilities and library.

class nitrogen.pes.CFOUR(atomic_symbols, params, work_dir='./scratch', cleanup=True, units='angstrom', sections=[])

A simple CFOUR interface for accessing single point energies and derivatives.

params

The CFOUR keyword values used for the calculation.

Type

dict

natoms

The number of atoms. This equals nx/3.

Type

int

atomic_symbols

The atomic symbols of each atom.

Type

list

work_dir

The work directory.

Type

str

cleanup

The work directory clean-up flag.

Type

bool

Create a new CFOUR interface.

atomic_symbolslist

The atomic symbols of each atom.

paramsdict

CFOUR keyword and value pairs.

work_dirstr, optional

The path to the work directory. The default is ‘.scratch’.

cleanupbool, optional

If True (default), work directories will be deleted after use.

units{‘angstrom’, ‘bohr’}, optional

The Cartesian units.

sectionslist of str, optional

A list of additional keyword sections. Each section is a single string with appropriate newline characters. “

” is written to the ZMAT

file after each section.

All elements in params will be added as keywords to the *CFOUR() section of a CFOUR ZMAT input file. The COORD, UNITS, DERIV_LEV, PRINT, and VIB keywords are handled automatically. These should not be supplied by the user. Only methods for which VIB=ANALYTIC (i.e. analytic Hessians) is available are supported. For other methods, use lower CFOUR derivatives with NITROGEN finite difference drivers.

nitrogen.pes.curvVib(Q0, pes, cs, masses, mode='bodyframe', fidx=0)

Calculate curvilinear vibrational normal coordinates and frequencies at a stationary point.

Parameters

  • Q0 (array_like) – The stationary point coordinates.

  • pes (DFun) – Potential energy surface function.

  • cs (CoordSys) – Coordinate system.

  • masses (array_like) – Coordinate system masses.

  • mode ({'bodyframe'}) – Coordinate frame mode. ‘bodyframe’ treats the cs coordinate system as the body-fixed frame. ‘bodyframe’ is the default.

  • fidx (integer, optional) – The DFun function index to use. The default is 0

Returns

  • omega (ndarray) – The harmonic frequencies (times \(\hbar\)). Negative frequencies are returned for imaginary frequencies.

  • nctrans (LinearTrans) – The normal coordinate transformation. The columns of nctrans.T are the (right) eigenvectors of the GF matrix, scaled to the dimensionless normal coordinates.

Notes

The frequency calculation requires the inverse of the molecular moment of inertia tensor. At linear geometries, this tensor is singular. Instead, its pseudo-inverse is calculated, which gives the correct results at linear geometries.

nitrogen.pes.loadpes(pesname)

Load built-in PES.

Parameters

pesname (str) – PES name.

Returns

pes – The PES as a DFun object.

Return type

DFun

nitrogen.pes.opt_bfgs(F, x0, Hinit=None, alpha=1.0, tol=1e-06, disp=False, fidx=0, var=None, maxiter=100)

Quasi-Newton optimization with a Broyden-Fletcher-Goldfarb-Shannon (BFGS) Hessian update.

Parameters

  • F (DFun) – A DFun object with second derivative support.

  • x0 (array_like) – The initial coordinate values.

  • Hinit (array_like, optional) – The initial Hamiltonian. If None, the initial Hessian will be calculated explicitly. If array_like, then a user-supplied array is used. Note that Hinit must be supplied with respect to only the variables in var and in that order.

  • alpha (float, optional) – The update step scaling coefficint. The default is 1.0.

  • tol (float, optional) – The gradient norm tolerance. The default is 1e-6.

  • disp (bool, optional) – Print detailed output. The default is False.

  • fidx (int, optional) – The DFun function index to optimize. The default is 0.

  • var (list of int, optional) – The variables to optimize. All others will be constrained to the values in x0. If None (default), all coordinates will be optimized

  • maxiter (int, optional) – The maximum number of steps. The default is 100.

Returns

  • xopt (ndarray) – The optimized coordinates

  • Fopt (scalar) – The value of the optimized function

nitrogen.pes.opt_newton(F, x0, alpha=1.0, tol=1e-06, disp=False, fidx=0, var=None, maxiter=100)

Simple Newton-Raphson optimization.

Parameters

  • F (DFun) – A DFun object with second derivative support.

  • x0 (array_like) – The initial coordinate values.

  • alpha (float, optional) – The update step scaling coefficint. The default is 1.0.

  • tol (float, optional) – The gradient norm tolerance. The default is 1e-6.

  • disp (bool, optional) – Print detailed output. The default is False.

  • fidx (int, optional) – The DFun function index to optimize. The default is 0.

  • var (list of int, optional) – The variables to optimize. All others will be constrained to the values in x0. If None (default), all coordinates will be optimized

  • maxiter (int, optional) – The maximum number of steps. The default is 100.

Returns

  • xopt (ndarray) – The optimized coordinates

  • Fopt (scalar) – The value of the optimized function

fit.py

Surface fitting tools

class nitrogen.pes.fit.BondVectorPIP(yfun, vfun, p)

Bases: nitrogen.dfun.DFun

A general PIP expansion for dipole or other vector functions.

Parameters

  • yfun (DFun) – The internuclear distance function.

  • vfun (DFun) – The bond vector function.

  • p ((nvec,n) ndarray) – The expansion coefficients of each bond vector function.

class nitrogen.pes.fit.BondVectorR(n, Vij)

Bases: nitrogen.dfun.DFun

Bond vector function for linear distance

\[v_{ij} = X_j - X_i\]
n

The number of atoms.

Type

integer

Vij

The bond pairs

Type

list of (2,)

Parameters

  • n (integer) – The number of atoms

  • Vij (list of (2,)) – The bond pairs

class nitrogen.pes.fit.ExpPIP(a, natoms, coeff, terms, input_fun=None)

Bases: nitrogen.dfun.DFun

A Cython implementation of a simple PIP w.r.t Morse-scaled internuclear distances.

See also

fitSimplePIP

linear least squares PIP fitting

Create a ExpPIP object

Parameters

  • a (float) – The Morse parameter

  • natoms (integer) – The number of atoms

  • coeff ((nt,) ndarray) – The expansion coefficients

  • terms ((nt,ny) ndarray) – The (unsorted) expansion terms, where ny = natoms`(`natoms-1)/2

Notes

coeff will be cast to np.float64 terms will be cast to np.int32

class nitrogen.pes.fit.ExpPIPTerms(a, natoms, terms, input_fun=None)

Bases: nitrogen.dfun.DFun

Term-by-term of ExpPIP

See also

ExpPIP

Create a ExpPIP object

Parameters

  • a (float) – The Morse parameter

  • natoms (integer) – The number of atoms

  • terms ((nt,ny) ndarray) – The (unsorted) expansion terms, where ny = natoms`(`natoms-1)/2

Notes

terms will be cast to np.int32

class nitrogen.pes.fit.InternuclearExp(n, a, r0=0.0, offset=0.0)

Bases: nitrogen.dfun.DFun

Internuclear coordinate function for exponential (Morse) scaled distance

\[y_{ij} = \exp[-(r_{ij} - r_0) / a] - \mathrm{offset}\]

The \(y_{ij}\) functions are ordered \(y_{12}, y_{13}, \ldots, y_{1n}, y_{23}, \ldots, y_{n-1,n}\).

a

The exponential length parameter

Type

float

n

The number of atoms.

Type

integer

r0

The reference length.

Type

float

offset

The offset of each coordinate function.

Type

float

Parameters

  • n (integer) – The number of atoms

  • a (float) – The exponential length parameter.

  • r0 (float, optional) – A reference length. The default is 0.

  • offset (float, optional) – Offset of each coordinate function. The default is 0.

class nitrogen.pes.fit.InternuclearR(n)

Bases: nitrogen.dfun.DFun

Internuclear coordinate function for linear distance

n

The number of atoms.

Type

integer

Parameters

n (integer) – The number of atoms

nitrogen.pes.fit.Sn(n, indices)

Return the permutations of identical particles.

Parameters

  • n (integer) – The total number of particles

  • indices (array_like) – A list of identical indices (0 through n-1).

Returns

P – A list of permutations

Return type

list

nitrogen.pes.fit.atom2vperm(P, Vij)

Convert a list of permutations P of atoms to a list of permutation of bond-pair vectors, including the sign of the permutation.

Parameters

  • P (list) – The atomic permutations.

  • Vij (list of (2,)) – The bond pairs.

Returns

  • Pv (list) – The bond vector permutations

  • Sv (list) – The sign of the permutations.

nitrogen.pes.fit.atom2yidx(n)

Calculate the 2d table mapping pairs of atomic indices to y-variable indices.

Parameters

n (int) – The atom count.

Returns

yidx – The map of atom indices to y indices. The map is symmetric, y[i,j] = y[j,i], and only the off-diagonal elements are valid.

Return type

ndarray

nitrogen.pes.fit.atom2yperm(P)

Convert a list of permutations P of atoms to a list of permutation of y variables.

Parameters

P (list) – The atomic permutations.

Returns

Py – The y-variable permutations

Return type

list

nitrogen.pes.fit.count_idx_bodies(idx, natoms)

Count the number of atoms involved in a single monomial

Parameters

  • idx ((ny,) array_like) – The single monomial powers

  • natoms (integer) – The total number of atoms

Returns

The number of atoms involved.

Return type

integer

Notes

The \(y_{ij}\) internuclear functions are ordered \(y_{12}, y_{13}, \ldots, y_{1n}, y_{23}, \ldots, y_{n-1,n}\), where \(n\) = natoms. #

nitrogen.pes.fit.fitFourier(x, y, max_freq, period=None, symmetry=None)

Fit a Fourier series.

Parameters

  • x ((N,) array_like) – The argument at N sampling points.

  • y ((N,) or (m,N) array_like) – The values of one or more (m) functions.

  • max_freq (integer) – The maximum frequency harmonic.

  • period (float, optional.) – The period of x. If None (default), period = \(2\pi\) is assumed.

  • symmetry (integer or (m,) array_like, optional) – The symmetry type of each function. If None (default), no symmetry is assumed. See Notes.

Returns

The Fourier coefficients

Return type

(n,) or (m,n) ndarray

Notes

The expansion coefficients are defined as

\[f(x) = c_0 + c_1 \sin \sigma x + c_2 \cos \sigma x + c_3 \sin 2 \sigma x + c_4 \cos 2 \sigma x + \cdots\]

where \(\sigma = 2\pi/\)period.

The symmetry keyword specifies a constraint on the Fourier series for each fitted function. A value of 0 fits all terms, 1 fits only cosine terms, 2 fits only sine terms, and -1 fixes all parameters to 0.

See also

nitrogen.dfun.FourierSeries()

A Fourier series DFun.

nitrogen.pes.fit.fitNBodyExpPIP(X, F, P, degree, a, n, Xscale)

Fit an exponential-permutationally invariant polynomial with n-body constraints.

Parameters

  • X ((3*n,N) ndarray) – The Cartesian coordinates of n particles at N sampling points.

  • F ((nd,N) ndarray) – The derivative array of the function values. If only the value is fitted, for example, nd is 1.

  • P (list) – The permutation elements. Each element of P is a permutation of [0,1,2,...,n-1].

  • degree (integer) – The polynomial degree.

  • a (float) – The Morse-exponential scale parameter, \(y = \exp(-r/a)\).

  • n (integer) – The maximum number of atoms per term.

  • Xscale (float) – The Cartesian length scale. Its powers will be used to scale derivatives to the value units.

Returns

  • p (ndarray) – The expansion coefficients.

  • terms (ndarray) – The expansion terms.

  • res (ndarray) – The scaled residuals.

See also

fitSimplePIP

nitrogen.pes.fit.fitSimplePIP(X, F, P, yfun, degree, Xscale)

Fit a permutationally invariant polynomial in terms of an internuclear distance function.

Parameters

  • X ((3*n,N) ndarray) – The Cartesian coordinates of n particles at N sampling points.

  • F ((nd,N) ndarray) – The derivative array of the function values. If only the value is fitted, for example, nd is 1.

  • P (list) – The permutation elements. Each element of P is a permutation of [0,1,2,...,n-1].

  • yfun (DFun) – The internuclear coordinate function. It takes as its argument the (3*n,…) Cartesian coordinates and returns the n*(n-1)//2 internuclear functions. These functions need to be consisten with the permutational symmetries implied by P.

  • degree (integer) – The polynomial degree.

  • Xscale (float) – The Cartesian length scale. Its powers will be used to scale derivatives to the value units.

Returns

  • p (ndarray) – The expansion coefficients

  • res (ndarray) – The scaled residuals.

  • Ffun (DFun) – The fitted function, F(X).

Notes

Derivatives up to arbitrary order can be fitted.

See also

InternuclearExp

internuclear function for exponentially scaled distance

InternuclearR

internuclear function for linear distance

nitrogen.pes.fit.fitSimplePIPDipole(X, D, P, Vij, yfun, vfun, degree, Xscale)

Fit a permutationally invariant polynomial surface for a dipole moment function.

Parameters

  • X ((3*n,N) ndarray) – The Cartesian coordinates of n particles at N sampling points.

  • D ((nd,3,N) ndarray) – The dipole moment derivative array in the same frame as X at each point.

  • P (list) – The permutation elements. Each element of P is a permutation of [0,1,2,...,n-1].

  • Vij (list of (2,)) – The list of bond vector atom-pairs, e.g., [(0,1), (0,2)].

  • yfun (DFun) – The internuclear coordinate function.

  • vfun (DFun) – The bond vector function.

  • degree (integer) – The polynomial degree.

  • Xscale (float) – The Cartesian length scale. Its powers will be used to scale derivatives to the value units.

Returns

  • p ((nvec, nterms) ndarray) – The expansion coefficients of each bond-vector expansion.

  • res (ndarray) – The scaled residuals.

  • D_function (DFun) – The fitted dipole function, D(X).

Notes

Derivatives up to arbitrary order can be fitted.

See also

InternuclearExp

internuclear function for exponentially scaled distance

InternuclearR

internuclear function for linear distance

nitrogen.pes.fit.ny2natoms(ny)

Calculate the number of atoms given the number of internuclear functions

Parameters

ny (integer) – The number of internuclear functions

Returns

The number of atoms

Return type

integer

nitrogen.pes.fit.productP(P1, P2)

Return the direct product of all permutations in P1 and P2.

Parameters

  • P1 (list) – A list of permutations

  • P2 (list) – A list of permutations

Returns

P – The direct products

Return type

list

nitrogen.pes.fit.trimNBody(idx_list, n)

Trim a list of monomial terms of internuclear functions to those involving n or few atoms,

Parameters

  • idx_list ((nt,ny) array_like) – The list of monomial powers

  • n (integer) – The maximum number of simultaneous bodies.

Returns

trimmed_list – The trimmed list

Return type

(ns,ny) array_like

rxnpath.py

Reaction path routines

class nitrogen.pes.rxnpath.CoriolisBlock(cs, masses=None)

Bases: nitrogen.dfun.DFun

The rovibrational (Coriolis) block of the inverse metric.

mode = ‘bodyframe’ is assumed.

Parameters

  • cs (CoordSys) – The coordinate system.

  • masses (array_like, optional) – The masses. If None, unit masses are assumed.

class nitrogen.pes.rxnpath.InverseInertia(cs, masses=None)

Bases: nitrogen.dfun.DFun

An inverse effective inertia tensor.

The derivatives of the lower triangle in packed storage of the inverse effective inertia tensor is calculated.

mode = ‘bodyframe’ is assumed.

Parameters

  • cs (CoordSys) – The coordinate system.

  • masses (array_like, optional) – The masses. If None, unit masses are assumed.

class nitrogen.pes.rxnpath.InverseMetric(cs, masses=None, mode='bodyframe', planar_axis=None)

Bases: nitrogen.dfun.DFun

An inverse (vibrational) metric function.

The derivatives of the lower triangle in packed storage of the inverse metric, or the vibrational block of the inverse metric for bodyframe embedding, is calculated.

Parameters

  • cs (CoordSys) – The coordinate system.

  • masses (array_like, optional) – The masses. If None, unit masses are assumed.

  • mode ({'bodyframe'}, optional) – The embedding mode.

  • planar_axis ({None,0,1,2}, optional) – The normal axis for linear/planar coordinate systems. If None, this is not used. See Notes for more details.

Notes

The planar_axis parameter is used to avoid indeterminancies in the inverse metric at linear geometries. In this case, only strictly planar coordinate systems should be used. The block of the metric tensor for the two in-plane axes (which is singular) decouples from the rest of the metric and can be ignored for calculating the vibrational block of the inverse metric.

nitrogen.pes.rxnpath.LQA_nonstationary(q0, V, G, arclength='massweighted', proxy_index=0)

Compute the covariant gradient (i.e., the reaction path tangent vector) within a local quadratic approximation at non-stationary points.

Parameters

  • q0 ((nq,...) array_like) – The non-stationary evaluation points.

  • V (DFun) – The potential energy surface

  • G (DFun) – The inverse metric. If None, a constant unit metric is assumed.

  • arclength ({'massweighted','gradient','proxy'}, optional) – The path parameterization convention. The default is ‘massweighted’.

  • proxy_index (integer, optional) – The proxy coordinate index.

Returns

  • w0 ((nq,…) ndarray) – The path tangent at q0.

  • W ((nq,nq,…) ndarray) – The path tangent-gradient at q0. W[i,j] is equal to \(\partial_j w_i \vert_0\).

  • G0 ((nq,nq,…) ndarray) – The inverse metric evaluated at q. The inner product of tangent vectors (not gradient vectors) with G0 yields their proper 2-norm.

Notes

The path tangent is defined as

\[\mathbf{w} = -G\mathbf{f}/h,\]

where \(\mathbf{f}\) is the gradient, \(G\) is the inverse metric tensor, and \(h\) defines the path parameterization normalization.

The LQA approximates the local path tangent as

\[\mathbf{w} \approx \mathbf{w}_0 + W(\mathbf{q} - \mathbf{q}_0)\]

For arclength == 'massweighted', the natural mass-weighted arc length is used, \(h = (\mathbf{f}^T G \mathbf{f})^{1/2}\).

For arclength == 'gradient', \(h = 1\).

For arclength == 'proxy', \(h = -(G\mathbf{f})_i\), where \(i\) is the coordinate index specified by proxy_index. In this case, the path is parameterized by one of the coordinates themselves – the ‘’proxy’’ coordinate – so \(w_i = 1\).

nitrogen.pes.rxnpath.christoffel_symbol(q, G, kind='first')

Calculate the Christoffel symbols, \(\Gamma_{ijk}\) or \(\Gamma^i{}_{jk}\).

Parameters

  • q ((nq,...) array_like) – The evaluation points.

  • G (DFun) – The inverse metric tensor, \(G^{ij}\).

  • kind ({'first','second'}, optional) – Calculate symbols of the first or second kind. ‘first’ is the default.

Returns

Gamma – The Christoffel symbols.

Return type

(nq,nq,nq,…) ndarray

Notes

The Christoffel symbols of the first kind are

\[\Gamma_{ijk} = \frac{1}{2}\left( \partial_k g_{ij} + \partial_j g_{ik} - \partial_i g_{jk} \right)\]

The Christoffel symbols of the second kind are

\[\Gamma^i{}_{jk} = G^{im} \Gamma_{mjk}\]

Either kind is symmetric in the last two indices.

nitrogen.pes.rxnpath.correct_vib_order(omega, T, Tl, hbar=None)

Correct the ordering and phase of reaction path vibrational modes.

Parameters

  • omega ((nq,N) ndarray) – The harmonic frequncies (in energy units.)

  • T ((nq,nm,N) ndarray) – The displacement vectors normalized as reduced dimensionless normal modes.

  • Tl ((nq,nm,N) ndarray) – The left-hand displacement vectors (i.e. T multiplied by the effective metric tensor).

  • hbar (float, optional) – The value of \(\hbar\). If None, the default NITROGEN units will be used.

Returns

  • omega_new ((nq,N) ndarray) – The ordered frequencies

  • T_new ((nq,nm,N) ndarray) – The ordered displacement vectors.

  • Tl_new ((nq,nm,N) ndarray) – The ordered left-hand vectors.

See also

pathvib_nonstationary

Calculates omega, T, and Tl.

nitrogen.pes.rxnpath.covariant_hessian(q, V, G)

Calculate the covariant Hessian tensor, \(\nabla_i \nabla_j V\).

Parameters

  • q ((nq,...)) – The evaluation points.

  • V (DFun) – The potential energy surface.

  • G (DFun) – The inverse metric tensor.

Returns

  • f ((nq,…)) – The gradient.

  • H ((nq,nq,…)) – The covariant Hessian.

Notes

The covariant gradient is equal to the regular gradient, \(\nabla_i V = \partial_i V\). The covariant Hessian is

\[ \begin{align}\begin{aligned}H_{ij} &= \nabla_i \nabla_j V \\\ &= \nabla_i f_j \\\ &= \partial_i \partial_j V - f_k \Gamma^k{}_{ij}\end{aligned}\end{align} \]

The covariant Hessian is symmetric, \(H_{ij} = H_{ji}\).

nitrogen.pes.rxnpath.cubic_spline(x, y, boundary='natural', boundary_value=(0.0, 0.0), return_jacobian=False)

Calculate the parameters for a cubic spline.

Parameters

  • x (1d array_like) – The node points.

  • y (1d array_like) – The function vales.

  • boundary ({'natural', 'notaknot', 'first', 'second'}, optional) – The boundary condition type. See Notes.

  • boundary_value ((2,) tuple, optional) – The boundary condition values, if applicable. See Notes.

  • return_jacobian (bool, optional) – Also return the derivatives of the spline parameters with respect to the function values y.

Returns

  • c ((n-1,4) ndarray) – The cubic spline parameters.

  • dc ((n,n-1,4) ndarray) – The parameter Jacobian. dc[i] is the derivative of the spline parameters with respect to y[i]. Only returned if return_jacobian is True.

Notes

The cubic spline consists of a cubic polynomial in each of the \(n-1\) regions between the \(n\) node points for a total of \(4n-4\) parameters. Matching the node values, as well as enforcing continuity of the first and second derivatives at internal nodes, yields only \(4n-6\) constraints and leaves two more to be chosen. The boundary parameter determines these two conditions.

A value of 'natural' sets the second derivatives of the spline at the endpoints to 0.

A value of 'notaknot' uses the not-a-knot condition. The third derivative at the first and last interior nodes is continuous.

A value of ‘first’ sets the first derivatives at the boundaries equal to the values passed in boundary_value.

A value of ‘second’ sets the second derivatives at the boundaries equal to the values passed in boundary_value.

nitrogen.pes.rxnpath.cubic_spline_derivative(x, c, x0)

Evaluate the derivative of a cubic spline

Parameters

  • x (array_like) – The evaluation points.

  • c ((n-1,4)) – The spline parameters.

  • x0 (array_like) – The ordered spline nodes.

Returns

dy – The spline derivatives at x.

Return type

ndarray

nitrogen.pes.rxnpath.cubic_spline_val(x, c, x0)

Evaluate a cubic spline

Parameters

  • x (array_like) – The evaluation points.

  • c ((n-1,4) ndarray) – The spline parameters.

  • x0 (array_like) – The ordered spline nodes.

Returns

y – The spline values at x.

Return type

ndarray

nitrogen.pes.rxnpath.invertderiv(f, x0=0.0)

Calculate the derivatives of the inverse function given the derivatives of a function.

Parameters

  • f ((deriv+1,...) ndarray) – The derivative array of the original function, \(f(x)\).

  • x0 ((...) ndarray or scalar, optional) – The expansion value of \(x\). The default is zero.

Returns

F – The deriative array of the inverse function, \(F(y)\).

Return type

(deriv+1,…) ndarray

nitrogen.pes.rxnpath.pathderivchain(A, X)

Calculate the derivative of a quantity with respect to the path parameter via the multivariate chain rule.

Parameters

  • A ((nd,...',...) ndarray) – The derivative array of \(A(x)\) with respect to nvar \(x\) coordinates.

  • X ((deriv+1,nvar,...) ndarray) – The derivative array of nvar :math:`x coordinates with respect to the path parameter \(s\).

Returns

B – The derivative array of \(B(s) = A(x(s))\)

Return type

(deriv+1,…’,…) ndarray

nitrogen.pes.rxnpath.pathvib_nonstationary(q, V, G, hbar=None, fproj=None)

Calculate reaction path normal modes using the orthogonally projected covariant Hessian.

Parameters

  • q ((nq,...) array_like) – The evaluation points at non-stationary geometries.

  • V (DFun) – The potential energy surface.

  • G (DFun) – The inverse metric tensor.

  • hbar (float, optional) – The value of \(\hbar\). If None, the default NITROGEN units will be used.

  • fproj ((nq,...) array_like, optional) – The projection vector at each point. If None (default), the gradient will be used.

Returns

  • omega ((nq-1,…) ndarray) – The orthogonal path frequencies in energy units.

  • T ((nq,nq-1,…) ndarray) – The displacement vectors of the nq-1 orthogonal modes, normalized as reduced dimensionless normal coordinates. T[i,j] is the displacement of coordinate i for unit amplitude of normal mode j.

  • Tl ((nq,nq-1,…) ndarray) – The lowered-index transformation of T.

nitrogen.pes.rxnpath.proxyderiv(qn, star_index)

Convert path derivatives to those with respect to one of the coordinates

Parameters

  • qn ((deriv+1,n,...)) – The derivative array of the n coordinates with respect to a path parameter

  • star_index (integer) – The index of the proxy coordinate

Returns

qn_star – The derivatives of the n coordinates with respect to the proxy coordinate

Return type

(deriv+1,n,…)

nitrogen.pes.rxnpath.qderiv_nonstationary(q0, deriv, V, G, direction='descend')

Calculate the reaction path deriatives evaluated at an arbitrary (non-stationary) point.

Parameters

  • q0 ((nq,...) array_like) – The evaluation points

  • deriv (integer) – The maximum derivative order to calculate.

  • V (DFun) – The potential energy surface.

  • G (DFun) – The inverse metric tensor.

  • direction ({'descend', 'ascend'}) – The direction of the path coordinate. If ‘descend’, then the path follows the negative gradient for increasing arc length.

Returns

q – The derivatives of the reaction path with respect to the arc length parameter.

Return type

(deriv + 1, nq, …)

nitrogen.pes.rxnpath.qderiv_stationary(q0, deriv, V, G, direction='normal')

Calculate the reaction path deriatives evaluated at a stationary point.

Parameters

  • q0 (array_like) – The evaluation point.

  • deriv (integer) – The maximum derivative order to calculate.

  • V (DFun) – The potential energy surface.

  • G (DFun) – The inverse metric tensor. If None, identity is assumed

  • direction ({'normal', 'reverse'}, optional) – The direction of the path coordinate. If ‘normal’, the sign of the path tangent is determined by making its largest element positive. If ‘reverse’, the sign is reversed.

Returns

q – The derivative of the reaction path with respect to the arc length parameter.

Return type

(deriv + 1, nq) ndarray

nitrogen.pes.rxnpath.singleLeibniz(a1, a2)

A simple single-variable Leibniz product rule

a1,a2(deriv+1,…) ndarray

Single variable deriative arrays

nitrogen.pes.rxnpath.spline_proxy_path(V, G, q0, q1, proxy_index, nodes, match_level=0, is_stat=(True, True), max_iter=20, deltarms=1e-06)

Compute a reaction path between two points as a spline function with respect to a proxy coordinate.

Parameters

  • V (DFun) – The potential energy surface.

  • G (DFun) – The inverse metric tensor.

  • q0 ((nq,)) – A path end-point.

  • q1 ((nq,)) – A path end-point.

  • proxy_index (integer) – The proxy coordinate index.

  • nodes (integer) – The number of interior spline nodes.

  • match_level (integer, optional) – The exact boundary condition constraint level. If match_level = 0, the path end-points are constrained. If match_level = 1, the path tangent is also constrained. The default is 0.

  • is_stat ((2,) tuple of boolean, optional) – Specifies whether each end-point is a stationary point. The default is (True,True). This only matters if match_level is greater than 0.

  • max_iter (integer, optional) – The maximum number of path updates.

  • deltarms (float) – The threshold change to the estimated path rms error per step.

Returns

  • path_list (list of (nq,nodes) ndarray) – The path nodes for each update.

  • spline_list (list of (nq, nodes-1, 4) ndarray) – The cubic spline parameters for each update.

  • rms_list (list of float) – The path rms for each update.

Notes

An initial guess of the spline nodes is formed by a minimial polynomial interpolantion that meets the matching conditions (a linear interpolant for match_level = 0, and a cubic interpolant for match_level = 1). The positions of the spline nodes are updated by a Newton optimization scheme based on a local quadratic approximation to the path tangent near each node point.

See also

cubic_spline_val

Evaluate the cubic spline functions.

nitrogen.pes.rxnpath.symfull_axis(A, axis=0)

Expand packed matrix to full symmetric matrix.

Parameters

  • A (ndarray) – The packed array

  • axis (integer, optional) – The packed axis. The default is 0.

Returns

Afull – A new array with the packed axis expanded into two symmetric axes.

Return type

ndarray

lapath.py

Least-action paths

nitrogen.pes.lapath.calcPathAction(qpath, V, cs, masses, band_action=1000.0, kinetic_energy=0.0, deriv=0, Vmin=0, ignore_g_deriv=False)

Calculate the path action and its derivatives.

Parameters

  • qpath ((nq,N) ndarray) – The coordinate path of nq coordinates along N nodes.

  • V (DFun) – The potential energy surface supporting deriv >= 1.

  • cs (CoordSys) – The coordinate system

  • masses (array_like) – The masses.

  • band_action (float, optional) – The band action force constant. The default is 1000.0. This parameter has units of [action] / [arc length]^2

  • kinetic_energy (float, optional) – An energy offset added to the potential. The default is 0. This parameter has units of [energy].

  • deriv ({0,1,2}, optional) – The derivative order to calculate. The default is 0.

  • Vmin (float, optional) – The minimum energy to subtract. The default is 0.

  • ignore_g_deriv (bool, optional) – If True, derivatives of the metric tensor are ignored. The default is False.

Returns

  • s (float) – The path arc length.

  • I (float) – The path action.

  • B (float) – The elastic band action.

  • If deriv >= 1

  • Ds ((nq,N) ndarray) – The derivative of s with respect to the path node parameters.

  • DI ((nq,N) ndarray) – The derivative of I with respect to the path node parameters.

  • DB ((nq,N) ndarray) – The derivative of B with respect to the path node parameters.

  • If deriv >= 2

  • D2s ((nq,N,nq,N) ndarray) – The Hessian of s with respect to the path node parameters.

  • D2I ((nq,N,nq,N) ndarray) – The Hessian of I with respect to the path node parameters.

  • D2B ((nq,N,nq,N) ndarray) – The Hessian of B with respect to the path node parameters.

Notes

The path action is approximated by a sum over sequential linear segments between nodes in coordinate space.

The length of each path segment is calculated as

\[\delta s_i = \sqrt{ \delta q_i \frac{\tilde{g}_i + \tilde{g}_{i+1}}{2} \delta q_i },\]

where \(\delta q_i = (q_{i+1} - q_i)\) is the path segment displacement and \(\tilde{g}\) is the effective vibrational metric tensor (i.e. the inverse of the vibrational block of the full inverse metric tensor \(G\).)

The effective momentum at each node is \(p_i = \sqrt{2(V(q_i) + E)}\), where \(E\) is the kinetic energy offset. The total action is then

\[I = \sum_i \frac{p_i + p_{i+1}}{2} \delta s_i\]

The arc length and action units are defined by the V, cs, and masses parameters.

An elastic band action is also calculated as

\[B = \sum_i b (\delta s_i - \delta s_{i+1})^2,\]

where b is the band action force constant.

Note that the minimum energy of V will be subtracted before computing the action and associated derivatives.

This implementation ignores the kinetic pseudo-potential contribution, which may become significant if the path approaches singular points of the metric tensor.

nitrogen.pes.lapath.calcWKBtunneling(qpath, V, cs, masses, disp=True, hbar=None)

Perform a simple WKB tunneling analysis on a 1D path.

Parameters

  • qpath ((nq,N)) – The path between two equivalent minima.

  • V (DFun or function) – The potential energy

  • cs (CoordSys) – The coordinate system.

  • masses (array_like) – The masses.

  • disp (bool, optional) – Display/print flag. The default is True.

  • hbar (float, optional) – The value of hbar. If None, default NITROGEN units are used.

Returns

  • omega (float) – The estimated harmonic energy

  • phi (float) – The imaginary action in units of \(\hbar\).

  • dE (float) – The tunneling energy splitting

nitrogen.pes.lapath.opt_path_bfgs(qpath, V, cs, masses, band_action=1000.0, kinetic_energy=0.0, Hinit=None, Vmin=0.0, alpha=1.0, tol=1e-06, maxiter=100, disp=True, streak=None)

Optimize a least-action path using BFGS Hessian updates.

Parameters

  • qpath ((nq,N) ndarray) – The initial coordinate path of nq coordinates along N nodes.

  • V (DFun) – The potential energy surface supporting deriv >= 1.

  • cs (CoordSys) – The coordinate system

  • masses (array_like) – The masses.

  • band_action (float, optional) – The band action force constant. The default is 1000.0.

  • kinetic_energy (float, optional) – An energy offset added to the potential. The default is 0.

  • Vmin (float, optional) – The minimum potential energy. The default is 0.

  • alpha (float, optional) – The update step scaling coefficint. The default is 1.0.

  • tol (float, optional) – The gradient norm tolerance. The default is 1e-6.

  • maxiter (int, optional) – The maximum number of steps. The default is 100.

  • disp (bool, optional) – Print detailed output. The default is True.

Returns

  • qopt ((nq,N) ndarray) – The optimized path

  • I (float) – The optimized path action

  • s (float) – The optimized path length

  • B (float) – The optimized band action