       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

(
{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:= Timing[Reduce[CompleteEqnTab2, Flatten[xvec], Reals]]

Out= {18.527, x1 == 0 && y1 == 0 && z1 == 1/Sqrt && x2 == 0 &&
y2 == 1/Sqrt &&
z2 == 0 && x3 == 1/Sqrt && y3 == 0 && z3 == 0 && x4 == 1/
Sqrt &&
y4 == 1/Sqrt && z4 == 1/Sqrt && x5 == -(1/(3 Sqrt)) &&
y5 == (2 Sqrt)/3 && z5 == (2 Sqrt)/3 && x6 == -(11/(9
Sqrt)) &&
y6 == 5/(9 Sqrt) && z6 == 5/(9 Sqrt)}

Similar issues hold for larger clusters; Mathematica solves the less-
symmetric systems and hangs on the more-symmetric ones.  Is there
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

```

• Prev by Date: Re: More /.{I->-1} craziness
• Next by Date: Re: Re: More /.{I->-1} craziness
• Previous by thread: Re: Labeling points on ListPlot - Transform string date element in list of tripes
• Next by thread: Re: Trouble with coupled quadratic equations where the