Hamiltonian Computational Chemistry: Geometrical Structures in Chemical Dynamics and Kinetics

TL;DR

This paper explores the shared geometrical structures in classical mechanics, quantum mechanics, and thermodynamics to study chemical dynamics and kinetics using Hamiltonian methods.

Contribution

The paper introduces a unified Hamiltonian framework for chemical dynamics and kinetics using geometrical structures and high-order solvers.

Findings

Physical states in integrable systems are represented by Lagrangian submanifolds in phase space.

High-order finite-difference algorithms are accurate but resource-intensive for systems with many degrees of freedom.

Hamiltonian neural networks show promise for solving Hamilton’s equations in complex systems.

Abstract

The common geometrical (symplectic) structures of classical mechanics, quantum mechanics, and classical thermodynamics are unveiled with three pictures. These cardinal theories, mainly at the non-relativistic approximation, are the cornerstones for studying chemical dynamics and chemical kinetics. Working in extended phase spaces, we show that the physical states of integrable dynamical systems are depicted by Lagrangian submanifolds embedded in phase space. Observable quantities are calculated by properly transforming the extended phase space onto a reduced space, and trajectories are integrated by solving Hamilton’s equations of motion. After defining a Riemannian metric, we can also estimate the length between two states. Local constants of motion are investigated by integrating Jacobi fields and solving the variational linear equations. Diagonalizing the symplectic fundamental matrix, eigenvalues equal to one reveal the number of constants of motion. For conservative systems, geometrical quantum mechanics has proved that solving the Schrödinger equation in extended Hilbert space, which incorporates the quantum phase, is equivalent to solving Hamilton’s equations in the projective Hilbert space. In classical thermodynamics, we take entropy and energy as canonical variables to construct the extended phase space and to represent the Lagrangian submanifold. Hamilton’s and variational equations are written and solved in the same fashion as in classical mechanics. Solvers based on high-order finite differences for numerically solving Hamilton’s, variational, and Schrödinger equations are described. Employing the Hénon–Heiles two-dimensional nonlinear model, representative results for time-dependent, quantum, and dissipative macroscopic systems are shown to illustrate concepts and methods. High-order finite-difference algorithms, despite their accuracy in low-dimensional systems, require substantial computer resources when they are applied to systems with many degrees of freedom, such as polyatomic molecules. We discuss recent research progress in employing Hamiltonian neural networks for solving Hamilton’s equations. It turns out that Hamiltonian geometry, shared with all physical theories, yields the necessary and sufficient conditions for the mutual assistance of humans and machines in deep-learning processes.

Figures

Full text

1. Introduction

Chemistry is now deeply rooted in the two fundamental physical theories, quantum and classical mechanics. Quantum chemistry and molecular dynamics computer programs are indispensable devices in almost all kinds of chemical research. Thermodynamics, also a vital theory for chemistry that lacked a genuine mathematical foundation for years, has acquired the necessary mathematical framework in the last two decades, which brings it to the same level as classical and quantum mechanics. On the other hand, the advancements in mathematics that occurred in the 20th century in the fields of differential geometry and topology have revealed common geometrical structures in all physical theories. These discoveries provide a deeper understanding of the physical theories per se and pave the way for their application in other scientific fields, such as chemistry.

The purpose of this article is to render an introductory and graphical presentation of common geometrical structures of the principal physical theories and the consequences they may have in chemistry, especially via their numerical applications. Specifically, after almost two centuries of evolution of Hamiltonian theory, a modern geometrical description emerged, and significant theorems and techniques for locating invariant structures in phase space, and thus constants of motion, have been found [1,2,3]. Chemical dynamics and spectroscopy have tremendously benefited from applying these methods to comprehend the behaviors of highly excited molecules [4,5].

It is essential to underline that the stage of action for molecular dynamics, within a classical mechanical approach, is the phase space and its tangent bundle, which means that both generalized coordinates and conjugate momenta should be taken into account. We also emphasize that in the 1980s, it was found that quantum and classical mechanics share some common geometrical properties, which explain similarities as well as differences. It was shown that instead of the usual algebraic linear quantum theory formulated either in the Schrödinger or Heisenberg picture, we can take the quantum analog of phase space to be the projective Hilbert space of the extended Hilbert space [6,7,8].

Earlier, in the 1970s, it was recognized that the mathematical framework for thermodynamics is the contact geometry [9,10] of the physical states embedded in an odd-dimensional state space [11,12,13,14,15,16,17]. However, at the beginning of the twenty-first century, Balian and Valentin [18] made a significant contribution by publishing a Hamiltonian theory of thermodynamics in an extended phase space. Based on Callen’s [19] formulation of thermodynamics, they produced a Hamiltonian theory for reversible and irreversible processes equivalent to classical mechanics. They studied homogeneous Hamiltonians with generalized coordinates, the complete set of extensive properties (i.e., entropy, internal energy, volume, particle numbers, etc.) and conjugate momenta proportional to intensive properties (i.e., temperature, pressure, chemical potential, etc.), either in the entropy or energy representation. Gibbs’s fundamental equation was employed to describe the physical state manifold. This work triggered numerous studies, the results of which have revealed several geometrical properties common to classical and quantum mechanics [20,21,22,23,24].

It turns out that the mathematical abstraction and generalization of geometrical Hamiltonian theory in extended phase space lead to a common computational platform for working in both microscopic and macroscopic worlds. The aim of the present article is to highlight that Hamiltonian theories share some common geometrical properties in extended phase space with classical mechanics, quantum mechanics, and thermodynamics, which manifest the foundations for formulating and comprehending chemical dynamics and chemical kinetics.

In Section 2, we present, with the help of two pictures, the Hamiltonian theories of classical and quantum mechanics. Similarly, in Section 3, we discuss Hamiltonian thermodynamics in contemporary mathematical language that unveils the geometrical properties of this theory. In Section 4, we describe numerical algorithms for solving Hamilton’s equations of motion and variational equations, common to the three principal theories. We mainly focus on high-order finite-difference (FD) methods and their relation to pseudospectral methods (PS) [25]. Furthermore, to illustrate the mathematical concepts introduced in the previous sections, we have performed numerical calculations with a rather simple two-dimensional nonlinear model, that of Hénon–Heiles [26]. We deduce that high-order finite-difference methods based on the Lagrangian interpolation polynomials are appropriate for solving initial value problems, as well as partial differential equations, such as the Schrödinger equation, necessary in chemical theories. Finally, the conclusions are summarized in Section 5, where recent research on Hamiltonian neural networks (HNN) is discussed. Proofs for some equations, tables that summarize results in Hamiltonian thermodynamics, and an example of formulating a chemical kinetic model with thermodynamic Hamiltonian theory are presented in Appendix A.

2. Geometrical Structures in Chemical Dynamics

Molecules are sets of $[eqn]$ nuclei and $[eqn]$ electrons interacting with Coulomb forces. Usually, their quantum mechanical treatment is obtained by solving the Schrödinger equation in the Born–Oppenheimer approximation [27] that separates the electronic from the nuclear motion. The electronic energies are calculated by freezing the nuclei at specific configurations, which produces the adiabatic potential energy surface for each electronic state.

These molecular potentials, named also Potential Energy Surfaces, are employed to solve the nuclear equations of motion either in quantum mechanics or in classical mechanics. It is, thus, important to investigate common geometrical structures in the foundations of the two basic theories of physics, which in turn may assist in the numerical solutions of the corresponding equations of motion. In the following two subsections, we examine the topological and geometrical properties of classical and quantum mechanics, whereas in Section 3, the relatively new Hamiltonian formulation of thermodynamics is reviewed, all of them in extended phase spaces and at the non-relativistic approximation.

2.1. Canonical Classical Mechanics

2.1.1. Manifolds and Maps

