1 package ffx.numerics.multipole;
2
3 import java.util.HashMap;
4
5 import static ffx.numerics.math.ScalarMath.binomial;
6 import static java.lang.String.format;
7
8 public class MultipoleUtilities {
9
10
11 /**
12 * Returns the number of tensors for derivatives to the given order.
13 *
14 * @param order maximum number of derivatives.
15 * @return the number of tensors.
16 * @since 1.0
17 */
18 public static int tensorCount(int order) {
19 long ret = binomial(order + 3, 3);
20 assert (ret < Integer.MAX_VALUE);
21 return (int) ret;
22 }
23
24 /**
25 * Convenience method for writing out tensor indices.
26 *
27 * @param l number of d/dx partial derivatives.
28 * @param m number of d/dx partial derivatives.
29 * @param n number of d/dx partial derivatives.
30 * @return a String of the form <code>Rlmn</code>.
31 */
32 protected static String rlmn(int l, int m, int n) {
33 return format("R%d%d%d", l, m, n);
34 }
35
36 /**
37 * Convenience method for writing out intermediate terms in the recursion.
38 *
39 * @param l number of d/dx partial derivatives.
40 * @param m number of d/dx partial derivatives.
41 * @param n number of d/dx partial derivatives.
42 * @return a String of the form <code>termLMN</code>.
43 */
44 protected static String term(int l, int m, int n) {
45 return format("term%d%d%d", l, m, n);
46 }
47
48 /**
49 * Convenience method for writing out intermediate terms in the recursion.
50 *
51 * @param l number of d/dx partial derivatives.
52 * @param m number of d/dx partial derivatives.
53 * @param n number of d/dx partial derivatives.
54 * @param j the jth intermediate term.
55 * @return a String of the form <code>termLMNJ</code>.
56 */
57 protected static String term(int l, int m, int n, int j) {
58 return format("term%d%d%d%d", l, m, n, j);
59 }
60
61 /**
62 * The index is based on the idea of filling tetrahedron.
63 * <p>
64 * 1/r has an index of 0.
65 * <br>
66 * derivatives of x are first; indices from 1..o for d/dx..(d/dx)^o
67 * <br>
68 * derivatives of x and y are second; base triangle of size (o+1)(o+2)/2
69 * <br>
70 * derivatives of x, y and z are last; total size (o+1)*(o+2)*(o+3)/6
71 * <br>
72 * <p>
73 * This function is useful to set up masking constants:
74 * <br>
75 * static int Tlmn = ti(l,m,n,order)
76 * <br>
77 * For example the (d/dy)^2 (1/R) storage location:
78 * <br>
79 * static int T020 = ti(0,2,0,order)
80 *
81 * @param dx int The number of d/dx operations.
82 * @param dy int The number of d/dy operations.
83 * @param dz int The number of d/dz operations.
84 * @param order int The maximum tensor order (0 .LE. dx + dy + dz .LE. order).
85 * @return int in the range (0..binomial(order + 3, 3) - 1)
86 */
87 protected static int ti(int dx, int dy, int dz, int order) {
88 if (dx < 0 || dy < 0 || dz < 0 || dx + dy + dz > order) {
89 return -1;
90 }
91
92 int size = (order + 1) * (order + 2) * (order + 3) / 6;
93 /*
94 We only get to the top of the tetrahedron if dz = order, otherwise
95 subtract off the top, including the level of the requested tensor
96 index.
97 */
98 int top = order + 1 - dz;
99 top = top * (top + 1) * (top + 2) / 6;
100 int zIndex = size - top;
101 /*
102 Given the "dz level", dy can range from (0 .. order - dz).
103 To get to the row for a specific value of dy, dy*(order + 1) - dy*(dy-1)/2 indices are skipped.
104 This is an operation that looks like the area of rectangle, minus the area of an empty triangle.
105 */
106 int yIndex = dy * (order - dz) - (dy - 1) * (dy - 2) / 2 + 1;
107 /*
108 Given the dz level and dy row, dx can range from (0..order - dz - dy)
109 The dx index is just walking down the dy row for "dx" steps.
110 */
111 return dx + yIndex + zIndex;
112 }
113
114 /**
115 * Convenience method for writing out tensor indices.
116 *
117 * @param l number of d/dx partial derivatives.
118 * @param m number of d/dx partial derivatives.
119 * @param n number of d/dx partial derivatives.
120 * @return a String of the form <code>lmn</code>.
121 */
122 protected static String lmn(int l, int m, int n) {
123 return format("%d%d%d", l, m, n);
124 }
125
126 /**
127 * Load a tensor element into a SIMD register.
128 * <p>
129 * For l=m=n=0, the returned String is:
130 * DoubleVector t000 = fromArray(SPECIES, t, T000);
131 *
132 * @param l The number of d/dx operations.
133 * @param n The number of d/dy operations.
134 * @param m The number of d/dz operations.
135 * @return a String of the form <code>DoubleVector t000 = fromArray(SPECIES, t, T000);</code>
136 */
137 protected static String loadTensor(int l, int n, int m, HashMap<Integer, String> tensorMap) {
138 String s = lmn(l, n, m);
139 String offset = "T" + lmn(l, n, m);
140 String name = "t" + s;
141 if (tensorMap.containsValue(name)) {
142 return "";
143 }
144 int size = tensorMap.size();
145 tensorMap.put(size, name);
146 return format("\tDoubleVector %s = fromArray(SPECIES, t, %s);\n", name, offset);
147 }
148
149 /**
150 * Code to store an electrostatic potential element into an array.
151 * <p>
152 * For l=m=n=0, the returned String is (including a preceding tab):
153 * term000.intoArray(to, T000);
154 *
155 * @param to variable name of the array to store into.
156 * @param l number of d/dx partial derivatives.
157 * @param n number of d/dy partial derivatives.
158 * @param m number of d/dz partial derivatives.
159 * @return a String of the form <code>term000.intoArray(to, T000);</code>
160 */
161 protected static String storePotential(String to, int l, int n, int m) {
162 String from = term(l, n, m);
163 String offset = "T" + lmn(l, n, m);
164 return format("\t%s.intoArray(%s, %s);\n", from, to, offset);
165 }
166
167 /**
168 * Code to store a negated electrostatic potential element into an array.
169 * <p>
170 * For l=m=n=0, the returned String is (including a preceding tab):
171 * term000.neg().intoArray(to, T000);
172 *
173 * @param to variable name of the array to store into.
174 * @param l number of d/dx partial derivatives.
175 * @param n number of d/dy partial derivatives.
176 * @param m number of d/dz partial derivatives.
177 * @return a String of the form <code>term000.intoArray(to, T000);</code>
178 */
179 protected static String storePotentialNeg(String to, int l, int n, int m) {
180 String from = term(l, n, m);
181 String offset = "T" + lmn(l, n, m);
182 return format("\t%s.neg().intoArray(%s, %s);\n", from, to, offset);
183 }
184 }