Skip to main content

differential equations - Minimum energy path of a potential energy surface


I have calculated the energies of a molecule by varying two parameters (Data points are here data). I want to find the minimum energy path from the point (180.0, 179.99, 21.132) to (124.5, 124.49, 0) using Mathematica.



Answer




Let φ:RdR denote a potential function (e.g., from OP's data file).


We attempt to solve the system {γ(0)=p,grad(φ)|γ(t)=λ(t)γ(t)for all t[0,1],grad(φ)|γ(1)=0,|γ(t)|=μfor all t[0,1],μ=L(γ).


Here γ:[0,1]Rd is the curve that we are looking for, pRd is the starting point, L(γ) denotes arc length of the curve, and μR and λ:[0,1]R are dummy variables. For the interested reader who wonders why the last two equations are not joined to the single equation |γ(t)|=L(γ): This is done in order to enforce that the linearization of the discretized system (see below) is a sparse matrix. Moreover, these two equations are added for stability and in order to have as many degrees of freedom as equations in the discretized setting.



Next, we discretize the system. We replace the curve γ by a polygonal line with n edges and vertex coordinates given by x1,x2,,xn+1 in Rd and the function λ by a sequence λ1,,λn (one may think of this as function that is piecewise constant on edges). Replacing derivatives γ(t) by finite differences xi+1xin and by evaluating grad(φ) on the first n1 edge midpoints, we arrive at the discretized system



{x1=p,grad(φ)|(xi+xi+1)/2=λixi+1xinfor all i{1,,n1},grad(φ)|xn+1=0,|xi+1xin|=μfor all i{1,,n1},μ=ni=1|xi+1xi|.



These are routines for the system and its linearization. While ϕ denotes the potential, and DDϕ denote its first and second derivative, respectively. n is the number of edges and d is the dimension of the Euclidean space at hand.


The code is already optized for performance (direct initialization of SparseArray from CRS data), so hard to read. Sorry.


ClearAll[F, DF];
F[X_?VectorQ] :=
Module[{x, λ, μ, p1, p2, tangents, midpts, speeds},
x = Partition[X[[1 ;; (n + 1) d]], d];
λ = X[[(n + 1) d + 1 ;; -2]];
μ = X[[-1]];

p1 = Most[x];
p2 = Rest[x];
tangents = (p1 - p2) n;
midpts = 0.5 (p1 + p2);
speeds = Sqrt[(tangents^2).ConstantArray[1., d]];
Join[
x[[1]] - p,
Flatten[Dϕ @@@ Most[midpts] - λ Most[tangents]],
Dϕ @@ x[[-1]],
Most[speeds] - μ,

{μ - Total[speeds]/n}
]
];

DF[X_?VectorQ] :=
Module[{x, λ, μ, p1, p2, tangents, unittangents, midpts, id, A1x, A2x, A3x, A4x, A5x, buffer, A2λ, A, A4μ, A5μ},
x = Partition[X[[1 ;; (n + 1) d]], d];
λ = X[[(n + 1) d + 1 ;; -2]];
μ = X[[-1]];
p1 = Most[x];

p2 = Rest[x];
tangents = (p1 - p2) n;
unittangents = tangents/Sqrt[(tangents^2).ConstantArray[1., d]];
midpts = 0.5 (p1 + p2);
id = IdentityMatrix[d, WorkingPrecision -> MachinePrecision];
A1x = SparseArray[Transpose[{Range[d], Range[d]}] -> 1., {d, d (n + 1)}, 0.];
buffer = 0.5 DDϕ @@@ Most[midpts];
A2x = With[{
rp = Range[0, 2 d d (n - 1), 2 d],
ci = Partition[Flatten[Transpose[ConstantArray[Partition[Range[d n], 2 d, d], d]]], 1],

vals = Flatten@Join[
buffer - λ ConstantArray[id n, n - 1],
buffer + λ ConstantArray[id n, n - 1],
3
]
},
SparseArray @@ {Automatic, {d (n - 1), d (n + 1)}, 0., {1, {rp, ci}, vals}}];
A3x = SparseArray[Flatten[Table[{i, n d + j}, {i, 1, d}, {j, 1, d}], 1] -> Flatten[DDϕ @@ x[[-1]]], {d, d (n + 1)}, 0.];
A = With[{
rp = Range[0, 2 d n, 2 d],

ci = Partition[Flatten[Partition[Range[ d (n + 1)], 2 d, d]], 1],
vals = Flatten[Join[unittangents n, -n unittangents, 2]]
},
SparseArray @@ {Automatic, {n, d (n + 1)}, 0., {1, {rp, ci}, vals}}];
A4x = Most[A];
A5x = SparseArray[{-Total[A]/n}];
A2λ = With[{
rp = Range[0, d (n - 1)],
ci = Partition[Flatten[Transpose[ConstantArray[Range[n - 1], d]]], 1],
vals = Flatten[-Most[tangents]]

},
SparseArray @@ {Automatic, {d (n - 1), n - 1}, 0., {1, {rp, ci}, vals}}
];
A4μ = SparseArray[ConstantArray[-1., {n - 1, 1}]];
A5μ = SparseArray[{{1.}}];
ArrayFlatten[{
{A1x, 0., 0.},
{A2x, A2λ, 0.},
{A3x, 0., 0.},
{A4x, 0., A4μ},

{A5x, 0., A5μ}
}]
];


Let's load the data and make an interpolation function from it. Because DF assumes that ϕ is twice differentiable, we use spline interpolation of order 5.


data = Developer`ToPackedArray[N@Rest[
Import[FileNameJoin[{NotebookDirectory[], "data-e.txt"}], "Table"]
]];


p = {180.0, 179.99};
q = {124.5, 124.49};

Block[{x, y},
ϕ = Interpolation[data, Method -> "Spline", InterpolationOrder -> 5];
Dϕ = {x, y} \[Function] Evaluate[D[ϕ[x, y], {{x, y}, 1}]];
DDϕ = {x, y} \[Function] Evaluate[D[ϕ[x, y], {{x, y}, 2}]];
];

Now, we create some initial data for the system: The straight line from p to q as initial guess for x and some educated guesses for λ and μ. I also perturb the potential minimizer q a bit, because that seems to help FindRoot.



d = Length[p];
n = 1000;
tlist = Subdivide[0., 1., n];
xinit = KroneckerProduct[(1 - tlist), p] + KroneckerProduct[tlist , q + RandomReal[{-1, 1}, d]];

p1 = Most[xinit];
p2 = Rest[xinit];
tangents = (p1 - p2) n;
midpts = 0.5 (p1 + p2);


λinit = Norm /@ (Dϕ @@@ Most[midpts])/(Norm /@ Most[tangents]);
μinit = Norm[p - q];
X0 = Join[Flatten[xinit], λinit, {μinit}];

Throwing FindRoot onto it:


sol = FindRoot[F[X] == 0. F[X0], {X, X0}, Jacobian -> DF[X]]; // AbsoluteTiming //First
Xsol = (X /. sol);
Max[Abs[F[Xsol]]]



0.120724


6.15813*10^-10



Partitioning the solution and plotting the result:


xsol = Partition[Xsol[[1 ;; (n + 1) d]], d];
λsol = Xsol[[(n + 1) d + 1 ;; -2]];
μsol = Xsol[[-1]];
curve = Join[xsol, Partition[ϕ @@@ xsol, 1], 2];

{x0, x1} = MinMax[data[[All, 1]]];

{y0, y1} = MinMax[data[[All, 2]]];

Show[
Graphics3D[{Thick, Line[curve]}],
Plot3D[ϕ[x, y], {x, x0, x1}, {y, y0, y1}]
]

enter image description here



I'd like to express my scepticism about curves of steepest descent being a meaningful concept for the application in chemistry. Contrary to the notion of derivative Dφ, the notion of a gradient grad(φ) requires the additional structure of a Riemannian metric. In a nutshell, for detemining the slope of a function φ at a point x, you do not only need a way to measure height differences (those are provided by φ itself), but you need also a way of measuring (infinitesimal) lengths in the x-space; this is because of the elementary formula slope=Δφ|Δx|.



So far, we used the Euclidean metric on the parameter space Rd as Riemannian metric. This is however very arbitrary and depends on the choice of coordinates for the parameters. That means: Whenever different coordinates ˜x=Ψ(x) for parameters are chosen, the curve of steepst descent will change in a noncovariant way (see covariance principle). The covariance principle requires that after a change of coordinates, the resulting curve ˜γ satisfies the following equation (maybe up to reparameterization of the curves):


˜γ(t)=Ψ(γ(t)).


One can easily check that the curves of steepest descent with respect to the Euclidean metric does violate the covariance principle. For example, we may use the simple linear transformation Ψ(x1,x2)=(x1,2x2) (this amounts to using ψ(y1,y2)=φ(y1,y2/2) instead of φ in the algorithms described above and transforming the resulting curve back with Ψ) results in the following red curve of steepest descent, which is very different to the original one (black)


enter image description here


Once more, I'd like to refer to Quapp - Analysis of the concept of minimum energy path on the potential energy surface of chemically reacting systems where the principle of covariance is discussed more thoroughly.


Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...