We introduce the basic geometrical concepts of classical mechanics for time-dependent systems by elucidating the graphics of Figure 1. Starting from the bottom and moving upwards, we denote the set of configurations of the system with $[eqn]$ degrees of freedom (DOF) with the column vector (The symbol $[eqn]$ denotes the transpose of a matrix. Thus, a row vector becomes a column vector, and vice versa, a column vector is converted into a row vector.), $[eqn]$ , and the parameter (time) with $[eqn]$ . Capital letters designate the set of $[eqn]$ coordinates as a column vector, $[eqn]$ . These coordinates parametrize the Extended Configuration Manifold, $[eqn]$ . Generally, this is a smooth (differentiable) nonlinear manifold.

Manifolds generalize the geometrical objects of curves and surfaces in three-dimensional Euclidean space into $[eqn]$ dimensional ( $[eqn]$ ) spaces. It took more than two centuries to develop the current mathematical definition of a manifold since it required the parallel development of various other branches of Mathematics, such as topology, geometry, and algebra [28].

The extended configuration space $[eqn]$ can be described locally by a coordinate system (chart) Q, i.e., a homeomorphism

$$[eqn]$$

of an open set U of $[eqn]$ onto an open set $[eqn]$ of a Euclidean space $[eqn]$ of dimension $[eqn]$ .

In Euclidean space, we understand the definition of a coordinate system in $[eqn]$ as

$$[eqn]$$

for every point $[eqn]$ , and $[eqn]$ are the canonical projections taken to be differentiable functions. Transition maps provide the transformation from one coordinate system to another for points that belong to the intersection of two different open subsets.

The tangent space of $[eqn]$ at a point $[eqn]$ is a vector space, and the union of all tangent spaces for all points s of $[eqn]$ form the tangent bundle  $[eqn]$ with the extended configuration space $[eqn]$ to be the base space

$$[eqn]$$

The tangent bundle contains both the configuration manifold $[eqn]$ and its tangent spaces $[eqn]$ called the fibers, and it is a smooth manifold of dimension $[eqn]$ . Since $[eqn]$ is also a smooth manifold, a chart is defined by the diffeomorphism

$$[eqn]$$

Each coordinate system $[eqn]$ from the atlas of $[eqn]$ induces a coordinate system $[eqn]$ for $[eqn]$ . This chart is said to be the bundle chart associated with $[eqn]$ . The velocities $[eqn]$ live in this space.

The potential function, $[eqn]$ , is a function on the configuration manifold to real numbers. The Lagrangian, $[eqn]$ , is a function on the tangent bundle to real numbers.

The dual space of $[eqn]$ (the set of all linear maps on tangent bundle to real numbers) is the cotangent bundle ( $[eqn]$ ), also named phase space. The phase space is a differentiable manifold of $[eqn]$ dimension for which the tangent bundle can also be defined with charts described by the generalized coordinates ( $[eqn]$ ), the conjugate momenta (Notice that we use superscripts for coordinates and subscripts for momenta.) $[eqn]$ , and their velocities, where

$$[eqn]$$

The Extended Hamiltonian, $[eqn]$ , is a function on the phase space to real numbers, $[eqn]$ , obtained by a Legendre transformation ( $[eqn]$ ) of the Lagrangian. We may consider that the Legendre transformation generates a differentiable map between the tangent and cotangent bundles of $[eqn]$ , $[eqn]$ . $[eqn]$ and $[eqn]$ are canonical projections to extended configuration manifold of tangent and cotangent bundles, respectively. The tangent bundle of phase space is denoted by $[eqn]$ and $[eqn]$ is the canonical projection to phase space.

For a system of particles with $[eqn]$ DOF, we define the Lagrangian as the difference between kinetic energy $[eqn]$ and potential energy

$$[eqn]$$

with $[eqn]$ to be the kinetic metric tensor, written as a function of coordinates q and the particle masses m. The metric is the non-degenerate, symmetric, covariant tensor rank-2 that defines the kinetic energy. The momentum $[eqn]$ is the covector of the velocity $[eqn]$ ,

$$[eqn]$$

which is a map from the tangent bundle ( $[eqn]$ ) to the cotangent bundle ( $[eqn]$ ). Obviously, for a diagonal unit metric, $[eqn]$ , momenta are equal to velocities.

The velocities $[eqn]$ can also be extracted from the momenta by the inverse tensor $[eqn]$

$$[eqn]$$

where $[eqn]$ (The components of Kronecker delta tensor, $[eqn]$ , are equal to 1 for $[eqn]$ and 0 for $[eqn]$ .).

Momenta may also be considered to be $[eqn]$ forms which act on the tangent bundle to real numbers, $[eqn]$ , by taking the interior product (contraction) of $[eqn]$ forms with vector fields. Employing Dirac’s notation, we write

$$[eqn]$$

If $[eqn]$ is the Hamiltonian of a system of particles with $[eqn]$ DOF, then the extended Hamiltonian, $[eqn]$ , is defined with the Legendre transformation (Notice that in chemical thermodynamics the Legendre transformation is defined as the difference between the function and the sum of products of conjugate variables).

$$[eqn]$$

The physical states are obtained by imposing the two constraints

$$[eqn]$$

which result in

$$[eqn]$$

Hence, the time-extended system is a conservative system with Hamiltonian $[eqn]$ .

2.1.2. Equations of Motion

We collect the generalized coordinates $[eqn]$ and their conjugate momenta $[eqn]$ to a single column vector $[eqn]$ of $[eqn]$ dimension, where we use $[eqn]$ . Hamilton’s principle of stationary action leads to Hamilton’s equations of motion. Then, Hamilton’s equations with a Hamiltonian $[eqn]$ are written in the form

$$[eqn]$$

or

$$[eqn]$$

where $[eqn]$ is the gradient of Hamiltonian function, and J the symplectic matrix. J is the map on the tangent bundle of phase space M

$$[eqn]$$

$[eqn]$ and $[eqn]$ are the zero and unit $[eqn]$ matrices, respectively. It is proved that J satisfies the relations,

$$[eqn]$$

Thus, the Hamiltonian vector field is

$$[eqn]$$

or using the coordinate base in the tangent bundle of phase space $[eqn]$ , we write

$$[eqn]$$

that lives in the tangent space of phase space.

In a more general approach, we can extract the Hamiltonian vector field as follows. In the cotangent space of extended coordinate manifold, let us denote with $[eqn]$ the $[eqn]$ form acting on the phase space manifold $[eqn]$

$$[eqn]$$

and with $[eqn]$ the $[eqn]$ forms on the configuration manifold $[eqn]$

$$[eqn]$$

Since $[eqn]$ is a linear map from $[eqn]$ to M and $[eqn]$ an $[eqn]$ form on M we can pull-back $[eqn]$ to $[eqn]$ to produce the $[eqn]$ form $[eqn]$ , which lives on the base manifold $[eqn]$ . Then, the canonical Poincaré  $[eqn]$ form satisfies the relation (tautological  $[eqn]$ form)

$$[eqn]$$

Hence, we can expand $[eqn]$ as

$$[eqn]$$

$[eqn]$ is invariant under coordinate transformations

$$[eqn]$$

which we assume to be invertible

$$[eqn]$$

This is proved by arguing as follows. The velocities are given by

$$[eqn]$$

The new momenta are

$$[eqn]$$

Hence,

$$[eqn]$$

The canonical Symplectic  $[eqn]$ form is extracted by taking the exterior derivative (Notice the negative sign in our formulation) of $[eqn]$

$$[eqn]$$

This is a non-degenerate, skew-symmetric, closed  $[eqn]$ form ( $[eqn]$ ). In local coordinates $[eqn]$ , $[eqn]$ is expressed by the wedge products

$$[eqn]$$

If we introduce $[eqn]$ , the symplectic $[eqn]$ form (Equation (29)) is written

$$[eqn]$$

A pair $[eqn]$ , i.e., the phase space with the symplectic $[eqn]$ form, is said to be a symplectic manifold. Those local coordinates which satisfy, $[eqn]$ , are said to be canonical and symplectic. In the following, we shall see that Hamiltonian mechanics and its geometrical properties can be formulated by $[eqn]$ .

