This script does not implement any known saddle point search schemes such as the Dimer method or Nudged Elastic Band (NEB) method. It employs a unique approach based on the Euler-Lagrange equation to find the MEP directly by minimizing the action integral.
This repository contains a Python script to compute the Minimum Energy Path (MEP) between two points on a 2D potential energy surface. Using principles from classical mechanics, specifically the Euler-Lagrange equation, this method finds the path of least action through a saddle point.
Unlike traditional methods (e.g., Nudged Elastic Band or Dimer), this approach directly minimizes the action integral, avoiding artificial forces or constraints.
The script allows users to:
- Define their own potential energy function or use a default.
- Compute the MEP by minimizing the action integral.
- Ensure the path passes through a saddle point, where the gradient of the potential energy is zero, and the Hessian has both positive and negative eigenvalues.
The Potential Energy Surface (PES) represents the energy of a system as a function of its spatial coordinates. In two dimensions, it is expressed as V(x, y). Finding the MEP on a PES is crucial for understanding reaction pathways, transition states, and energy barriers.
The Lagrangian is defined as the difference between kinetic and potential energy:
L = T - V
For this problem:
- Kinetic Energy:
T = 0.5 * (dy/dx)^2 - Potential Energy:
V(x, y)
Thus, the Lagrangian becomes:
L(x, y, y') = 0.5 * (dy/dx)^2 + V(x, y)
The action S is the integral of the Lagrangian:
S = ∫[x0 to xf] L(x, y, y') dx
The path y(x) that minimizes S is obtained by solving the Euler-Lagrange equation:
d/dx (∂L/∂y') - ∂L/∂y = 0
Substituting L(x, y, y'), this simplifies to:
d^2y/dx^2 - ∂V/∂y = 0
A saddle point satisfies:
-
Gradient is zero:
∇V(xs, ys) = (∂V/∂x, ∂V/∂y) = (0, 0) -
Determinant of the Hessian matrix is negative:
det(H) = (∂^2V/∂x^2)(∂^2V/∂y^2) - (∂^2V/∂x∂y)^2 < 0
- Discretize the path between start
(x0, y0)and end(xf, yf). - Initial Guess for
y-coordinates (including saddle point) is provided. - Use spline interpolation to ensure smoothness and differentiability of
y(x).
Minimize the action S while enforcing:
- Gradient Constraint:
∇V(xs, ys) = (0, 0) - Hessian Determinant Constraint:
det(H) < 0
- Lagrangian:
L(x, y, y') = 0.5 * (dy/dx)^2 + V(x, y) - Action:
S = ∫[x0 to xf] L(x, y, y') dx - Euler-Lagrange Equation:
d^2y/dx^2 - ∂V/∂y = 0
- Python 3.x
- NumPy
- SciPy
- Matplotlib
- SymPy
Install dependencies:
pip install numpy scipy matplotlib sympy
git clone https://github.com/yourusername/mep-euler-lagrange.git
cd mep-euler-lagrange
python mep_finder.pyWhen prompted, enter the potential energy function or press Enter to use the default:
V(x, y) = -sin(πx) * sin(πy)Running the script with the default potential:
python mep_finder.pyOutput:
A contour plot of the PES with the MEP and saddle point. Console details: saddle point coordinates, gradient, and Hessian determinant.
Direct Minimization: No artificial forces or constraints. No Spring Forces: Avoids constructs like NEB's spring forces. Flexibility: Adapts to any user-defined potential energy function. Saddle Point Identification: Explicitly searches for saddle points.
Computational Cost: High for complex potentials or many discretization points. Dimensionality: Current implementation is limited to 2D. Initial Guess Sensitivity: Optimization success depends on initial guess.
Contributions are welcome! Feel free to submit pull requests or open issues.
Number of iterations: 368, function evaluations: 13545, CG iterations: 1037, optimality: 4.83e+00, constraint violation: 6.99e-21, execution time: 2.9 s. Saddle Point Coordinates: x = -7.077349822562493e-22, y = -7.077350130900838e-22 Gradient at saddle point: [6.9850646e-21 6.9850643e-21] Determinant of Hessian at saddle point: -97.40909103400243