(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 4.0,
MathReader 4.0, or any compatible application. The data for the notebook
starts with the line containing stars above.
To get the notebook into a Mathematica-compatible application, do one of
the following:
* Save the data starting with the line of stars above into a file
with a name ending in .nb, then open the file inside the application;
* Copy the data starting with the line of stars above to the
clipboard, then use the Paste menu command inside the application.
Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode. Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).
NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing the
word CacheID, otherwise Mathematica-compatible applications may try to
use invalid cache data.
For more information on notebooks and Mathematica-compatible
applications, contact Wolfram Research:
web: http://www.wolfram.com
email: info@wolfram.com
phone: +1-217-398-0700 (U.S.)
Notebook reader applications are available free of charge from
Wolfram Research.
***********************************************************************)
(*CacheID: 232*)
(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[ 25340, 659]*)
(*NotebookOutlinePosition[ 26006, 683]*)
(* CellTagsIndexPosition[ 25962, 679]*)
(*WindowFrame->Normal*)
Notebook[{
Cell["\<\
Random walk in a 3d hexagonal close-packed (cannonball) \
lattice.\
\>", "Text"],
Cell[TextData[{
"Let each integer position on the ",
Cell[BoxData[
\(TraditionalForm\`x\)]],
" axis of a cartesian ",
Cell[BoxData[
\(TraditionalForm\`x, y, z\)]],
" coordinate system be a lattice point. In the ",
Cell[BoxData[
\(TraditionalForm\`x, y\)]],
" plane pack additional rows so that the new lattice points are eactly one \
unit from all their neighbors. An angle of 60 degrees and radius one, or ",
Cell[BoxData[
\(TraditionalForm\`\((1\/2, \@3\/2)\)\)]],
" in rectangular notation."
}], "Text"],
Cell[TextData[{
"In the third dimension we need to find the point in the ",
Cell[BoxData[
\(TraditionalForm\`z > 0\)]],
" half of the ",
Cell[BoxData[
\(TraditionalForm\`x = 1\/2\)]],
" plane that is exactly one unit from the origin and the point ",
Cell[BoxData[
\(TraditionalForm\`\((1\/2, \[Sqrt]3\/2)\)\)]],
"."
}], "Text"],
Cell[BoxData[{
\(\(p\_0 = {0, 0, 0};\)\), "\[IndentingNewLine]",
\(\(p\_x = {1, 0, 0};\)\), "\[IndentingNewLine]",
\(\(p\_y = {1\/2, y\_0, 0};\)\), "\[IndentingNewLine]",
\(\(p\_z = {1\/2, y\_z, z\_0};\)\), "\[IndentingNewLine]",
\(\(norm[v_] := Sqrt[v . v];\)\ (*\
real\ valued\ *) \), "\[IndentingNewLine]",
\(comp =
Last[Solve[{norm[p\_y - p\_x] \[Equal] 1, norm[p\_z - p\_x] \[Equal] 1,
norm[p\_z - p\_y] \[Equal] 1}, {y\_0, , y\_z, z\_0}]]\)}], "Input"],
Cell["Store the vectors explicitly ...", "Text"],
Cell[BoxData[{
\(p\_x = {1, 0, 0} /. comp\), "\[IndentingNewLine]",
\(p\_y = {1\/2, y\_0, 0} /. comp\), "\[IndentingNewLine]",
\(p\_z = {1\/2, y\_z, z\_0} /. comp\)}], "Input"],
Cell["\<\
These vectors are not orthogonal, so we cannot simply treat it as a \
coordinate system. It will be sufficient to have an algorithm for converting \
a path description in terms of integer multiples of these vectors to \
Cartesian coordinates. Vector addition will do, but we will develop an \
integer arithemtic-oriented algorithm.\
\>", "Text"],
Cell[TextData[{
"Let us construct a synthetic orthogonal grid coordinate system. In the ",
Cell[BoxData[
\(TraditionalForm\`x - y\)]],
" plane, the row of lattice points along the ",
Cell[BoxData[
\(TraditionalForm\`x\)]],
" axis is a one-dimensional unit. Along the synthetic ",
Cell[BoxData[
\(TraditionalForm\`y = 1\)]],
" line is a copy of this lattice row with each node displaced by one ",
Cell[BoxData[
\(TraditionalForm\`p\_y\)]],
" from the corresponding node in the base row. In the ",
Cell[BoxData[
\(TraditionalForm\`y = 2\)]],
" row the nodes are again aligned with those in the base row: the \
appropriate offset is:"
}], "Text"],
Cell[BoxData[
\(2 p\_y - p\_x\)], "Input"],
Cell[TextData[{
"A similar treatment applies in the ",
Cell[BoxData[
\(TraditionalForm\`z\)]],
" direction. The base ",
Cell[BoxData[
\(TraditionalForm\`x - y\)]],
" lattice plane is replicated at mutliple of ",
Cell[BoxData[
\(TraditionalForm\`p\_z\)]],
", with 3 steps along ",
Cell[BoxData[
\(TraditionalForm\`p\_z\)]],
" bringing us to a plane aligned exactly with the base plane:"
}], "Text"],
Cell[BoxData[
\(3 p\_z - p\_x - p\_y\)], "Input"],
Cell[TextData[{
"We can exploit the above to efficiently compute the location given a path \
description. It is sufficient to record the net motion in each of the \
lattice bond directions. The Cartesian coordinates for a lattice point at ",
Cell[BoxData[
\(TraditionalForm\`i\ p\_x + j\ p\_y + k\ p\_z\)]],
" can be efficiently computed as:"
}], "Text"],
Cell[BoxData[{
\(\(pathtoloc[i0_, j0_, k0_] :=
Module[{i, j, k}, \[IndentingNewLine]k =
k0/3; \[IndentingNewLine]j = \((j0 + k)\)/
2; \[IndentingNewLine]i = i0 + j + k; \[IndentingNewLine]{i,
j*\@3, k*\@6}\[IndentingNewLine]];\)\), "\[IndentingNewLine]",
\(\(pathtoloc[{i_, j_, k_}] := pathtoloc[i, j, k];\)\)}], "Input"],
Cell[BoxData[{
\(\(Clear[i, j, k];\)\), "\[IndentingNewLine]",
\(Simplify[
pathtoloc[i, j, k] \[Equal] i\ p\_x + j\ p\_y + k\ p\_z]\)}], "Input"],
Cell["\<\
Each lattice point has twelve neighbors at a distance of exactly \
one unit. These can be generated by taking the differences of all ordered \
pairs selected without replacement from the set consisting of the three bond \
direction vectors and the zero vector:\
\>", "Text"],
Cell[BoxData[{
\(\(nlist = {p\_x - 0, 0 - p\_x, p\_y - 0, 0 - p\_y, p\_z - 0,
0 - p\_z, \[IndentingNewLine]p\_x - p\_y, p\_y - p\_x, p\_y - p\_z,
p\_z - p\_y, p\_x - p\_z,
p\_z - p\_x};\)\), "\[IndentingNewLine]",
\(Simplify[Map[norm, nlist]]\)}], "Input"],
Cell[TextData[{
"We can denote one of these twelve neighbor vectors by a string of length \
three from the alphabet ",
Cell[BoxData[
\(TraditionalForm\`{\(-1\), 0, \(+1\)}\)]],
" where no two members of the string may have a positive product. The \
vector corresponding to the string ",
Cell[BoxData[
\(TraditionalForm\`{i, j, k}\)]],
" is then ",
Cell[BoxData[
\(TraditionalForm\`{i\ p\_x + j\ p\_y + k\ p\_z}\)]],
", which is necessarily then one of the vectors constructed above."
}], "Text"],
Cell[TextData[{
"A random walk in this latice can be generated by a Markov chain over these \
strings. Suppose we wish to impose a stochastic limit on the curvature of \
the walk in 3-space. Given the step vector that brought us to the current \
node, consider the projection (dot product) of each of the neighbor vectors \
onto this arrival vector. The vector pointing to the previous node will have \
a projection of ",
Cell[BoxData[
\(TraditionalForm\`\(-1\)\)]],
" while the arrival vector will of course have projection one. There are \
two vectors orthogonal to the arrival vector. Of the remaining eight \
vectors, four will have projection ",
Cell[BoxData[
\(TraditionalForm\`\(-\(1\/2\)\)\)]],
" and the other four will have projection ",
Cell[BoxData[
\(TraditionalForm\`1\/2\)]],
"."
}], "Text"],
Cell[BoxData[{
\(\(a = nlist[\([9]\)];\)\), "\[IndentingNewLine]",
\(Simplify[Map[a . # &, nlist]]\)}], "Input"],
Cell[TextData[{
"It is a relatively simple matter to tabulate the projection relation in \
terms of the vector description string and to assign transition probabilities \
based on projection equivalence classes. If angles more acute than ",
Cell[BoxData[
\(TraditionalForm\`120 \[Degree]\)]],
" are to be avoided one may assign a probability to going straight and \
divide the complement equally among the remaining four forward vectors."
}], "Text"],
Cell[TextData[{
"It is noteworthy that the four forward curve directions are not uniformly \
distributed about the arrival vector. If the observer is situated at the \
current node and looking in the direction of the arrival vector and oriented \
so that one of the forward curve vectors points toward the hoizon ",
Cell[BoxData[
\(TraditionalForm\`60 \[Degree]\)]],
" left of center, then another points to the horizon ",
Cell[BoxData[
FormBox[
FormBox[\(60 \[Degree]\),
"TraditionalForm"], TraditionalForm]]],
" right of center. Of the other two, one is above the horizon on one side \
while the last is below the horizon on the opposite side. This picture is \
not changed when stepping along the arrival vector and can be thought of as a \
kind of bias in the lattice; a bias against the underrepresented quadrants of \
the local sky and ground. This should have no macroscopic significance."
}], "Text"],
Cell[BoxData[
\(<< Graphics`Shapes`\)], "Input"],
Cell[BoxData[{
\(\(s = Sphere[1/2];\)\), "\[IndentingNewLine]",
\(\(t[g_, v_] := TranslateShape[g, v];\)\), "\[IndentingNewLine]",
\(\(uball = Map[t[s, #] &, nlist];\)\), "\[IndentingNewLine]",
\(\(Show[Graphics3D[uball], AspectRatio \[Rule] 1,
PlotLabel \[Rule]
StyleForm["\", Subsubsection]];\)\)}], "Input"],
Cell[BoxData[{
\(\(td[g_, d_] :=
TranslateShape[g, pathtoloc[d]];\)\), "\[IndentingNewLine]",
\(\(cstack =
Flatten[Table[
td[s, {i, j, k}], {i, 0, 3}, {j, 0, 3 - i}, {k, 0, 3 - i - j}],
2];\)\), "\[IndentingNewLine]",
\(\(Show[Graphics3D[cstack], AspectRatio \[Rule] 1,
PlotLabel \[Rule] "\"];\)\)}], "Input"],
Cell["\<\
Here is the view of the forward vector targets described \
above.\
\>", "Text"],
Cell[BoxData[{
\(\(sky =
Map[t[s, #] &,
Select[nlist, # . p\_x > 0 &]];\)\), "\[IndentingNewLine]",
\(\(Show[Graphics3D[sky], ViewPoint \[Rule] \(-p\_x\),
PlotLabel \[Rule]
StyleForm["\", Subsubsection]];\)\)}], "Input"],
Cell["\<\
In the simulation of vorticity filaments we need to know the \
midpoints of the vectors, so we will have our walk make two steps along each \
vector. We will parameterize our walk by the probability of going directly \
forward, with the complementary events equally divided among the four \
off-forward directions.\
\>", "Text"],
Cell[TextData[{
"We begin by tabulating the off-forward vectors. If our arrival vector is \
desctibed by the string ",
Cell[BoxData[
\(TraditionalForm\`{i, j, k}\)]],
" then the simplest way to find the foward vectors is to generate a list of \
all possible descritions and select those that have the desired projection."
}], "Text"],
Cell[BoxData[{
\( (*\ "\"\ will\ return\ a\ structure\ containing\ lists\ of\ \
vector\ descriptions\ organized\ by\ their\ dot\ product\ with\ the\ vector\ \
described\ by\ the\ parameter\ *) \[IndentingNewLine]\(strdot[a_, b_] :=
pathtoloc[a] . pathtoloc[b];\)\), "\[IndentingNewLine]",
\(\(alldirs =
Select[Flatten[
Table[{i, j, k}, {i, \(-1\), 1}, {j, \(-1\), 1}, {k, \(-1\), 1}],
2], strdot[#, #] \[Equal] 1 &];\)\), "\[IndentingNewLine]",
\(\(forward[
d_] := \[IndentingNewLine]{d, \[IndentingNewLine]Select[alldirs,
strdot[d, #] \[Equal] 1/2 &], \[IndentingNewLine]Select[alldirs,
strdot[d, #] \[Equal] 0 &], \[IndentingNewLine]Select[alldirs,
strdot[d, #] \[Equal] \(-1\)/
2 &], \[IndentingNewLine]\(-d\)};\)\), \
"\[IndentingNewLine]",
\(forward[{1, 0, 0}]\)}], "Input"],
Cell[BoxData[{
\(\(Clear[ftab];\)\), "\[IndentingNewLine]",
\(\(Do[
ftab[alldirs[\([i]\)]] = \(forward[alldirs[\([i]\)]]\)[\([2]\)], {i,
12}];\)\), "\[IndentingNewLine]",
\(ftab[{1, 0, 0}]\)}], "Input"],
Cell["\<\
Given an arrival vector in string form we generate the departure \
vector description as follows:\
\>", "Text"],
Cell[BoxData[
\(<< Statistics`DiscreteDistributions`\)], "Input"],
Cell[BoxData[{
\(\(u[lo_, hi_] :=
Random[DiscreteUniformDistribution[hi - lo + 1]] - 1 +
lo;\)\), "\[IndentingNewLine]",
\(\(nextdir[d_, p_] :=
If[Random[] \[LessEqual] p,
d, \(ftab[d]\)[\([u[1, 4]]\)]];\)\), "\[IndentingNewLine]",
\(nextdir[{1, 0, 0}, 1]\), "\[IndentingNewLine]",
\(nextdir[{1, 0, 0}, 0]\)}], "Input"],
Cell[TextData[{
"Given a starting position and direction, both in path description form, we \
generate a walk of length ",
Cell[BoxData[
\(TraditionalForm\`n\)]],
" as follows:"
}], "Text"],
Cell[BoxData[{
\( (*\
Each\ element\ is\ a\ center\ position\ and\ a\ \(\(direction\)\(.\)\)\ \
*) \[IndentingNewLine]\(stepfunc[{loc_, dir_}, p_] :=
Module[{newdir = nextdir[dir, p]}, \[IndentingNewLine]{loc + dir +
newdir, newdir}];\)\), "\[IndentingNewLine]",
\(\(hexwalk[sloc_, sdir_, p_, n_] :=
NestList[stepfunc[#, p] &, {sloc + sdir, sdir}, n];\)\)}], "Input"],
Cell[BoxData[{
RowBox[{\(<< Graphics`Shapes`\), "\[IndentingNewLine]",
RowBox[{"(*", " ",
RowBox[{
"build", " ", "an", " ", "arrow", " ", "two", " ", "units", " ",
"long", " ", "at", " ", Cell["(0,0,0)"], " ", "and", " ",
"oriented", " ", "along", " ", "the", " ", "vector", " ",
"described", " ", "by", " ", "d_"}], " ",
"*)"}]}], "\[IndentingNewLine]", \(vpic[d_] :=
Module[{x, y, z, arrow}, \[IndentingNewLine]{x, y, z} =
pathtoloc[
d]; \[IndentingNewLine]arrow = {TranslateShape[
Cylinder[ .08, .8, 8], {0, 0, \(- .2\)}],
TranslateShape[
Cone[ .3, .4, 10], {0,
0, .6}]}; \[IndentingNewLine]RotateShape[arrow, 0,
Pi/2 - ArcTan[norm[{x, y}], z],
Pi/2 - ArcTan[x,
y]]\[IndentingNewLine]];\), "\[IndentingNewLine]", \(vwalk[
w_] := Map[
TranslateShape[vpic[#[\([2]\)]], pathtoloc[#[\([1]\)]]] &,
w];\), "\[IndentingNewLine]", \(Show[
Graphics3D[vwalk[hexwalk[{0, 0, 0}, {1, 0, 0}, .2, 6]]],
Axes \[Rule] True];\)}], "Input"],
Cell["\<\
We will generate our walks by starting at the origin and going in a \
randomly selected initial direction, then offsetting the entire walk to a \
randomly selected inital point near the origin.\
\>", "Text"],
Cell[BoxData[{
\(\(rwalk[p_, n_, s_] :=
Module[{w, loc}, \[IndentingNewLine]w =
hexwalk[{0, 0, 0}, alldirs[\([u[1, 12]]\)], p,
n]; \[IndentingNewLine]For[
loc = {u[\(-s\), s], u[\(-s\), s], u[\(-s\), s]},
norm[pathtoloc[loc]] \[GreaterEqual] s,
loc = {u[\(-s\), s], u[\(-s\), s],
u[\(-s\),
s]}]; \[IndentingNewLine]loc -= \(w[\([Floor[
n/2]]\)]\)[\([1]\)]; \[IndentingNewLine]Map[{#[\([1]\)] +
loc, #[\([2]\)]} &,
w]\[IndentingNewLine]];\)\), "\[IndentingNewLine]",
\(rwalk[ .1, 3, 10]\)}], "Input"],
Cell["\<\
To visualize a pair of walks embedded in the lattice it is \
sufficient to concatenate their list representations and feed the result to \
vwalk:\
\>", "Text"],
Cell[BoxData[{
\(\(w1 = rwalk[ .1, 6, 6];\)\), "\[IndentingNewLine]",
\(w2 = rwalk[ .1, 6, 6];
Show[Graphics3D[vwalk[Join[w1, w2]]]];\)}], "Input"],
Cell[TextData[{
"With all vectors having length two the dot product of any pair will be one \
of ",
Cell[BoxData[
\(TraditionalForm\`0, \(\[PlusMinus]2\)\ or\ \[PlusMinus] 4\)]],
". We now tabulate these dot products by descriptions string:"
}], "Text"],
Cell[BoxData[{
\(\(Do[
dottab[alldirs[\([i]\)], alldirs[\([j]\)]] =
strdot[2 alldirs[\([i]\)], 2 alldirs[\([j]\)]], {i, 12}, {j,
12}];\)\), "\[IndentingNewLine]",
\(dottab[{1, 0, 0}, {1, 0, 0}]\)}], "Input"],
Cell["\<\
Given a pair of filaments (walks) in the lattice, we compute a \
statistic correspondiing to the interaction energy of the filaments. For \
each pair of vectors with one member chosen from each filament, the energy is \
the dot product of the vectors divided by the distance between their \
midpoints. We therefore begin our computation by partitioning the vectors of \
each filament into subsets by orientation. In the resulting representation \
the orientation information will be implicit.\
\>", "Text"],
Cell[BoxData[{
\(\(setrep[w_] :=
Table[Cases[w, {pos_, alldirs[\([i]\)]} \[Rule] pos], {i,
12}];\)\), "\[IndentingNewLine]",
\(setrep[{{{0, 0, 0}, {1, 0, 0}}, {{1, 0, 0}, {0, 1, 0}}}]\)}], "Input"],
Cell["\<\
Given two walks in set representation, we can compute the total \
energy by constructing all pairs of vectors grouped by orientation. There \
are 144 ways to choose one orientation from each filament, and we need a list \
of the distances between all pairs of vectors having the corresponding \
orientations. We will avoid constructing this list for those orientation \
pairs that have are orthogonal, since they contribute nothing to the \
energy.\
\>", "Text"],
Cell["\<\
The contribution of each vector pair within a set having identical \
orientations is inversely proportional to the distance between center points. \
To keep everything in integer form we will store the norm squared rather \
than the vector or its norm in for each pair.\
\>", "Text"],
Cell[BoxData[{
\(Remove[difnorm2]\), "\[IndentingNewLine]",
\(\(difnorm2[v1_, v2_] :=
Module[{v = pathtoloc[v2 - v1]}, v . v];\)\[IndentingNewLine] (*\
Note : \ this\ could\ be\ a\ lot\ faster\ if\ we\ expanded\ pathtoloc\ \
inline\ and\ compiled\ it\ *) \)}], "Input"],
Cell[BoxData[
\(\(spectrum[s1_, s2_] :=
Module[{d1, d2, dot, spec}, \[IndentingNewLine]spec =
Table[0, {i, 200^2}]; \[IndentingNewLine]Do[
d1 = alldirs[\([i]\)]; \[IndentingNewLine]d2 =
alldirs[\([j]\)]; \[IndentingNewLine]dot =
dottab[d1, d2]; \[IndentingNewLine]If[dot != 0\ ,
Scan[\((spec[\([#]\)] += dot)\) &,
Outer[difnorm2, s1[\([i]\)], s2[\([j]\)], 1], {2}]], {i,
12}, {j,
12}\[IndentingNewLine]]; \[IndentingNewLine]spec\
\[IndentingNewLine]];\)\)], "Input"],
Cell["\<\
Let's have a look at the spectrum for the pair of small filaments \
in the figure above.\
\>", "Text"],
Cell[BoxData[{
\(\(s1 = setrep[w1];\)\), "\[IndentingNewLine]",
\(s2 = setrep[w2]\)}], "Input"],
Cell[BoxData[
\(Take[spectrum[s1, s2], 100]\)], "Input"],
Cell[BoxData[
\(s = spectrum[s1, s2];
ListPlot[Table[Sum[s[\([i]\)]/Sqrt[i], {i, 1, n}], {n, 1, 300, 10}],
PlotJoined \[Rule] True];\)], "Input"],
Cell["\<\
To compute the total energy we need to divide the bin counts in \
this spectrum by the reciprocal of the corresponding norm. Since the table \
is indexed by the square of the norm, let's build a table of the reciprocal \
square roots:\
\>", "Text"],
Cell[BoxData[
\(\(rnormtab = Table[N[1/Sqrt[i]], {i, 200^2}];\)\)], "Input"],
Cell["\<\
The energy then is given by the sum of the elements of the \
term-by-term product of these two lists. Unfortunately this list is sorted \
with the largest elements first, statistically speaking, so it would be more \
numerically sensible to add it up in reverse order.\
\>", "Text"],
Cell[BoxData[{
\(\(energy[s_] :=
Plus @@ Reverse[N[s]*rnormtab];\)\), "\[IndentingNewLine]",
\(\(energy[w1_, w2_] :=
energy[spectrum[setrep[w1],
setrep[w2]]];\)\), "\[IndentingNewLine]",
\(energy[w1, w2] // Timing\)}], "Input"],
Cell["\<\
Let's see how long it takes to compute the energy for a larger \
scale example:\
\>", "Text"],
Cell[BoxData[{
\(\(lw1 = rwalk[ .6, 201, 25];\)\), "\[IndentingNewLine]",
\(\(lw2 = rwalk[ .6, 201, 25];\)\), "\[IndentingNewLine]",
\(energy[lw1, lw2] // Timing\)}], "Input"],
Cell["Let's see where the time is going:", "Text"],
Cell[BoxData[{
\(\({ls1, ls2} = {setrep[lw1], setrep[lw2]};\) //
Timing\), "\[IndentingNewLine]",
\(\(s = spectrum[ls1, ls2];\) // Timing\), "\[IndentingNewLine]",
\(energy[s] // Timing\)}], "Input"],
Cell["\<\
Let's try to speed this up by building an optimized version of the \
difnorm2 function. Here is the core of pathtoloc and the norm squared of the \
result:\
\>", "Text"],
Cell[BoxData[{
\(\(k = k0/3;\)\), "\[IndentingNewLine]",
\(\(j = \((j0 + k)\)/2;\)\), "\[IndentingNewLine]",
\(\(i = i0 + j + k;\)\), "\[IndentingNewLine]",
\({i, j*\@3, k*\@6}\), "\[IndentingNewLine]",
\(Simplify[% . %]\)}], "Input"],
Cell["As we observed, this is an integer expression.", "Text"],
Cell[BoxData[
\(\(\( (*\ \(difnorm2[v1_, v2_] :=
Module[{v = pathtoloc[v2 - v1]},
v . v];\)\ *) \)\(\[IndentingNewLine]\)\(difnorm2 =
Compile[{{v1, _Integer, 1}, {v2, _Integer, 1}},
Module[{i0, j0, k0}, \[IndentingNewLine]{i0, j0, k0} =
v2 - v1; \[IndentingNewLine]i0^2 + j0^2 + j0\ k0 + k0^2 +
i0 \((j0 +
k0)\)\[IndentingNewLine]], {{i0, _Integer}, {j0, \
_Integer}, {k0, _Integer}}];\)\)\)], "Input"],
Cell[BoxData[{
\(\(s = spectrum[ls1, ls2];\) // Timing\), "\[IndentingNewLine]",
\(energy[s] // Timing\)}], "Input"],
Cell[TextData[{
"To speed this up further will require a different approach to accumulating \
the spectrum or directly computing the sum without the intermediate spectral \
representation. I have not yet found a better algorithm in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
". Note that the construction of the array of squared distances is only \
half the problem:"
}], "Text"],
Cell[BoxData[
\(Do[foo =
Sort[Flatten[Outer[difnorm2, ls1[\([i]\)], ls2[\([j]\)], 1]]], {i,
12}, {j, 12}] // Timing\)], "Input"],
Cell[BoxData[{
\(\(allpairs =
Flatten[Array[List, {12, 12}], 1];\)\), "\[IndentingNewLine]",
\(\(dotpairs =
Table[\[IndentingNewLine]If[
dot == 0, {}, \[IndentingNewLine]Select[
allpairs, \((dottab[alldirs[\([#[\([1]\)]]\)],
alldirs[\([#[\([2]\)]]\)]] \[Equal]
dot)\) &]\[IndentingNewLine]], \[IndentingNewLine]{dot, \
\(-4\), 4}];\)\), "\[IndentingNewLine]",
\(\(energy[s1_, s2_] :=
Module[{huge, bydot, pos, neg}, \[IndentingNewLine]huge =
Outer[Flatten[Outer[difnorm2, #1, #2, 1], 1] &, s1, s2,
1]; \[IndentingNewLine]bydot =
Map[Flatten[#, 1] &,
Map[huge[\([#[\([1]\)], #[\([2]\)]]\)] &, dotpairs, {2}],
1]; \[IndentingNewLine]2 \((Plus @@
rnormtab[\([bydot[\([7]\)]]\)] -
Plus @@ rnormtab[\([bydot[\([3]\)]]\)])\) +
4 \((Plus @@ rnormtab[\([bydot[\([9]\)]]\)] -
Plus @@
rnormtab[\([bydot[\([1]\)]]\)])\)\[IndentingNewLine]];\)\)\
}], "Input"],
Cell[BoxData[
\(energy[ls1, ls2] // Timing\)], "Input"],
Cell["\<\
This version is the speed champion, but it does not create a \
spectrum. Let's see if we can get a spectrum with reasonable additional \
overhead:\
\>", "Text"],
Cell[BoxData[
\(\(spectrum[s1_, s2_] :=
Module[{spec, huge, bydot, pos, neg}, \[IndentingNewLine]spec =
Table[0, {i, 200^2}]; \[IndentingNewLine]huge =
Outer[Flatten[Outer[difnorm2, #1, #2, 1], 1] &, s1, s2,
1]; \[IndentingNewLine]bydot =
Map[Flatten[#, 1] &,
Map[huge[\([#[\([1]\)], #[\([2]\)]]\)] &, dotpairs, {2}],
1]; \[IndentingNewLine]pos =
Split[Sort[
Join @@ bydot[\([{7, 9, 9}]\)]]]; \[IndentingNewLine]neg =
Split[Sort[
Join @@ bydot[\([{1, 1,
3}]\)]]]; \
\[IndentingNewLine]Scan[\((spec[\([#[\([1]\)]]\)] += Length[#])\) &,
pos]; \[IndentingNewLine]Scan[\((spec[\([#[\([1]\)]]\)] -=
Length[#])\) &, neg]; \[IndentingNewLine]2*
spec\[IndentingNewLine]];\)\)], "Input"],
Cell[BoxData[{
\(\(s = spectrum[ls1, ls2];\) // Timing\), "\[IndentingNewLine]",
\(energy[s] // Timing\)}], "Input"],
Cell["\<\
Well, at least it's correct. It isn't appreciably faster than the \
original, though.\
\>", "Text"]
},
FrontEndVersion->"4.0 for X",
ScreenRectangle->{{0, 800}, {0, 600}},
WindowSize->{601, 544},
WindowMargins->{{30, Automatic}, {-24, Automatic}},
CellLabelAutoDelete->True
]
(***********************************************************************
Cached data follows. If you edit this Notebook file directly, not using
Mathematica, you must remove the line containing CacheID at the top of
the file. The cache data will then be recreated when you save this file
from within Mathematica.
***********************************************************************)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[1717, 49, 90, 3, 32, "Text"],
Cell[1810, 54, 550, 15, 82, "Text"],
Cell[2363, 71, 359, 11, 58, "Text"],
Cell[2725, 84, 515, 9, 180, "Input"],
Cell[3243, 95, 48, 0, 32, "Text"],
Cell[3294, 97, 190, 3, 98, "Input"],
Cell[3487, 102, 357, 6, 68, "Text"],
Cell[3847, 110, 698, 18, 86, "Text"],
Cell[4548, 130, 47, 1, 27, "Input"],
Cell[4598, 133, 440, 14, 50, "Text"],
Cell[5041, 149, 54, 1, 27, "Input"],
Cell[5098, 152, 368, 7, 68, "Text"],
Cell[5469, 161, 379, 6, 136, "Input"],
Cell[5851, 169, 160, 3, 43, "Input"],
Cell[6014, 174, 285, 5, 68, "Text"],
Cell[6302, 181, 298, 5, 59, "Input"],
Cell[6603, 188, 528, 13, 86, "Text"],
Cell[7134, 203, 843, 18, 124, "Text"],
Cell[7980, 223, 121, 2, 43, "Input"],
Cell[8104, 227, 464, 8, 86, "Text"],
Cell[8571, 237, 955, 17, 140, "Text"],
Cell[9529, 256, 52, 1, 27, "Input"],
Cell[9584, 259, 412, 7, 123, "Input"],
Cell[9999, 268, 393, 8, 107, "Input"],
Cell[10395, 278, 89, 3, 32, "Text"],
Cell[10487, 283, 287, 6, 59, "Input"],
Cell[10777, 291, 339, 6, 68, "Text"],
Cell[11119, 299, 346, 7, 68, "Text"],
Cell[11468, 308, 918, 16, 235, "Input"],
Cell[12389, 326, 235, 5, 59, "Input"],
Cell[12627, 333, 121, 3, 32, "Text"],
Cell[12751, 338, 69, 1, 27, "Input"],
Cell[12823, 341, 378, 8, 75, "Input"],
Cell[13204, 351, 201, 6, 50, "Text"],
Cell[13408, 359, 416, 7, 107, "Input"],
Cell[13827, 368, 1191, 23, 209, "Input"],
Cell[15021, 393, 217, 4, 50, "Text"],
Cell[15241, 399, 671, 13, 139, "Input"],
Cell[15915, 414, 169, 4, 50, "Text"],
Cell[16087, 420, 161, 3, 59, "Input"],
Cell[16251, 425, 267, 6, 50, "Text"],
Cell[16521, 433, 248, 5, 59, "Input"],
Cell[16772, 440, 519, 8, 104, "Text"],
Cell[17294, 450, 228, 4, 43, "Input"],
Cell[17525, 456, 474, 8, 86, "Text"],
Cell[18002, 466, 294, 5, 68, "Text"],
Cell[18299, 473, 294, 5, 75, "Input"],
Cell[18596, 480, 596, 11, 171, "Input"],
Cell[19195, 493, 112, 3, 32, "Text"],
Cell[19310, 498, 104, 2, 43, "Input"],
Cell[19417, 502, 60, 1, 27, "Input"],
Cell[19480, 505, 161, 3, 59, "Input"],
Cell[19644, 510, 259, 5, 68, "Text"],
Cell[19906, 517, 80, 1, 27, "Input"],
Cell[19989, 520, 293, 5, 68, "Text"],
Cell[20285, 527, 272, 6, 59, "Input"],
Cell[20560, 535, 103, 3, 32, "Text"],
Cell[20666, 540, 189, 3, 59, "Input"],
Cell[20858, 545, 50, 0, 32, "Text"],
Cell[20911, 547, 220, 4, 59, "Input"],
Cell[21134, 553, 180, 4, 50, "Text"],
Cell[21317, 559, 258, 5, 98, "Input"],
Cell[21578, 566, 62, 0, 32, "Text"],
Cell[21643, 568, 501, 9, 107, "Input"],
Cell[22147, 579, 125, 2, 43, "Input"],
Cell[22275, 583, 397, 8, 68, "Text"],
Cell[22675, 593, 154, 3, 43, "Input"],
Cell[22832, 598, 1122, 22, 235, "Input"],
Cell[23957, 622, 59, 1, 27, "Input"],
Cell[24019, 625, 171, 4, 50, "Text"],
Cell[24193, 631, 902, 17, 187, "Input"],
Cell[25098, 650, 125, 2, 43, "Input"],
Cell[25226, 654, 110, 3, 32, "Text"]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)