Let $[eqn]$ be a symplectic manifold of dimension $[eqn]$ with $[eqn]$ a canonical symplectic $[eqn]$ form. The Hamiltonian function $[eqn]$ is a smooth function on the phase space. The Hamiltonian vector field, $[eqn]$ , (Equations (17) and (18)) is then defined via the relationship

$$[eqn]$$

$[eqn]$ symbolizes interior product (contraction) and the triple $[eqn]$ is a Hamiltonian system. In particular, for the variable $[eqn]$ we extract the equation

$$[eqn]$$

We can see that with this formulation of time-dependent systems and taking into account the constraints Equations (11) and (12), the trajectories are described at each time t in the physical phase space of the system of $[eqn]$ dimension. With $[eqn]$ Hamilton’s equations of motion with the time-dependent Hamiltonian are written

$$[eqn]$$

For a smooth function $[eqn]$ on phase space, the Poisson bracket is defined as

$$[eqn]$$

Since

$$[eqn]$$ $$[eqn]$$

or employing the coordinate base in the tangent space of phase space

$$[eqn]$$ $$[eqn]$$

we write the Poisson brackets in a coordinate system as

$$[eqn]$$

We have used

$$[eqn]$$

where $[eqn]$ is the Kronecker’s delta tensor.

The Lie derivative of a dynamical quantity $[eqn]$ with respect to the Hamiltonian vector field $[eqn]$ is defined as the directional derivative of $[eqn]$ along the vector $[eqn]$

$$[eqn]$$

So,

$$[eqn]$$ $$[eqn]$$

From the above equation, we infer that for conserved quantities, $[eqn]$ , the Poisson’s bracket commutes, i.e., $[eqn]$ .

Applying the above equation to coordinates and conjugate momenta, $[eqn]$ , we extract Hamilton’s equations

$$[eqn]$$ $$[eqn]$$

Poisson brackets defined on a set of smooth functions $[eqn]$ on M satisfy the properties; For $[eqn]$ , then

$[eqn]$ is bilinear, $[eqn]$ antisymmetric, $[eqn]$ , and $[eqn]$ (Jacobi identity).

2.1.3. Integrable Hamiltonian Systems

In classical mechanics, a finite $[eqn]$ dimensional Hamiltonian system, $[eqn]$ is completely integrable if it admits n independent constants of motion whose Poisson brackets are in involution, i.e., pairwise commute.

Let the integrals of motion are $[eqn]$ and are assigned to values $[eqn]$ . Then, the corresponding level set, $[eqn]$ is a Lagrangian submanifold.

In a more formal wording, for an even-dimensional phase space, M, the symplectic $[eqn]$ form $[eqn]$ , which satisfies the condition (volume form)

$$[eqn]$$

the condition $[eqn]$ determines the $[eqn]$ Lagrangian submanifold $[eqn]$ .

For compact phase spaces, the Lagrangian submanifolds have the structure of a $[eqn]$ torus, $[eqn]$ . Moreover, in a neighborhood of every such invariant torus, one can find angle-action coordinates  $[eqn]$ ,

$$[eqn]$$

In such a coordinate system, the Hamiltonian function depends only on $[eqn]$ , i.e., $[eqn]$ , and Hamilton’s equations give

$$[eqn]$$

$[eqn]$ are the normal frequencies, and the action variables are constants. The Poisson brackets of the action coordinates, I, pairwise commute, $[eqn]$ , and also satisfy $[eqn]$ .

2.1.4. Complexification of Classical Hamilton’s Equations

Hamilton’s equations in classical mechanics can also be cast in a complex manifold by complexification of phase space, i.e., by introducing the complex transformation

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Similarly, in Section 2.2.3 by realification of the quantum Hilbert space, we bring the Schrödinger equation into the form of Hamilton’s equations.

To make the complex transformation canonical, we define the complex variables

$$[eqn]$$

The inverse transformation gives the real functions

$$[eqn]$$

Because of the symmetry of Hamiltonians under the unitary group $[eqn]$ , i.e., the Hamiltonian $[eqn]$ should be invariant with a unitary transformation. We write the canonical Poincaré $[eqn]$ form as

$$[eqn]$$ $$[eqn]$$

and thus, the canonical symplectic $[eqn]$ form becomes

$$[eqn]$$

For example, a harmonic oscillator in scaled normal coordinates and introducing these canonical complex coordinates results in

$$[eqn]$$ $$[eqn]$$

Since the transformation to complex coordinates and conjugate momenta is symplectic, Hamilton’s equations are also written as

$$[eqn]$$

where $[eqn]$ is the Hamiltonian in complex coordinates.

2.1.5. Jacobi Fields and Variational Equations

Geodesic curves are obtained by searching for the critical points of the length functional, i.e., by requiring

$$[eqn]$$

$[eqn]$ is the infinitesimal length of the curve given by the norm of the velocity vector.

The same critical points are obtained by varying the integrand $[eqn]$ instead

$$[eqn]$$

It is worth noting that the above integral depends on the parameter t. This integral is related to the action or the energy of a physical system. Indeed, Hamilton’s principle of stationary action results in Hamilton’s equations, the solutions of which are geodesics on the phase space manifold.

Let C be a trajectory in phase space with kinetic metric g

$$[eqn]$$

From the above equations, we deduce that the coefficients $[eqn]$ result from the action of the metric on the base coordinates of the tangent space

$$[eqn]$$

The length of a trajectory starting from the point $[eqn]$ and finishing at the point $[eqn]$ is calculated as the line integral

$$[eqn]$$

In case we want to investigate the behavior of neighboring trajectories $[eqn]$ to a reference one $[eqn]$ in phase space, we examine the time evolution of the variation vector $[eqn]$ , (Figure 2). The time derivative of the vector field $[eqn]$ with respect to the vector field $[eqn]$ (Lie derivative) is equal to

$$[eqn]$$

J is the symplectic matrix (Equation (15)). The higher order term (h.o.t.) is a function of the displacement $[eqn]$ at time t, which contains all the terms in the Taylor expansion larger than the first order. The derivatives are computed at the reference trajectory with initial conditions $[eqn]$ . $[eqn]$ is a Jacobi field and the equations

$$[eqn]$$

the variational equations.

The fundamental matrix  $[eqn]$ satisfies the variational equations with initial condition $[eqn]$ . It is a symplectic matrix, and therefore, the following equations are valid [3]

$$[eqn]$$

It is proved that if $[eqn]$ is an eigenvalue of a real symplectic matrix, so are $[eqn]$ and $[eqn]$ . $[eqn]$ are complex conjugate numbers. Also, for every constant of motion, two eigenvalues of the fundamental matrix are equal to one. For proofs and applications see references [4,5,29,30].

2.2. Geometrical Quantum Mechanics

2.2.1. Manifolds and Maps

Moving to the quantum world, the states of a system are described by complex vectors, $[eqn]$ , with s to be the configuration coordinates and t the time instead of the real coordinates and momenta (or velocities) in classical mechanics. The square of the norm of state vectors is interpreted as a probability density, with $[eqn]$ to be the probability of finding the system in the coordinate intervals $[eqn]$ at time t. The observable quantities, which are real-valued functions on phase space in classical mechanics, are replaced by operators that designate linear transformations on the complex state vector space, named Hilbert space of dimension n ( $[eqn]$ ) (Figure 3). Infinite-dimensional Hilbert spaces for quantum systems also exist.

The time evolution of the quantum states is given by Schrödinger’s equation

$$[eqn]$$

where $[eqn]$ , ℏ the reduced Heisenberg constant, $[eqn]$ is the Hamiltonian operator of the system that corresponds to its energy, and $[eqn]$ is the Hamiltonian–Schrödinger vector field. To any observable $[eqn]$ , we assign the operator $[eqn]$ and the Schrödinger vector field  $[eqn]$ . Hence, we may consider the vectors $[eqn]$ as vector fields, and thus, the Hilbert space $[eqn]$ to be also the tangent space $[eqn]$ at the state $[eqn]$ .

