Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
- To: mathgroup at smc.vnet.net
- Subject: [mg106476] Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
- From: Robert Hoy <robert.hoy at yale.edu>
- Date: Thu, 14 Jan 2010 05:46:47 -0500 (EST)
Hi, all! I am trying to solve for the structures of small atomic
clusters using Mathematica. Using adjacency (contact) matrices as
input, my the desired output is 3-D Cartesian coordinates. I
formulate the problem as coupled quadratic equations and
inequalities. For example, for clusters of M = 6 atoms I introduce
the variables
xvec = {{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}, {x5, y5,
z5}, {x6, y6, z6}}
Then, for an adjacency matrix A
{
{0, 1, 0, 1, 1, 1},
{1, 0, 1, 1, 0, 1},
{0, 1, 0, 1, 1, 1},
{1, 1, 1, 0, 1, 0},
{1, 0, 1, 1, 0, 1},
{1, 1, 1, 0, 1, 0}
}
there are M(M-1)/2 coupled equations and inequalities:
EqnTab = {(x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 ==
1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 >
1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 >
1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 ==
1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 ==
1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}
In these equations and inequalities, the left hand sides are
Dot[xvec[[m]] - xvec[[n]], xvec[[m]] - xvec[[n]] (i. e. the squared
Euclidean distance between xvec[[m]] and xvec[[n]] and the 'right hand
sides' are '==1' if A[[m,n]]==1 and '> 1' if A[[m,n]]=0. Only terms
above the diagonal of the adjacency matrix (n > m) are necessary since
A is symmetric.
I know the solution to these equations; the six elements of xvec form
the vertices of an octahedron.
Mathematica hangs when I try to solve this set. This isn't
surprising, because the above equations alone are invariant under
translations and rotations, so the solutions are continuously
degenerate. To eliminate the continuous degeneracies I add equations
constraining the first three 'atoms' to a plane (I can do this because
they are noncollinear). There is another discrete mirror degeneracy
(the other atoms can lie on either side of the defined plane), so I
add another inequality saying the fourth atom is on one side. This
adds ten equations/inequalities and gives the system
CompleteEqnTab = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 > 0,
y3 == 0, z3 == 0, x3 > 0,
x4/x3 + y4/y2 + z4/z1 > 1, [the 15 terms of EqnTab]}
Now, the problem is, I have been unable (after a considerable amount
of effort) to solve this system with Mathematica.
Reduce[CompleteEqnTab, Flatten[xvec], Reals] fails. So does
FindInstance[CompleteEqnTab, Flatten[xvec], Reals]. The memory used
by the kernel just increases (well past 1GB) until the program hangs.
I think the problem is that the solution for xvec (i. e. the vertices
of an octahedron) has many rotational symmetries - its point group is
O_h. Perhaps these symmetries give Mathematica's solver trouble?
In contrast, the exact same approach for a less symmetric cluster with
adjacency matrix
(
{0, 1, 1, 1, 1, 1},
{1, 0, 1, 1, 1, 1},
{1, 1, 0, 1, 0, 0},
{1, 1, 1, 0, 1, 0},
{1, 1, 0, 1, 0, 1},
{1, 1, 0, 0, 1, 0}
)
and equation/inequality set
CompleteEqnTab2 = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 > 0,
y3 == 0, z3 == 0, x3 > 0,
x4/x3 + y4/y2 + z4/z1 > 1, (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 ==
1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 ==
1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 ==
1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 >
1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 >
1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}
gives a quick solution
In[32]:= Timing[Reduce[CompleteEqnTab2, Flatten[xvec], Reals]]
Out[32]= {18.527, x1 == 0 && y1 == 0 && z1 == 1/Sqrt[2] && x2 == 0 &&
y2 == 1/Sqrt[2] &&
z2 == 0 && x3 == 1/Sqrt[2] && y3 == 0 && z3 == 0 && x4 == 1/
Sqrt[2] &&
y4 == 1/Sqrt[2] && z4 == 1/Sqrt[2] && x5 == -(1/(3 Sqrt[2])) &&
y5 == (2 Sqrt[2])/3 && z5 == (2 Sqrt[2])/3 && x6 == -(11/(9
Sqrt[2])) &&
y6 == 5/(9 Sqrt[2]) && z6 == 5/(9 Sqrt[2])}
Similar issues hold for larger clusters; Mathematica solves the less-
symmetric systems and hangs on the more-symmetric ones. Is there
anything I can do about this? I thought FindInstance was designed for
this sort of thing but it doesn't seem to help.
As you can see from the above, the form of the equations and
inequalities is simple. I assume Mathematica can handle the symmetric
cases and I'm just unaware of the necessary functionality. Any help
would be greatly appreciated. I am using Mathematica 7.0.1.0 on Mac
OS 10.5.8.
Thanks,
Robert Hoy
- Follow-Ups:
- Re: SymbolicPolynomialMod
- From: danl@wolfram.com
- Re: Trouble with coupled quadratic equations where the
- From: Daniel Lichtblau <danl@wolfram.com>
- Re: SymbolicPolynomialMod