The cotangent space of the Hilbert space, $[eqn]$ , at state $[eqn]$ is $[eqn]$ , with covectors to be the complex conjugate functionals $[eqn]$ , and $[eqn]$ .

$[eqn]$ (We use Dirac’s (bra, ket) notation.) denotes the Hermitian inner product, ( $[eqn]$ ).

Separating the real and imaginary parts of the Hermitian inner product, we write

$$[eqn]$$

$[eqn]$ are the real and $[eqn]$ the imaginary components, respectively, of the kth component of the complex vectors. G is a Riemannian metric, i.e., it is real, positive definite, and strongly non-degenerate

$$[eqn]$$

$[eqn]$ is a symplectic (antisymmetric) closed $[eqn]$ form, and strongly non-degenerate

$$[eqn]$$

Both $[eqn]$ act on the tangent bundle $[eqn]$ .

Since the physical interpretation of a quantum state is probabilistic, i.e., for a pure state $[eqn]$ (the states of an isolated system), the following normalization condition should be satisfied

$$[eqn]$$

The observables, or measurable quantities of the system, $[eqn]$ , are represented by self-adjoint linear operators, $[eqn]$ , on $[eqn]$ , and are thus vector fields. The expectation value of an observable with operator $[eqn]$ at the state $[eqn]$ is the real-valued function

$$[eqn]$$

In particular, for the Hamiltonian operator of the system, we write

$$[eqn]$$

2.2.2. Projective Hilbert Space

For any nonzero factor $[eqn]$ , such as $[eqn]$ , $[eqn]$ yields the same expectation value for the observable $[eqn]$ as does $[eqn]$ . $[eqn]$ are the elements of the unitary Lie group, $[eqn]$ . Thus, the inclusion of these phases to the initial Hilbert space results in the extended Hilbert space of $[eqn]$ dimension, $[eqn]$ , which is isomorphic to $[eqn]$ with a Hermitian inner product.

The set of vectors obtained by multiplying the state $[eqn]$ with $[eqn]$ consists of a $[eqn]$ dimensional subspace of $[eqn]$ , called ray, which we symbolize as $[eqn]$ . A ray is an equivalence class of vectors in $[eqn]$ : two vectors are equivalent if and only if one is a nonzero complex scalar multiple of the other. Also, adopting normalized vectors (Equation (70)), the physical quantum states are elements of the complex Projective Hilbert space, $[eqn]$ of $[eqn]$ dimension obtained by the canonical projection map $[eqn]$ of the extended Hilbert space $[eqn]$ .

Hence, a quantum system may be described with what mathematicians call principal bundle

$$[eqn]$$

with state vectors the rays $[eqn]$ in $[eqn]$ , and physical states $[eqn]$ in $[eqn]$ obtained by the projection map $[eqn]$ . The inverse projection is

$$[eqn]$$

To recapitulate, the inclusion of the unitary group $[eqn]$ extends the $[eqn]$ dimensional Hilbert space $[eqn]$ to the $[eqn]$ extended Hilbert space. For normalized state vectors, the extended Hilbert space is mapped to the unit sphere $[eqn]$ of dimension $[eqn]$ embedded in the real Euclidean $[eqn]$ space. Finally, by the projection $[eqn]$ we produce the Projective Hilbert space isomorphic to the complex space $[eqn]$ or to real of $[eqn]$ dimensional space, $[eqn]$ ,

$$[eqn]$$

The tangent space of the projective space at point $[eqn]$ , $[eqn]$ , is isomorphic to the kernel of the ray, $[eqn]$ , i.e.,

$$[eqn]$$

and the push-forward of the projection map, $[eqn]$ , is a complex conjugate linear isomorphism onto $[eqn]$

$$[eqn]$$

An illustrative figure is shown in the article of Dorje C. Brody and Lane P. Hughston [7].

For a normalized ray and its kernel, one can argue that [6,7,8]

$$[eqn]$$

which gives a well-defined Hermitian inner product on $[eqn]$ . Moreover, we deduce that the equation

$$[eqn]$$

provides a strong symplectic form on $[eqn]$ , and the equation

$$[eqn]$$

defines a strong Riemannian metric on $[eqn]$ called the Fubini–Study metric. ℜ and ℑ are the real and imaginary parts of the Hermitian inner product, respectively, in the extended complex Hilbert space $[eqn]$ . Both $[eqn]$ and $[eqn]$ are invariant under all transformations $[eqn]$ , for all unitary operators $[eqn]$ on $[eqn]$ .

2.2.3. Realification of Hilbert Space and Kähler Manifolds

Expanding the vectors of Hilbert space in a $[eqn]$ dimensional basis, as well as its dual basis, we write the $[eqn]$ coefficients $[eqn]$ with the real $[eqn]$ and imaginary $[eqn]$ components. Similarly, the covectors are written as $[eqn]$ for the real and $[eqn]$ for the imaginary components. Thus, working with the real vectors $[eqn]$ we can describe quantum states in a real vector space $[eqn]$ .

An almost complex structure is introduced in $[eqn]$ by replacing the imaginary number $[eqn]$ with the symplectic matrix $[eqn]$ , Equation (15). We also derive

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The triple $[eqn]$ attributes a Kähler structure, and thus, $[eqn]$ is a Kähler manifold [2].

By realification of the Hilbert space, we may express the Hamiltonian–Schrödinger vector field as

$$[eqn]$$

Also, for an observable, $[eqn]$ , the corresponding Schrödinger vector field is

$$[eqn]$$

For a time-independent Hamiltonian, $[eqn]$ , we solve the Schrödinger equation by employing the unitary propagator $[eqn]$ , and write the Hamiltonian flow as

$$[eqn]$$

$[eqn]$ generates a one-parameter group of transformations on $[eqn]$ , which preserve the metric G and the symplectic $[eqn]$ form $[eqn]$ .

The triple $[eqn]$ is a Hamiltonian system with Hamiltonian function the Expectation value of $[eqn]$ at the state $[eqn]$

$$[eqn]$$

We can prove (see Appendix A.1) that for an observable $[eqn]$ and normalized states $[eqn]$ , the differential $[eqn]$ form of the expectation value of $[eqn]$ is the interior product of the symplectic $[eqn]$ form $[eqn]$ with the Schrödinger vector field $[eqn]$ (Equation (86))

$$[eqn]$$

where $[eqn]$ . The above equation is in accordance with Equation (31) of classical mechanics. Hence, Hermitian operators $[eqn]$ give rise to quadratic real-valued functions $[eqn]$ . The vector field $[eqn]$ generated by the expectation value function of the operator $[eqn]$ is a Schrödinger (Hamiltonian) field.

If $[eqn]$ are the expectation value functions of two observables, the Poisson bracket (Equation (39)) is defined by the equation (see Appendix A.2)

$$[eqn]$$

$[eqn]$ is the commutator of the two operators $[eqn]$ .

The uncertainty (dispersion) of a quantum observable $[eqn]$ in a normalized state $[eqn]$ is defined as [8,31]

$$[eqn]$$

If $[eqn]$ is an eigenvalue of the observable $[eqn]$ , then, $[eqn]$ .

For two quantum observables $[eqn]$ the product of the uncertainties of the expectation value functions $[eqn]$ is written as

$$[eqn]$$

and the covariance or correlation function is expressed as

$$[eqn]$$

where the anticommutator is denoted by $[eqn]$ . It is proved that

$$[eqn]$$

The Schrödinger vector fields are expressed as usually

$$[eqn]$$

If we use Equation (90), then we can express the commutator

$$[eqn]$$

Schwartz inequality implies

$$[eqn]$$

Since

$$[eqn]$$

and

$$[eqn]$$

the Robertson-Schrödinger uncertainty relation is written as (Notice that to be consistent with the literature in discussing the covariance of two observables, we introduce the factor $[eqn]$ in the commutator and anticommutator).

$$[eqn]$$ $$[eqn]$$

2.2.4. Equations of Motion

We have seen that in the Extended Hilbert space, the vector states are the rays $[eqn]$ , which for simplicity we symbolize as $[eqn]$ . The projection of the Extended Hilbert space on the $[eqn]$ dimensional tangent Projective Hilbert phase makes the solutions of the time-dependent Schrödinger equation to be equivalent to the solutions of Hamilton’s equations in the Projective Hilbert space, with the expectation value of the Hamiltonian operator playing the role of Hamiltonian function. Furthermore, the phase space is a Kähler manifold, both in the extended Hilbert space as well as in the Projective Hilbert space.

If we use as basis set the coordinate basis $[eqn]$ the representation of the state $[eqn]$ is

$$[eqn]$$

The coefficient $[eqn]$ is called wavefunction and it is a complex function, $[eqn]$ , which we assume normalized to one

$$[eqn]$$

With an orthonormal coordinate basis, the completeness relation is written as

$$[eqn]$$

In the following, we adopt the coordinate representation of the state vectors using wavefunctions in the Schrödinger picture, $[eqn]$ . We expand $[eqn]$ of a dynamical system in an arbitrary orthonormal basis set, $[eqn]$

$$[eqn]$$

In this expansion, the basis functions $[eqn]$ are time-independent, whereas the coefficients $[eqn]$ depend on time. $[eqn]$ are solutions of the Schrödinger equation evolving in the extended Hilbert space

$$[eqn]$$

and their complex conjugate solutions of the equation

$$[eqn]$$

We assume $[eqn]$ to be normalized, $[eqn]$ , at any time. Then, by substituting Equation (104) in the expectation value of the Hamiltonian, $[eqn]$ , we obtain

$$[eqn]$$

The Schrödinger equation in the extended Hilbert space is mapped into Hamilton’s equations in the Projective Hilbert space $[eqn]$ with Hamiltonian function the expectation value of the Hamiltonian operator $[eqn]$ . For simplicity, we again denote the quantum states $[eqn]$ in the Projective Hilbert space with $[eqn]$ . Then, we differentiate Equation (107) with respect to $[eqn]$

$$[eqn]$$

Similarly, we take

$$[eqn]$$

We define the complex variables ( $[eqn]$ ) by introducing the real functions $[eqn]$ and $[eqn]$ to correspond to real and imaginary parts of the complex variables, respectively,

$$[eqn]$$

Equations (108) and (109) are the quantum equivalent of Hamilton’s equations of motion

$$[eqn]$$ $$[eqn]$$

Hence,

$$[eqn]$$

$[eqn]$ are the transition frequencies. It is also worth noting that the quantum Hamiltonian is similar to that of a harmonic oscillator after complexification (Equation (56)).

We can take the inverse of Equations (110) and write the real functions $[eqn]$ and $[eqn]$ as

$$[eqn]$$

Adopting the realification of the complex Projective Hilbert space, we may transform Hamilton’s equations to

$$[eqn]$$ $$[eqn]$$

or

$$[eqn]$$

$[eqn]$ is the representation of Hamiltonian operator in the basis set $[eqn]$ , and thus, a Hermitian matrix (The Hamiltonian operator is Hermitian, and its representation in a basis results in a Hermitian matrix, $[eqn]$ ).

The Hamiltonian vector field is extracted from an equation similar to that of classical mechanics (Equation (31)). Hence, in the realified phase space with coordinates $[eqn]$ the canonical Poincaré  $[eqn]$ form

$$[eqn]$$

provides the symplectic $[eqn]$ form

$$[eqn]$$

Then, the interior product (contraction) of $[eqn]$ gives

$$[eqn]$$

and the Hamiltonian vector field is extracted as

$$[eqn]$$

We define the length of a curve $[eqn]$ in the Projective Hilbert space by introducing the Riemannian metric $[eqn]$ as

$$[eqn]$$

For real wavefunction, the length is the Euclidean distance

$$[eqn]$$

2.2.5. Quantum Systems as Totally Integrable Hamiltonian Systems

It has been shown that geometrical quantum mechanics naturally describes a completely integrable system, as follows. Let $[eqn]$ be a complex separable Hilbert space of dimension $[eqn]$ , and view the triple $[eqn]$ as a Hamiltonian dynamical system on the phase space $[eqn]$ , equipped with the symplectic $[eqn]$ form $[eqn]$ . Let $[eqn]$ be the self-adjoint Hamiltonian operator for the system, and assume that each eigenspace of $[eqn]$ is one-dimensional (non-degenerate). Choose an orthonormal basis $[eqn]$ for $[eqn]$ consisting of eigenvectors of $[eqn]$ . Define the projection operators $[eqn]$ by

$$[eqn]$$

where we have used Equation (104).

Without loss of generality, we set the lowest eigenvalue ( $[eqn]$ ) to be 0 so that we can expand the Hamiltonian operator as

$$[eqn]$$

Observe that the projectors $[eqn]$ form a mutually commuting set of n operators on $[eqn]$ , and we can define the corresponding expectation value functions $[eqn]$ in the Projective Hilbert space $[eqn]$ as

$$[eqn]$$

Hence,

$$[eqn]$$

$[eqn]$ are the real and imaginary parts of the wavefunction. Thus, the Hamiltonian in the projective Hilbert space resembles that of a sum of n harmonic oscillators in scaled canonical coordinates (Equation (55)).

Therefore, we may conclude that there are n constants of motion, $[eqn]$ , and the trajectories lie on $[eqn]$ dimensional tori (Lagrangian tori).

3. Hamiltonian Chemical Thermodynamics

3.1. Manifolds and Maps

To establish an analogous theoretical framework for classical thermodynamics to that of Hamiltonian classical mechanics, we consider the macroscopic properties of the system to consist of the configuration manifold (Figure 4). For a typical chemical system we usually take the extensive physical properties of entropy (S), internal energy (U), volume (V), and the number of molecules (or moles) ( $[eqn]$ ) of M chemical species to consist of the coordinates of the extended configuration manifold $[eqn]$ , where $[eqn]$ . If we attribute entropy to be a homogeneous first-degree function of $[eqn]$ , then according to Euler’s theorem for homogeneous functions we can write

$$[eqn]$$

where we have introduced the conjugate intensive variables $[eqn]$ as the partial derivatives of entropy

$$[eqn]$$

T is the absolute temperature, P the pressure, and $[eqn]$ the chemical potentials of M compounds.

$[eqn]$ and entropy comprise of a $[eqn]$ coordinate system for the contact manifold $[eqn]$ , where the physical states of the system live. The extended configuration manifold and the contact manifold are generally nonlinear, and the maps $[eqn]$ and $[eqn]$ determine local coordinate systems in Euclidean space.

The total differential of S is Gibbs’s fundamental equation that describes the physical thermodynamic submanifold (PTS) of a thermodynamic system

$$[eqn]$$

If we assign the entropy S to coordinate $[eqn]$ , we can also write

$$[eqn]$$

From Equations (129) and (130) we deduce that the $[eqn]$ form $[eqn]$ satisfies

$$[eqn]$$

This equation expresses the physical state submanifold of a thermodynamic system, named Legendrian submanifold, and it is what we represent with the equation of states [32].

In a more formal wording, in the contact space $[eqn]$ the locally defined $[eqn]$ form $[eqn]$ that satisfies the condition, i.e., there is the volume form

$$[eqn]$$

ascertains the $[eqn]$ Legendrian submanifold $[eqn]$ of $[eqn]$ by requiring $[eqn]$ . In other words, the kernel of $[eqn]$ provides the maximal dimension hyperplanes tangent to $[eqn]$ .

To construct a Hamiltonian system, we need to define the canonical conjugate momenta. Following Balian and Valentine [18], we denote the conjugate momentum of entropy with $[eqn]$ , and we take it as a non-vanishing free parameter or a gauge variable. Then, the canonical conjugate momenta to the other coordinates of $[eqn]$ are defined as

$$[eqn]$$

Hence, $[eqn]$ act as new variables, which replace the physical intensive variables $[eqn]$ . We write $[eqn]$ and Equation (131) determines the Poincaré $[eqn]$ form

$$[eqn]$$

which produces the symplectic $[eqn]$ form

$$[eqn]$$

This equation sets the Thermodynamic Extended Physical State Submanifold (TEPSS) in the extended phase space, $[eqn]$ , named Lagrangian submanifold.

In a more formal wording, for an even-dimensional phase space, $[eqn]$ , the symplectic $[eqn]$ form $[eqn]$ (see also Equation (29)), which satisfies the condition (volume form)

$$[eqn]$$

determines the $[eqn]$ Lagrangian submanifold $[eqn]$ in $[eqn]$ by requiring $[eqn]$ .

The thermodynamic contact space $[eqn]$ is the projective space of thermodynamic extended phase space $[eqn]$ (see Appendix A.4.1)

$$[eqn]$$

where $[eqn]$ is the projection map.

It is proved that the submanifold $[eqn]$ is a Legendrian submanifold, if and only if

$$[eqn]$$

is a homogeneous in momenta Lagrangian submanifold (see Table A1). This means that $[eqn]$ as well as

$$[eqn]$$

We say that all $[eqn]$ , which satisfy Equation (136), belong to the ray $[eqn]$ . Then, it is valid $[eqn]$ for every vector field tangent to $[eqn]$ . Also, every homogeneous Lagrangian submanifold originates from a Legendrian submanifold of the form $[eqn]$ . We can also state that a submanifold is homogeneous if its generating function is homogeneous.

In the entropy representation, the n extensive independent variables $[eqn]$ together with the $[eqn]$ parameterize the Lagrangian submanifold. If the entropy $[eqn]$ is the generating function of the Legendrian submanifold, then for the Lagrangian submanifold the generating function is $[eqn]$ . This is a property of homogeneous functions of first-degree. Similarly, in energy representation we take $[eqn]$ as generating function for the Lagrangian submanifold [29]. With the generating function in the entropy representation, we extract the remaining $[eqn]$ variables as

$$[eqn]$$

The extended Hamiltonian $[eqn]$ is a function on phase space to real numbers, whereas the tangent bundle of phase space, $[eqn]$ , is described by coordinates and momenta, $[eqn]$ , as well as their time derivatives. The extended Hamiltonian is a homogeneous function of first-degree in momenta, i.e.,

$$[eqn]$$

$[eqn]$ is the projection map of the tangent bundle of phase space to phase space.

In Appendix A.4, we summarize in Table A2 and Table A3 several formulae of Hamiltonian thermodynamics in entropy- and energy-representations, respectively.

3.2. Equations of Motion

The Hamiltonian function $[eqn]$ acting in the extended phase space is a homogeneous function of first-degree in momenta [18] and according to Euler’s theorem for homogeneous functions of first-degree, we can write

$$[eqn]$$

For a simple thermodynamic system, the extended Hamiltonian is $[eqn]$ . To extract the physical states, we impose the constrain $[eqn]$ , which implies $[eqn]$ on thermodynamic extended physical state submanifold, also written as [29]

$$[eqn]$$

We point out that a Hamiltonian system is the triple $[eqn]$ , where $[eqn]$ is the Hamiltonian vector field defined on the tangent bundle of phase space, and $[eqn]$ a skew-symmetric, non-degenerate, closed differential $[eqn]$ form. Then, we can use the machinery of geometrical classical mechanics as was developed in Section 2.1 to write down equations for equilibrium and non-equilibrium processes. Thus, the Hamiltonian vector field can be extracted from the equation

$$[eqn]$$

given the Hamiltonian function $[eqn]$ on the phase space $[eqn]$ . This equation provides the Hamiltonian vector field $[eqn]$

$$[eqn]$$

From the definition of the extended Poincaré $[eqn]$ form (Equation (133)) we prove that

$$[eqn]$$

for homogeneous in momenta of first-degree Hamiltonian functions.

Hence, Hamilton’s equations take their usual form (Equation (33))

$$[eqn]$$

Moreover, $[eqn]$ is a canonical coordinate system in the corresponding phase space, and thus, satisfy the Poisson brackets Equation (34)

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

3.2.1. Contact Equations of Motion

Since the Poisson bracket (Equation (34)) is invariant under gauge (canonical) transformations, we obtain

$$[eqn]$$

where F is a generating function and

$$[eqn]$$

we extract the contact vector field (see Appendix A.3)

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

For a generating function F, which is independent of $[eqn]$ variable, as well as a homogeneous function of first-degree in momenta, the contact equations of motion are reduced to those of Hamilton’s equations in $[eqn]$ phase space, i.e.,

$$[eqn]$$

3.2.2. Riemannian Metric on Lagrangian Submanifold

We can define a Riemannian metric [28] on a Lagrangian $[eqn]$ submanifold, $[eqn]$ , with generating function $[eqn]$ . In canonical coordinates $[eqn]$ it takes the form (see Table A2)

$$[eqn]$$

$[eqn]$ is the metric introduced by Ruppeiner [33,34] and it is usually called Ruppeiner metric.

Weinhold [35,36] has extracted a metric, $[eqn]$ , in the energy representation of thermodynamics with generating function $[eqn]$ [29] (see Table A3).

3.3. Embedding Systems in Homogeneous Media

We examine the case of a simple system with internal energy U and constant volume and number of particles. The embedded mechanical system is described by the Hamiltonian $[eqn]$ . Here, we consider the augmented system, system plus environment, isolated with total energy E constant, $[eqn]$

$$[eqn]$$

The entropy of the thermodynamic system is given by the equation

$$[eqn]$$

i.e., a function of $[eqn]$ independent variables. In the extended thermodynamic phase space, we include entropy as an independent variable as well

$$[eqn]$$

The canonical conjugate momenta are $[eqn]$ and $[eqn]$ . Hence, the dimension of thermodynamic extended phase space is $[eqn]$ .

In the entropy representation, we take as the generating function of the Legendrian submanifold in the thermodynamic contact space the entropy $[eqn]$ . Hence, the generating function of the corresponding Lagrangian submanifold in the extended phase space is $[eqn]$ . For $[eqn]$ the canonical momenta become

$$[eqn]$$

Therefore, we have attained the important conclusion that the thermodynamic extended state manifold in thermodynamic extended phase space is the Lagrangian submanifold, $[eqn]$ ,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

From Equation (140) the Hamiltonian function in the extended phase space is written as

$$[eqn]$$

Hamilton’s equations are then determined

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where the velocities $[eqn]$ are determined from the definition of the embedded system.

3.4. Chemical Kinetics

The kinetics of chemical reaction networks can also be formulated using a geometric Hamiltonian theory as an embedded system. However, classical chemical kinetics is better studied by introducing temperature, pressure, and the elementary chemical reaction coordinates $[eqn]$ as the independent variables.

3.4.1. Thermodynamics of Chemical Reactions

A general chemical reaction with $[eqn]$ reactant compounds and $[eqn]$ product compounds can be expressed by the formula

$$[eqn]$$

where $[eqn]$ . The stoichiometric integers $[eqn]$ are negative for the reactant molecules $[eqn]$ and positive $[eqn]$ , for the product molecules, $[eqn]$ .

For K elementary chemical reactions, we must employ K independent reaction coordinates $[eqn]$ . Then, the time variation of the number of molecules is related to the velocities of reactions by the equation

$$[eqn]$$

or simply as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ $[eqn]$ , being the stoichiometric matrix for the K elementary reactions.

The reaction Gibbs free energy $[eqn]$ is defined by introducing the chemical potentials of the $[eqn]$ compounds

$$[eqn]$$

where $[eqn]$ is Gibbs free energy, T the temperature and P the pressure. From Gibbs’s fundamental equation

$$[eqn]$$

and for a general chemical reaction at constant temperature and pressure, Equation (169), we have

$$[eqn]$$

Therefore, the above equation becomes

$$[eqn]$$

Usually, the derivative of Gibbs free energy with the reaction coordinate is denoted as $[eqn]$ and is called molar reaction Gibbs free energy.

For K elementary reactions, we write

$$[eqn]$$

employing the rates of forward ( $[eqn]$ ) and backward ( $[eqn]$ ) reactions, as expressed by the law of Mass Action [37,38,39,40].

Also, for an elementary chemical reaction De Donder [41,42] introduced the quantity of Affinity

$$[eqn]$$

which can be written as [30]

$$[eqn]$$

3.4.2. Thermodynamic Hamiltonian in Massieu-Gibbs Representation

Chemical reactions are better studied at constant temperature and pressure. Hence, we first transform the entropy from a function of the internal energy and volume to a function in temperature and pressure by a Legendre transformation

$$[eqn]$$ $$[eqn]$$

$[eqn]$ is Massieu function related to Gibbs free energy $[eqn]$ [19].

The above Legendre transformation is a diffeomorphism on the Legendrian and Lagrangian submanifolds. This means that we can define Affinity with the chemical potentials $[eqn]$ of M chemical compounds (reactants and products) given by Gibbs free energy

$$[eqn]$$

The rate of entropy production for the kth-irreversible reaction is written with the molar reaction Gibbs free energy $[eqn]$

$$[eqn]$$

where $[eqn]$ is the Affinity of kth-reaction and $[eqn]$ is the reaction coordinate that describes the progress (extent) of the reaction. The total entropy rate for K reactions becomes

$$[eqn]$$

or employing Equations (176) and (178), as

$$[eqn]$$

where R is the ideal gas constant. Obviously, for reactions at equilibrium in a closed thermodynamic system, we have $[eqn]$ , which means that $[eqn]$ , and thus, the total entropy is conserved. The inequality in Equation (183), which is consistent with the second thermodynamic law, holds because ln is a monotonous function in the domain of its definition.

The thermodynamic extended state manifold in thermodynamic extended phase space is a Lagrangian submanifold, $[eqn]$ embedding in $[eqn]$ phase space with $[eqn]$ . The generating function is $[eqn]$

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

$[eqn]$ is the conjugate momentum of $[eqn]$ .

The Hamiltonian function in the extended cotangent bundle (phase space) is written as

$$[eqn]$$

At constant temperature and pressure and using concentrations instead of the number of molecules for reactants and products, $[eqn]$ with conjugate momenta also denoted by $[eqn]$ , we take

$$[eqn]$$

For K elementary reactions in the reaction coordinate space, $[eqn]$ , and under the transformation of stoichiometric matrix $[eqn]$ we have

$$[eqn]$$

and the Hamiltonian (Equation (189)) becomes

$$[eqn]$$

Then, Hamilton’s equations are written

$$[eqn]$$

We can also define a metric on the Lagrangian submanifold in Massieu-Gibbs representation, and this is given as

$$[eqn]$$

Therefore, to measure the distance between initial and equilibrium states, we compute the length of a path on the thermodynamic extended state manifold

$$[eqn]$$

with $[eqn]$ to be the maximum integration time and substituting $[eqn]$ . The distance of two states has been utilized to define a better low bound of the entropy production or dissipated work (availability) than zero for finite time irreversible processes [43,44].

In Appendix A.5, an example of two elementary consecutive chemical reactions is presented and described in detail.

4. Numerical Implementations

4.1. High-Order Finite-Difference and Pseudospectral Methods

High-order finite-difference methods are widespread for solving the differential equations encountered in classical and quantum mechanics. They are based on polynomial approximations to the solution functions, and thus, finite-difference methods (FD) are suited for problems that can be solved by expanding the unknown functions to Taylor series [25,45]. In the past, we developed and tested high-order multivariable finite-difference methods in configuration space and time (herein they are denoted with x), formulated by Lagrange interpolating polynomials [46,47]. Given that variable order finite-difference algorithms are among the most suitable for modern high-performance computing, the computer technology to which computational chemistry is mainly addressed, as well as their programming simplicity, finite-difference methods emerge as one of the best choices for studying chemical dynamics.

There is a huge literature concerning finite-difference approximations to initial value problems (Notice that Hamilton’s equations, Jacobi fields, and Poisson brackets can be written as initial value problems, and thus, ordinary differential equations can simultaneously be integrated.), such as Hamilton’s and variational equations [45], as well as calculating the action of Hamiltonian operator on wavefunctions in quantum mechanics [25].

For ordinary differential equations (ODEs) and partial differential equations, we approximate the solution functions by expanding them in an appropriate basis set $[eqn]$

$$[eqn]$$

Different global basis functions produce different pseudospectral methods. From such so-called finite basis representation, we can transform to a cardinal set of basis functions, $[eqn]$ , also called discrete variable representation, by choosing N grid points, $[eqn]$ , at which the function is calculated. The cardinal functions obey the $[eqn]$ Kronecker property

$$[eqn]$$

so that the wavefunction is expressed by the set of grid points

$$[eqn]$$

The transformation from $[eqn]$ to the cardinal basis set is unitary, and the new basis is given in terms of the old one by the equation

$$[eqn]$$

where the grid points $[eqn]$ and the corresponding weights $[eqn]$ depend on the chosen quadrature rule.

The mth-derivative of the approximate solution is written as

$$[eqn]$$

where

$$[eqn]$$

The differentiation matrix  $[eqn]$ contains the coefficients necessary for calculating the mth- derivative at the collocation points, and T is the column vector of dimension N containing the basis functions.

Usually, the functions are interpolated by Lagrange cardinal basis set of order $[eqn]$

$$[eqn]$$

where

$$[eqn]$$

It is proved that by increasing the order of Lagrange polynomials, the series converges to the corresponding pseudospectral limits. For uniform equidistant grids, $[eqn]$ , and $[eqn]$ , it is shown that the pseudospectral limit is a $[eqn]$ function

$$[eqn]$$

Also, for periodic functions, Fourier series are associated with the cardinal functions

$$[eqn]$$

A relation of Fourier cardinal functions and sinc functions can be seen by the formula

$$[eqn]$$

Since sinc functions are the infinite order limit of an equispaced FD (Equation (203)), the correspondence now is that periodically repeated FD stencils will tend to the PS limit of Fourier functions as the number of grid points in the stencil approaches the total number of grid points in one period. Also, because equispaced FD can be considered to be a robust sum acceleration scheme of a sinc function series, we expect the same convergence properties of the FD approximation to the Fourier series as the one we find for radial variables [46].

To solve Hamilton’s and variational differential equations [5], we use methods for solving ordinary differential equations, such as those described in the book of Shampine and Gordon [45]. The algorithms are based on predictor-corrector methods with variable order finite-difference approximations, as well as variable time step and backward difference formulae. The accuracy of solutions is controlled by prespecified values.

The action of the Hamiltonian operator on a wavefunction can also be approximated by high-order finite-difference schemes or pseudospectral methods. We have demonstrated that the algorithm developed by Fornberg [25] to construct Lagrange interpolating polynomials for solving the Schrödinger equation is robust and fast. Applications can be found in references [46,47].

4.2. The Hénon–Heiles Model

In this section, we show results that demonstrate the effectiveness of solving the Hamiltonian ODEs encountered in the three principal theories: classical, quantum, and thermodynamic. We use a test model that of Hénon–Heiles [26] and mainly the software POMULT [48]. The Hénon–Heiles potential was initially proposed as a model for studying the dynamics of a galaxy. The potential function introduces cubic nonlinearities, and soon, it became the model of choice for numerically investigating the nonlinear behavior of two degrees of freedom systems.

Numerous articles have been published that explore the phase space structures of this nonlinear dynamical system in detail [49]. The Hamiltonian function for a Hénon–Heiles system is written as

$$[eqn]$$

$[eqn]$ denote the coordinates and $[eqn]$ the canonical conjugate momenta. In Figure 5, we present 3D graphs and isocontours from the projections of the potential in the coordinate plane. The potential has one minimum at $[eqn]$ and three saddles at the energy of $[eqn]$ . Trajectories with energy above D may escape to infinity from the three exit channels. The harmonic normal mode frequencies are $[eqn]$ , and thus, the system shows a 1:1 resonance. Notice that the spatial symmetry of the Hénon–Heiles potential is $[eqn]$ , the same symmetry as triatomic molecules with the same atom.

4.2.1. A Classical Time-Dependent Hénon–Heiles System

A time-dependent Hénon–Heiles system is produced by adding the $[eqn]$ term coupled to $[eqn]$ coordinate (see Figure 1)

$$[eqn]$$

with $[eqn]$ , and canonical conjugate momentum, $[eqn]$ . We solve Hamilton’s equations as given by Equations (14), and with the extended Hamiltonian, $[eqn]$ .

In Figure 6, we depict a dissociating trajectory with parameters of the time-dependent term, $[eqn]$ and $[eqn]$ . The initial conditions are selected from the stable periodic orbit of the principal family along the $[eqn]$ coordinate and at energy 0.1575. The red thick line in Figure 6b is the projection of this periodic orbit on the coordinate plane. The 3D plot in Figure 6a shows the evolution of the system with the absorbed energy that leads to dissociation.

4.2.2. The Quantum Hénon–Heiles System

We solve the time-dependent Schrödinger equation with an initial Gaussian wavepacket centered at the minimum of the Hénon–Heiles potential, and widths $[eqn]$ and $[eqn]$ . The average energy is 0.124.

The autocorrelation function is the overlap integral of the initial Gaussian with the evolving wavepacket

$$[eqn]$$

The Fourier transformation of the autocorrelation function depicts the power spectrum from which one can extract eigenenergies and eigenfunctions. Numerical technologies of accurately calculating the energy levels, as well as the corresponding eigenfunctions, have extensively been investigated [50].

The outputs of a typical run are depicted in Figure 7 and Figure 8.

4.2.3. Energy Dissipation of Hénon–Heiles System in Homogeneous Media

For a simple thermodynamic system, such as that of an inert atomic gas with constant volume and particle number (1 mole), the entropy function is given by the equation

$$[eqn]$$

with reference the energy $[eqn]$ at $[eqn]$ and $[eqn]$ to be the specific heat. $[eqn]$ denotes the Boltzmann constant and $[eqn]$ the Avogadro constant.

We assume a diagonal friction parameter matrix $[eqn]$ that couples the Hénon–Heiles nonlinear oscillator with the homogeneous environment. We assign the nonzero elements $[eqn]$ to the values, $[eqn]$ . This is an example of loosely coupling one DOF of the dynamical system to the environment, and we investigate the effect it may have on the entropy production and length of the path [29].

For zero frictions, the dynamical system remains uncoupled to the environment and thus, is a conserved Hamiltonian system. When the friction parameters are turned on, energy dissipates from the dynamical system to the environment. In the extended thermodynamic space, the number of DOF is 6 with coordinates the entropy, S, the total energy, E, and the four variables of the dynamical system, $[eqn]$ . Hence, the dimension of the extended phase space is 12.

The extended Hamiltonian $[eqn]$ (Equation (162)) is equal to zero on the Lagrangian thermodynamic extended state manifold. We integrate trajectories up to 50 time units with a precision in the value of extended Hamiltonian to about $[eqn]$ .

In Figure 9 and Figure 10, the results for two trajectories are depicted by treating the Hénon–Heiles as a dissipating system. We plot the produced entropy (Equation (163)) and the lengths of trajectories, computed via the Ruppeiner metric (Rlength), Equation (154), and with initial energy of the dynamical system $[eqn]$ .

Since one DOF of the dynamical system is weakly coupled to the environment, we find that for the integration times, trajectories are trapped around reduced tori related to the approximately decoupled mode. The thermodynamic length provides an estimate of how close the initial state is to the physical equilibrium state. Hence, we expect trapped trajectories in the approximately decoupled modes to have large lengths and to produce low entropies [29].

5. Conclusions

All cardinal physical theories have acquired a Lagrangian or Hamiltonian formalism. The advantages of employing a Hamiltonian framework stem from the fact that global and local constants of motion that dictate the dynamics of the system are respected by solving Hamilton’s equations. In particular, the geometrical structures of Hamiltonian theory were mainly investigated in the second half of the twentieth century with the development of contemporary differential geometry. The purpose of this article is to introduce and unveil common geometrical structures of the three most frequently applied theories in computational chemistry, those of classical mechanics, quantum mechanics, and classical thermodynamics, all of them at the non-relativistic approximation. We have shown that working in extended phase space, the physical states of the system in the three theories are described by Lagrangian submanifolds. Observables are calculated by canonically projecting the extended phase space in a reduced dimensional space. In classical mechanics, integrable systems guarantee the existence of $[eqn]$ constants of motion, and thus $[eqn]$ dimensional Lagrangian submanifolds embedded into $[eqn]$ dimensional phase space. Quantum systems can also be considered to be integrable systems in the projective Hilbert space. Finally, classical thermodynamics have also been given a Hamiltonian formalism in an extended phase space, similar to classical mechanics, and capable of formulating irreversible processes. We can obtain representations of Lagrangian manifolds when we use Hamiltonians in action-angle variables in classical mechanics and Gibbs’s fundamental equation in classical thermodynamics.

Noether’s theorem [51] proves that constants of motion are associated with symmetries of the system, which leave the Hamiltonian function invariant. Global continuous symmetries, such as time-reversal, translational, and rotational transformations, are theorized by Lie groups and their associated Lie algebras. Advanced methods to extract the constants of motion related to symmetries and to reduce the dimensionality of the problem have been developed in recent decades [2].

Here, we have not explored the important role of symmetry in dynamical systems. However, we have emphasized the importance of approximating local constants of motion by solving the variational equations of a Hamiltonian system, such as around equilibria and stable periodic orbits. This requires the diagonalization of the fundamental matrix (see Section 2.1.5). It is proved that each constant of motion corresponds to a pair of eigenvalues equal to one.

Numerically solving the equations of motion (classical, quantum, thermodynamic) is the main task of computational chemistry. Finite-difference methods have been broadly adopted for finding solutions of ordinary differential equations, as well as partial differential equations, such as the Schrödinger equation. Nevertheless, in spite of their many successes, they are restricted to low-dimensional systems since they require substantial computational resources for many degrees of freedom systems. Furthermore, non-symplectic integrators for Hamiltonian ODEs inevitably introduce numerical instability in long time integrations.

Molecules are generally complex systems with many degrees of freedom. Their phase space is entangled with regular and chaotic regions. Therefore, it is not a surprise that computational chemistry has always tried to exploit new computer advances, such as parallel computing (including grid and cloud), as well as GPU technology. As far as the integration of the trajectories of large molecules for a long time is concerned, the good performance of several symplectic integrators has been recognized. Among the most popular ones are the simple low-order symplectic algorithms, such as the leapfrog or velocity Verlet algorithms [52].

Presently, AI (artificial intelligence) methods, particularly machine-learning techniques, are under exploration. In the last few years, an intense interest has been directed at finding algorithms that incorporate physics knowledge into neural networks. Hamiltonian neural networks label techniques that try to solve Hamilton’s equations with deep neural networks [53,54,55]. Solving the Schrödinger equation has also attracted the interest of researchers in this field [56]. However, we must note that at present, physics neural networks (PNN) algorithms are tested with low-dimensional systems, and their extension to systems with many degrees of freedom remains to be proved.

Irrespective of the final result of this endeavor, what is worth pursuing is developing computational methods that take into account multiscale physical theories that cover extended temporal/spatial scales. What we have shown in this review is that Hamiltonian (symplectic) geometry is the foundation of principal physical theories. This should be taken into account in building PNN algorithms [57], and it signals the need for new projects. Hamiltonian geometry yields the necessary and sufficient conditions for the mutual assistance of humans and machines in deep-learning processes.