Java tutorial
/** * Title: Force Field X. * * Description: Force Field X - Software for Molecular Biophysics. * * Copyright: Copyright (c) Michael J. Schnieders 2001-2017. * * This file is part of Force Field X. * * Force Field X is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License version 3 as published by * the Free Software Foundation. * * Force Field X is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more * details. * * You should have received a copy of the GNU General Public License along with * Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple * Place, Suite 330, Boston, MA 02111-1307 USA * * Linking this library statically or dynamically with other modules is making a * combined work based on this library. Thus, the terms and conditions of the * GNU General Public License cover the whole combination. * * As a special exception, the copyright holders of this library give you * permission to link this library with independent modules to produce an * executable, regardless of the license terms of these independent modules, and * to copy and distribute the resulting executable under terms of your choice, * provided that you also meet, for each linked independent module, the terms * and conditions of the license of that module. An independent module is a * module which is not derived from or based on this library. If you modify this * library, you may extend this exception to your version of the library, but * you are not obligated to do so. If you do not wish to do so, delete this * exception statement from your version. */ package ffx.numerics; import java.util.logging.Level; import java.util.logging.Logger; import static java.lang.Math.PI; import static java.lang.String.format; import static org.apache.commons.math3.util.FastMath.exp; import static org.apache.commons.math3.util.FastMath.pow; import static org.apache.commons.math3.util.FastMath.sqrt; import static ffx.numerics.Erf.erfc; import static ffx.numerics.VectorMath.binomial; import static ffx.numerics.VectorMath.diff; import static ffx.numerics.VectorMath.dot; import static ffx.numerics.VectorMath.doubleFactorial; import static ffx.numerics.VectorMath.norm; import static ffx.numerics.VectorMath.r; import static ffx.numerics.VectorMath.scalar; /** * The MultipoleTensor class computes derivatives of 1/|<b>r</b>| via recursion * to arbitrary order for Cartesian multipoles in either a global frame or a * quasi-internal frame. * * @author Michael J. Schnieders * * @see * <a href="http://doi.org/10.1142/9789812830364_0002" * target="_blank"> * Matt Challacombe, Eric Schwegler and Jan Almlof, Modern developments in * Hartree-Fock theory: Fast methods for computing the Coulomb matrix. * Computational Chemistry: Review of Current Trends. pp. 53-107, Ed. J. * Leczszynski, World Scientifc, 1996. * </a> * * @since 1.0 */ public class MultipoleTensor { private static final Logger logger = Logger.getLogger(MultipoleTensor.class.getName()); /** * Operators that are supported. */ public enum OPERATOR { COULOMB, SCREENED_COULOMB, THOLE_FIELD }; /** * Global and Quasi-Internal (QI) coordinate systems are supported. */ public enum COORDINATES { GLOBAL, QI }; private OPERATOR operator; private COORDINATES coordinates; private COORDINATES bufferCoordinates = COORDINATES.GLOBAL; private final int order; /** * These are the "source" terms for the recursion. Source terms exist for * for:<br> * 1) the Coulomb operator (1/R).<br> * 2) the screened Coulomb operator erfc(R)/R. * */ private double T000j[]; private final double coulomb[]; private final double screened[]; /** * Ewald parameter. */ private final double beta; /** * Thole damping parameters. */ private double damp; /** * 1/(alphai*alphak)^6 where alpha is polarizability. */ private double aiak; /** * Separation distance. */ private double R; /** * Separation distance squared. */ private double r2; /** * Xk - Xi. */ private double x; /** * Yk - Yi. */ private double y; /** * Yk - Yi. */ private double z; private final int o1; private final int il; private final int im; private final int in; private final int size; private final double sqrtPI = sqrt(PI); /** * Store the auxillary tensor memory to avoid memory consumption. */ private final double T000[]; /** * Store the work array to avoid memory consumption. Note that rather than * use an array for intermediate values, a 4D matrix was tried. It was * approximately 50% slower than the linear work array. */ private final double work[]; private double dEdF, d2EdF2; /** * <p> * Constructor for MultipoleTensor.</p> * * @param operator The tensor operator. * @param order The order of the tensor. * @param beta The screening parameter. */ public MultipoleTensor(OPERATOR operator, COORDINATES coordinates, int order, double beta) { assert (order > 0); o1 = order + 1; il = o1; im = il * o1; in = im * o1; size = (order + 1) * (order + 2) * (order + 3) / 6; work = new double[in * o1]; this.order = order; this.operator = operator; this.coordinates = coordinates; this.beta = beta; if (operator == OPERATOR.SCREENED_COULOMB && beta == 0.0) { // logger.warning("Tried beta of zero for screened coulomb tensor."); // Switch to the Coulomb operator. operator = OPERATOR.COULOMB; } // Auxillary terms for Coulomb and Thole Screening. coulomb = new double[o1]; for (int n = 0; n <= order; n++) { /** * Math.pow(-1.0, j) returns positive for all j, with -1.0 as the // * argument rather than -1. This is a bug? */ // Challacombe Eq. 21, first two factors. coulomb[n] = pow(-1, n) * doubleFactorial(2 * n - 1); } // Auxillary terms for screened Coulomb (Sagui et al. Eq. 2.28) screened = new double[o1]; double prefactor = 2.0 * beta / sqrtPI; double twoBeta2 = -2.0 * beta * beta; for (int n = 0; n <= order; n++) { screened[n] = prefactor * pow(twoBeta2, n); } setOperator(operator); T000 = new double[order + 1]; // l + m + n = 0 (1) t000 = ti(0, 0, 0, order); // l + m + n = 1 (3) 4 t100 = ti(1, 0, 0, order); t010 = ti(0, 1, 0, order); t001 = ti(0, 0, 1, order); // l + m + n = 2 (6) 10 t200 = ti(2, 0, 0, order); t020 = ti(0, 2, 0, order); t002 = ti(0, 0, 2, order); t110 = ti(1, 1, 0, order); t101 = ti(1, 0, 1, order); t011 = ti(0, 1, 1, order); // l + m + n = 3 (10) 20 t300 = ti(3, 0, 0, order); t030 = ti(0, 3, 0, order); t003 = ti(0, 0, 3, order); t210 = ti(2, 1, 0, order); t201 = ti(2, 0, 1, order); t120 = ti(1, 2, 0, order); t021 = ti(0, 2, 1, order); t102 = ti(1, 0, 2, order); t012 = ti(0, 1, 2, order); t111 = ti(1, 1, 1, order); // l + m + n = 4 (15) 35 t400 = ti(4, 0, 0, order); t040 = ti(0, 4, 0, order); t004 = ti(0, 0, 4, order); t310 = ti(3, 1, 0, order); t301 = ti(3, 0, 1, order); t130 = ti(1, 3, 0, order); t031 = ti(0, 3, 1, order); t103 = ti(1, 0, 3, order); t013 = ti(0, 1, 3, order); t220 = ti(2, 2, 0, order); t202 = ti(2, 0, 2, order); t022 = ti(0, 2, 2, order); t211 = ti(2, 1, 1, order); t121 = ti(1, 2, 1, order); t112 = ti(1, 1, 2, order); // l + m + n = 5 (21) 56 t500 = ti(5, 0, 0, order); t050 = ti(0, 5, 0, order); t005 = ti(0, 0, 5, order); t410 = ti(4, 1, 0, order); t401 = ti(4, 0, 1, order); t140 = ti(1, 4, 0, order); t041 = ti(0, 4, 1, order); t104 = ti(1, 0, 4, order); t014 = ti(0, 1, 4, order); t320 = ti(3, 2, 0, order); t302 = ti(3, 0, 2, order); t230 = ti(2, 3, 0, order); t032 = ti(0, 3, 2, order); t203 = ti(2, 0, 3, order); t023 = ti(0, 2, 3, order); t311 = ti(3, 1, 1, order); t131 = ti(1, 3, 1, order); t113 = ti(1, 1, 3, order); t221 = ti(2, 2, 1, order); t212 = ti(2, 1, 2, order); t122 = ti(1, 2, 2, order); // l + m + n = 6 (28) 84 t600 = ti(6, 0, 0, order); t060 = ti(0, 6, 0, order); t006 = ti(0, 0, 6, order); t510 = ti(5, 1, 0, order); t501 = ti(5, 0, 1, order); t150 = ti(1, 5, 0, order); t051 = ti(0, 5, 1, order); t105 = ti(1, 0, 5, order); t015 = ti(0, 1, 5, order); t420 = ti(4, 2, 0, order); t402 = ti(4, 0, 2, order); t240 = ti(2, 4, 0, order); t042 = ti(0, 4, 2, order); t204 = ti(2, 0, 4, order); t024 = ti(0, 2, 4, order); t411 = ti(4, 1, 1, order); t141 = ti(1, 4, 1, order); t114 = ti(1, 1, 4, order); t330 = ti(3, 3, 0, order); t303 = ti(3, 0, 3, order); t033 = ti(0, 3, 3, order); t321 = ti(3, 2, 1, order); t231 = ti(2, 3, 1, order); t213 = ti(2, 1, 3, order); t312 = ti(3, 1, 2, order); t132 = ti(1, 3, 2, order); t123 = ti(1, 2, 3, order); t222 = ti(2, 2, 2, order); } /** * Set the Operator. * * @param operator */ public final void setOperator(OPERATOR operator) { this.operator = operator; OPERATOR op = operator; if (operator == OPERATOR.SCREENED_COULOMB && beta == 0.0) { // logger.warning("Tried beta of zero for screened coulomb tensor."); // Switch to the Coulomb operator. op = OPERATOR.COULOMB; } switch (op) { case SCREENED_COULOMB: T000j = screened; break; default: case THOLE_FIELD: case COULOMB: T000j = coulomb; } } public void generateTensor() { switch (order) { case 5: generateTensor5(); break; case 4: generateTensor4(); break; default: double r[] = { x, y, z }; recursion(r, work); } } public void generateTensor4() { switch (coordinates) { case QI: order4QI(); break; case GLOBAL: default: order4(); } } public void generateTensor5() { switch (coordinates) { case QI: order5QI(); break; case GLOBAL: default: order5(); } } public void setCoordinateSystem(COORDINATES coordinates) { this.coordinates = coordinates; } /** * Set the Thole damping parameters. * * @param damp * @param aiak 1/(alphai*alphak)^6 where alpha is polarizability */ public void setTholeDamping(double damp, double aiak) { this.damp = damp; // == PME's pgamma this.aiak = aiak; // == 1/(alphai*alphak)^6 where alpha is polarizability } public boolean applyDamping() { double rdamp = R * aiak; double test = -damp * rdamp * rdamp * rdamp; if (test > -50.0) { return true; } return false; } public void setR(double r[]) { switch (coordinates) { case QI: x = 0.0; y = 0.0; z = r(r); setQIRotationMatrix(r[0], r[1], r[2]); break; default: case GLOBAL: x = r[0]; y = r[1]; z = r[2]; r2 = (x * x + y * y + z * z); if (r2 == 0.0) { throw new ArithmeticException(); } R = sqrt(r2); } } public void setR(double r[], double lambdaFunction) { switch (coordinates) { case QI: setR_QI(r, lambdaFunction); break; case GLOBAL: x = r[0]; y = r[1]; z = r[2] + lambdaFunction; r2 = (x * x + y * y + z * z); if (r2 == 0.0) { throw new ArithmeticException(); } R = sqrt(r2); } } public void setR_QI(double r[]) { x = 0.0; y = 0.0; r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2]; if (r2 == 0.0) { throw new ArithmeticException(); } z = sqrt(r2); R = z; setQIRotationMatrix(r[0], r[1], r[2]); } public void setR_QI(double r[], double lambdaFunction) { x = 0.0; y = 0.0; double zl = r[2] + lambdaFunction; r2 = r[0] * r[0] + r[1] * r[1] + zl * zl; if (r2 == 0.0) { throw new ArithmeticException(); } z = sqrt(r2); R = z; setQIRotationMatrix(r[0], r[1], r[2] + lambdaFunction); } public double multipoleEnergy(double Fi[], double Ti[], double Tk[]) { switch (coordinates) { case GLOBAL: default: return multipoleEnergyGlobal(Fi, Ti, Tk); case QI: return multipoleEnergyQI(Fi, Ti, Tk); } } public double polarizationEnergy(double scaleField, double scaleEnergy, double scaleMutual, double Fi[], double Ti[], double Tk[]) { switch (coordinates) { case GLOBAL: default: return polarizationEnergyGlobal(scaleField, scaleEnergy, scaleMutual, Fi, Ti, Tk); case QI: return polarizationEnergyQI(scaleField, scaleEnergy, scaleMutual, Fi, Ti, Tk); } } /** * Generate source terms for the Challacombe et al. recursion. * * @param r Cartesian vector. * @param T000 Location to store the source terms. * @param damp Thole damping for this interaction. * @param aiak Inverse of the polarizability product to the sixth, ie * 1/(site i * site k)^6. */ private void source(double T000[]) { OPERATOR op = operator; if (operator == OPERATOR.SCREENED_COULOMB && beta == 0.0) { // Generate tensors for the Coulomb operator. op = OPERATOR.COULOMB; } switch (op) { case SCREENED_COULOMB: // Sagui et al. Eq. 2.22 // "beta" here == "aewald" from PME which is *NOT* "beta" from PME // (What PME calls "beta" (now "lambdaBufferDist") is the lambda buffer, which stays out of this class.) double betaR = beta * R; double betaR2 = betaR * betaR; double iBetaR2 = 1.0 / (2.0 * betaR2); double expBR2 = exp(-betaR2); // Fnc(x^2) = Sqrt(PI) * erfc(x) / (2*x) // where x = Beta*R double Fnc = sqrtPI * erfc(betaR) / (2.0 * betaR); for (int n = 0; n < o1; n++) { T000[n] = T000j[n] * Fnc; // Generate F(n+1)c from Fnc (Eq. 2.24 in Sagui et al.) // F(n+1)c = [(2*n+1) Fnc(x) + exp(-x)] / 2x // where x = (Beta*R)^2 Fnc = ((2.0 * n + 1.0) * Fnc + expBR2) * iBetaR2; } break; case THOLE_FIELD: assert (order <= 4); double ir = 1.0 / R; double ir2 = ir * ir; for (int n = 0; n < o1; n++) { T000[n] = T000j[n] * ir; ir *= ir2; } /** * Add the Thole damping terms: edamp = exp(-damp*u^3). */ double u = R * aiak; double u3 = damp * u * u * u; double u6 = u3 * u3; double u9 = u6 * u3; double expU3 = exp(-u3); T000[0] = 0.0; // The zeroth order term is not calculated for Thole damping. T000[1] *= expU3; T000[2] *= (1.0 + u3) * expU3; T000[3] *= (1.0 + u3 + threeFifths * u6) * expU3; T000[4] *= (1.0 + u3 + (18.0 * u6 + 9.0 * u9) * oneThirtyFifth) * expU3; break; case COULOMB: default: // Challacombe et al. Equation 21, last factor. // == (1/r) * (1/r^3) * (1/r^5) * (1/r^7) * ... ir = 1.0 / R; ir2 = ir * ir; StringBuilder sb = new StringBuilder(); sb.append(format("Coulomb shit for R,ir = %g,%g: 000(j),001(j),...\n", R, ir)); for (int n = 0; n < o1; n++) { T000[n] = T000j[n] * ir; sb.append(format(" %g (%g)\n", T000[n], T000j[n])); ir *= ir2; } // logger.info(sb.toString()); } } private static final double threeFifths = 3.0 / 5.0; private static final double oneThirtyFifth = 1.0 / 35.0; /** * Log the tensors. * * @param tensor */ public void log(double tensor[]) { StringBuilder sb = new StringBuilder(); sb.append(String.format("\n %s Operator to order %d:", operator, order)); sb.append(String.format("\n%5s %4s %4s %4s %12s\n", "Index", "d/dx", "d/dy", "d/dz", "Tensor")); sb.append(String.format("%5d %4d %4d %4d %12.8f\n", 0, 0, 0, 0, tensor[0])); int count = 1; // Print (d/dx)^l for l = 1..order (m = 0, n = 0) for (int l = 1; l <= order; l++) { double value = tensor[ti(l, 0, 0)]; if (value != 0.0) { sb.append(String.format("%5d %4d %4d %4d %12.8f\n", ti(l, 0, 0), l, 0, 0, value)); count++; } } // Print (d/dx)^l * (d/dy)^m for l + m = 1..order (m >= 1, n = 0) for (int l = 0; l <= o1; l++) { for (int m = 1; m <= order - l; m++) { double value = tensor[ti(l, m, 0)]; if (value != 0.0) { sb.append(String.format("%5d %4d %4d %4d %12.8f\n", ti(l, m, 0), l, m, 0, value)); count++; } } } // Print (d/dx)^l * (d/dy)^m * (d/dz)^n for l + m + n = 1..o (n >= 1) for (int l = 0; l <= o1; l++) { for (int m = 0; m <= o1 - l; m++) { for (int n = 1; n <= order - l - m; n++) { double value = tensor[ti(l, m, n)]; if (value != 0.0) { sb.append(String.format("%5d %4d %4d %4d %12.8f\n", ti(l, m, n), l, m, n, value)); count++; } } } } sb.append(String.format("\n Total number of active tensors: %d\n", count)); logger.log(Level.INFO, sb.toString()); } /** * Returns the number of tensors for derivatives to the given order. * * @param order maximum number of derivatives. * @return the number of tensors. * @since 1.0 */ public static int tensorCount(int order) { long ret = binomial(order + 3, 3); assert (ret < Integer.MAX_VALUE); return (int) ret; } /** * The index is based on the idea of filling tetrahedron. * <p> * 1/r has an index of 0 * <br> * derivatives of x are first; indeces from 1..o for d/dx..(d/dx)^o * <br> * derivatives of x and y are second; base triangle of size (o+1)(o+2)/2 * <br> * derivatives of x, y and z are last; total size (o+1)*(o+2)*(o+3)/6 * <br> * <p> * This function is useful to set up masking constants: * <br> * static int Tlmn = ti(l,m,n,order) * <br> * For example the (d/dy)^2 (1/R) storage location: <br> static int T020 = * ti(0,2,0,order) * <p> * * @param dx int The number of d/dx operations. * @param dy int The number of d/dy operations. * @param dz int The number of d/dz operations. * @param order int The maximum tensor order (0 .LE. dx + dy + dz .LE. * order). * * @return int in the range (0..binomial(order + 3, 3) - 1) */ public static int ti(int dx, int dy, int dz, int order) { if (dx < 0 || dy < 0 || dz < 0 || dx + dy + dz > order) { return -1; } int size = (order + 1) * (order + 2) * (order + 3) / 6; /** * We only get to the top of the tetrahedron if dz = order, otherwise * subtract off the top, including the level of the requested tensor * index. */ int top = order + 1 - dz; top = top * (top + 1) * (top + 2) / 6; int zindex = size - top; /** * Given the "dz level", dy can range from 0..order - dz) To get to the * row for a specific value of dy, dy*(order + 1) - dy*(dy-1)/2 indeces * are skipped. This is an operation that looks like the area of * rectangle, minus the area of an empty triangle. */ int yindex = dy * (order - dz) - (dy - 1) * (dy - 2) / 2 + 1; /** * Given the dz level and dy row, dx can range from (0..order - dz - dy) * The dx index is just walking down the dy row for "dx" steps. */ int ret = dx + yindex + zindex; return ret; } /** * The index is based on the idea of filling tetrahedron. * <p> * 1/r has an index of 0 * <br> * derivatives of x are first; indeces from 1..o for d/dx..(d/dx)^o * <br> * derivatives of x and y are second; base triangle of size (o+1)(o+2)/2 * <br> * derivatives of x, y and z are last; total size (o+1)*(o+2)*(o+3)/6 * <br> * <p> * This function is useful to set up masking constants: * <br> * static int Tlmn = ti(l,m,n,order) * <br> * For example the (d/dy)^2 (1/R) storage location: <br> static int T020 = * ti(0,2,0,order) * <p> * * @param dx int The number of d/dx operations. * @param dy int The number of d/dy operations. * @param dz int The number of d/dz operations. * * @return int in the range (0..binomial(order + 3, 3) - 1) */ public int ti(int dx, int dy, int dz) { return ti(dx, dy, dz, order); } /** * This method is a driver to collect elements of the Cartesian multipole * tensor given the recursion relationships implemented by the method * "Tlmnj", which can be called directly to get a single tensor element. It * does not store intermediate values of the recursion, causing it to scale * O(order^8). For order = 5, this approach is a factor of 10 slower than * recursion. * * @param r double[] vector between two sites. * @param tensor double[] length must be at least binomial(order + 3, 3). */ public void noStorageRecursion(double r[], double tensor[]) { setR(r); source(T000); // 1/r tensor[0] = T000[0]; // Find (d/dx)^l for l = 1..order (m = 0, n = 0) for (int l = 1; l <= order; l++) { tensor[ti(l, 0, 0)] = Tlmnj(l, 0, 0, 0, r, T000); } // Find (d/dx)^l * (d/dy)^m for l + m = 1..order (m >= 1, n = 0) for (int l = 0; l <= o1; l++) { for (int m = 1; m <= order - l; m++) { tensor[ti(l, m, 0)] = Tlmnj(l, m, 0, 0, r, T000); } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l + m + n = 1..o (n >= 1) for (int l = 0; l <= o1; l++) { for (int m = 0; m <= o1 - l; m++) { for (int n = 1; n <= order - l - m; n++) { tensor[ti(l, m, n)] = Tlmnj(l, m, n, 0, r, T000); } } } } /** * This method is a driver to collect elements of the Cartesian multipole * tensor given the recursion relationships implemented by the method * "Tlmnj", which can be called directly to get a single tensor element. It * does not store intermediate values of the recursion, causing it to scale * O(order^8). For order = 5, this approach is a factor of 10 slower than * recursion. * * @param r double[] vector between two sites. r[0] and r[1] must equal 0.0. * @param tensor double[] length must be at least binomial(order + 3, 3). */ public void noStorageRecursionQI(double r[], double tensor[]) { assert (r[0] == 0.0 && r[1] == 0.0); setR_QI(r); source(T000); // 1/r tensor[0] = T000[0]; // Find (d/dx)^l for l = 1..order (m = 0, n = 0) for (int l = 1; l <= order; l++) { tensor[ti(l, 0, 0)] = TlmnjQI(l, 0, 0, 0, r, T000); } // Find (d/dx)^l * (d/dy)^m for l + m = 1..order (m >= 1, n = 0) for (int l = 0; l <= o1; l++) { for (int m = 1; m <= order - l; m++) { tensor[ti(l, m, 0)] = TlmnjQI(l, m, 0, 0, r, T000); } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l + m + n = 1..o (n >= 1) for (int l = 0; l <= o1; l++) { for (int m = 0; m <= o1 - l; m++) { for (int n = 1; n <= order - l - m; n++) { tensor[ti(l, m, n)] = TlmnjQI(l, m, n, 0, r, T000); } } } } /** * This routine implements the recurrence relations for computation of any * Cartesian multipole tensor in ~O(L^8) time, where L is the total order l * + m + n, given the auxiliary elements T0000. * <br> * It implements the recursion relationships in brute force fashion, without * saving intermediate values. This is useful for finding a single tensor, * rather than all binomial(L + 3, 3). * <br> * The specific recursion equations (41-43) and set of auxiliary tensor * elements from equation (40) can be found in Challacombe et al. * * @param l int The number of (d/dx) operations. * @param m int The number of (d/dy) operations. * @param n int The number of (d/dz) operations. * @param j int j = 0 is the Tlmn tensor, j .GT. 0 is an intermediate. * @param r double[] The {x,y,z} coordinates. * @param T000 double[] Initial auxiliary tensor elements from Eq. (40). * * @return double The requested Tensor element (intermediate if j .GT. 0). * * @since 1.0 */ public static double Tlmnj(final int l, final int m, final int n, final int j, final double[] r, final double[] T000) { if (m == 0 && n == 0) { if (l > 1) { return r[0] * Tlmnj(l - 1, 0, 0, j + 1, r, T000) + (l - 1) * Tlmnj(l - 2, 0, 0, j + 1, r, T000); } else if (l == 1) { // l == 1, d/dx is done. return r[0] * Tlmnj(0, 0, 0, j + 1, r, T000); } else { // l = m = n = 0. Recursion is done. return T000[j]; } } else if (n == 0) { // m >= 1 if (m > 1) { return r[1] * Tlmnj(l, m - 1, 0, j + 1, r, T000) + (m - 1) * Tlmnj(l, m - 2, 0, j + 1, r, T000); } return r[1] * Tlmnj(l, 0, 0, j + 1, r, T000); } else { // n >= 1 if (n > 1) { return r[2] * Tlmnj(l, m, n - 1, j + 1, r, T000) + (n - 1) * Tlmnj(l, m, n - 2, j + 1, r, T000); } return r[2] * Tlmnj(l, m, 0, j + 1, r, T000); } } /** * This routine implements the recurrence relations for computation of any * Cartesian multipole tensor in ~O(L^8) time, where L is the total order l * + m + n, given the auxiliary elements T0000. * <br> * It implements the recursion relationships in brute force fashion, without * saving intermediate values. This is useful for finding a single tensor, * rather than all binomial(L + 3, 3). * <br> * The specific recursion equations (41-43) and set of auxiliary tensor * elements from equation (40) can be found in Challacombe et al. * * @param l int The number of (d/dx) operations. * @param m int The number of (d/dy) operations. * @param n int The number of (d/dz) operations. * @param j int j = 0 is the Tlmn tensor, j .GT. 0 is an intermediate. * @param r double[] The {x,y,z} coordinates. * @param T000 double[] Initial auxiliary tensor elements from Eq. (40). * * @return double The requested Tensor element (intermediate if j .GT. 0). * * @since 1.0 */ public static double TlmnjQI(final int l, final int m, final int n, final int j, final double[] r, final double[] T000) { double z = r[2]; assert (r[0] == 0.0 && r[1] == 0.0); if (m == 0 && n == 0) { if (l > 1) { return (l - 1) * Tlmnj(l - 2, 0, 0, j + 1, r, T000); } else if (l == 1) { // l == 1, d/dx is done. return 0.0; } else { // l = m = n = 0. Recursion is done. return T000[j]; } } else if (n == 0) { // m >= 1 if (m > 1) { return (m - 1) * Tlmnj(l, m - 2, 0, j + 1, r, T000); } return 0.0; } else { // n >= 1 if (n > 1) { return z * Tlmnj(l, m, n - 1, j + 1, r, T000) + (n - 1) * Tlmnj(l, m, n - 2, j + 1, r, T000); } return z * Tlmnj(l, m, 0, j + 1, r, T000); } } /** * This function is a driver to collect elements of the Cartesian multipole * tensor. Collecting all tensors scales slightly better than O(order^4). * <p> * For a multipole expansion truncated at quadrupole order, for example, up * to order 5 is needed for energy gradients. The number of terms this * requires is binomial(5 + 3, 3) or 8! / (5! * 3!), which is 56. * <p> * The packing of the tensor elements for order = 1 * <br> * tensor[0] = 1/|r| <br> * tensor[1] = -x/|r|^3 <br> * tensor[2] = -y/|r|^3 <br> * tensor[3] = -z/|r|^3 <br> * <p> * @param r double[] vector between two sites. * @param tensor double[] length must be at least binomial(order + 3, 3). * @since 1.0 */ public void recursion(final double r[], final double tensor[]) { setR(r); source(work); tensor[0] = work[0]; // Find (d/dx)^l for l = 1..order (m = 0, n = 0) // Any (d/dx) term can be formed as // Tl00j = x * T(l-1)00(j+1) + (l-1) * T(l-2)00(j+1) // All intermediate terms are indexed as l*il + m*im + n*in + j; double current; double previous = work[1]; // Store the l=1 tensor T100 (d/dx) tensor[ti(1, 0, 0)] = x * previous; // Starting the loop at l=2 avoids an if statement. for (int l = 2; l < o1; l++) { // Initial condition for the inner loop is formation of T100(l-1). // Starting the inner loop at a=2 avoid an if statement. // T100(l-1) = x * T000(l) current = x * work[l]; int iw = il + l - 1; work[iw] = current; for (int a = 1; a < l - 1; a++) { // T200(l-2) = x * T100(l-1) + (2 - 1) * T000(l-1) // T300(l-3) = x * T200(l-2) + (3 - 1) * T100(l-2) // ... // T(l-1)001 = x * T(l-2)002 + (l - 2) * T(l-3)002 current = x * current + a * work[iw - il]; iw += il - 1; work[iw] = current; } // Store the Tl00 tensor (d/dx)^l // Tl00 = x * [[ T(l-1)001 ]] + (l - 1) * T(l-2)001 tensor[ti(l, 0, 0)] = x * current + (l - 1) * previous; previous = current; } // Find (d/dx)^l * (d/dy)^m for l+m = 1..order (m > 0, n = 0) // Any (d/dy) term can be formed as: // Tlm0j = y * Tl(m-1)00(j+1) + (m-1) * Tl(m-2)00(j+1) for (int l = 0; l < order; l++) { // Store the m=1 tensor (d/dx)^l *(d/dy) // Tl10 = y * Tl001 previous = work[l * il + 1]; tensor[ti(l, 1, 0)] = y * previous; for (int m = 2; m + l < o1; m++) { // Tl10(m-1) = y * Tl00m; int iw = l * il + m; current = y * work[iw]; iw += im - 1; work[iw] = current; for (int a = 1; a < m - 1; a++) { // Tl20(m-2) = Y * Tl10(m-1) + (2 - 1) * T100(m-1) // Tl30(m-3) = Y * Tl20(m-2) + (3 - 1) * Tl10(m-2) // ... // Tl(m-1)01 = Y * Tl(m-2)02 + (m - 2) * T(m-3)02 current = y * current + a * work[iw - im]; iw += im - 1; work[iw] = current; } // Store the tensor (d/dx)^l * (d/dy)^m // Tlm0 = y * Tl(m-1)01 + (m - 1) * Tl(m-2)01 tensor[ti(l, m, 0)] = y * current + (m - 1) * previous; previous = current; } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l+m+n = 1..order (n > 0) // Any (d/dz) term can be formed as: // Tlmnj = z * Tlm(n-1)(j+1) + (n-1) * Tlm(n-2)(j+1) for (int l = 0; l < order; l++) { for (int m = 0; m + l < order; m++) { // Store the n=1 tensor (d/dx)^l *(d/dy)^m * (d/dz) // Tlmn = z * Tlm01 final int lm = m + l; final int lilmim = l * il + m * im; previous = work[lilmim + 1]; tensor[ti(l, m, 1)] = z * previous; for (int n = 2; lm + n < o1; n++) { // Tlm1(n-1) = z * Tlm0n; int iw = lilmim + n; current = z * work[iw]; iw += in - 1; work[iw] = current; final int n1 = n - 1; for (int a = 1; a < n1; a++) { // Tlm2(n-2) = z * Tlm1(n-1) + (2 - 1) * T1m0(n-1) // Tlm3(n-3) = z * Tlm2(n-2) + (3 - 1) * Tlm1(n-2) // ... // Tlm(n-1)1 = z * Tlm(n-2)2 + (n - 2) * Tlm(n-3)2 current = z * current + a * work[iw - in]; iw += in - 1; work[iw] = current; } // Store the tensor (d/dx)^l * (d/dy)^m * (d/dz)^n // Tlmn = z * Tlm(n-1)1 + (n - 1) * Tlm(n-2)1 tensor[ti(l, m, n)] = z * current + n1 * previous; previous = current; } } } } /** * This function is a driver to collect elements of the Cartesian multipole * tensor for Quasi-Internal coordinate evaluation. * * Thus, we assume r[0] and r[1] (i.e. dx and dy) are zero. * * Collecting all tensors scales better than O(order^4). * <p> * For a multipole expansion truncated at quadrupole order, for example, up * to order 5 is needed for energy gradients. The number of terms this * requires is binomial(5 + 3, 3) or 8! / (5! * 3!), which is 56. * <p> * The packing of the tensor elements for order = 1<br> tensor[0] = 1/|r| * <br> * tensor[1] = -x/|r|^3 <br> * tensor[2] = -y/|r|^3 <br> * tensor[3] = -z/|r|^3 <br> * <p> * * @param r double[] vector between two sites (assumes r[0] and r[1] = 0.0). * @param tensor double[] length must be at least binomial(order + 3, 3). * @since 1.0 */ public void recursionQI(final double r[], final double tensor[]) { setR_QI(r); assert (x == 0.0 && y == 0.0); source(work); tensor[0] = work[0]; // Find (d/dx)^l for l = 1..order (m = 0, n = 0) // Any (d/dx) term can be formed as // Tl00j = x * T(l-1)00(j+1) + (l-1) * T(l-2)00(j+1) // All intermediate terms are indexed as l*il + m*im + n*in + j; // Store the l=1 tensor T100 (d/dx) // Starting the loop at l=2 avoids an if statement. double current; double previous = work[1]; tensor[ti(1, 0, 0)] = 0.0; for (int l = 2; l < o1; l++) { // Initial condition for the inner loop is formation of T100(l-1). // Starting the inner loop at a=2 avoid an if statement. // T100(l-1) = 0.0 * T000(l) current = 0.0; int iw = il + l - 1; work[iw] = current; for (int a = 1; a < l - 1; a++) { // T200(l-2) = 0.0 * T100(l-1) + (2 - 1) * T000(l-1) // T300(l-3) = 0.0 * T200(l-2) + (3 - 1) * T100(l-2) // ... // T(l-1)001 = 0.0 * T(l-2)002 + (l - 2) * T(l-3)002 current = a * work[iw - il]; iw += il - 1; work[iw] = current; } // Store the Tl00 tensor (d/dx)^l // Tl00 = 0.0 * [[ T(l-1)001 ]] + (l - 1) * T(l-2)001 tensor[ti(l, 0, 0)] = (l - 1) * previous; previous = current; } // Find (d/dx)^l * (d/dy)^m for l+m = 1..order (m > 0, n = 0) // Any (d/dy) term can be formed as: // Tlm0j = y * Tl(m-1)00(j+1) + (m-1) * Tl(m-2)00(j+1) for (int l = 0; l < order; l++) { // Store the m=1 tensor (d/dx)^l *(d/dy) // Tl10 = y * Tl001 previous = work[l * il + 1]; tensor[ti(l, 1, 0)] = 0.0; for (int m = 2; m + l < o1; m++) { // Tl10(m-1) = y * Tl00m; int iw = l * il + m; current = 0.0; iw += im - 1; work[iw] = current; for (int a = 1; a < m - 1; a++) { // Tl20(m-2) = 0.0 * Tl10(m-1) + (2 - 1) * T100(m-1) // Tl30(m-3) = 0.0 * Tl20(m-2) + (3 - 1) * Tl10(m-2) // ... // Tl(m-1)01 = 0.0 * Tl(m-2)02 + (m - 2) * Tl(m-3)02 current = a * work[iw - im]; iw += im - 1; work[iw] = current; } // Store the tensor (d/dx)^l * (d/dy)^m // Tlm0 = y * Tl(m-1)01 + (m - 1) * Tl(m-2)01 tensor[ti(l, m, 0)] = (m - 1) * previous; previous = current; } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l+m+n = 1..order (n > 0) // Any (d/dz) term can be formed as: // Tlmnj = z * Tlm(n-1)(j+1) + (n-1) * Tlm(n-2)(j+1) for (int l = 0; l < order; l++) { for (int m = 0; m + l < order; m++) { // Store the n=1 tensor (d/dx)^l *(d/dy)^m * (d/dz) // Tlmn = z * Tlm01 final int lm = m + l; final int lilmim = l * il + m * im; previous = work[lilmim + 1]; tensor[ti(l, m, 1)] = z * previous; for (int n = 2; lm + n < o1; n++) { // Tlm1(n-1) = z * Tlm0n; int iw = lilmim + n; current = z * work[iw]; iw += in - 1; work[iw] = current; final int n1 = n - 1; for (int a = 1; a < n1; a++) { // Tlm2(n-2) = z * Tlm1(n-1) + (2 - 1) * T1m0(n-1) // Tlm3(n-3) = z * Tlm2(n-2) + (3 - 1) * Tlm1(n-2) // ... // Tlm(n-1)1 = z * Tlm(n-2)2 + (n - 2) * Tlm(n-3)2 current = z * current + a * work[iw - in]; iw += in - 1; work[iw] = current; } // Store the tensor (d/dx)^l * (d/dy)^m * (d/dz)^n // Tlmn = z * Tlm(n-1)1 + (n - 1) * Tlm(n-2)1 tensor[ti(l, m, n)] = z * current + n1 * previous; previous = current; } } } } /** * This function is a driver to collect elements of the Cartesian multipole * tensor. Collecting all tensors scales slightly better than O(order^4). * <p> * For a multipole expansion truncated at quadrupole order, for example, up * to order 5 is needed for energy gradients. The number of terms this * requires is binomial(5 + 3, 3) or 8! / (5! * 3!), which is 56. * <p> * The packing of the tensor elements for order = 1<br> tensor[0] = 1/|r| * <br> * tensor[1] = -x/|r|^3 <br> tensor[2] = -y/|r|^3 <br> tensor[3] = -z/|r|^3 * <br> * <p> * * @param r double[] vector between two sites. * @param tensor double[] length must be at least binomial(order + 3, 3). * @return Java code for the tensor recursion. * * @since 1.0 */ public String codeTensorRecursion(final double r[], final double tensor[]) { setR(r); source(work); StringBuilder sb = new StringBuilder(); tensor[0] = work[0]; if (work[0] > 0) { sb.append(format("%s = %s;\n", rlmn(0, 0, 0), term(0, 0, 0, 0))); } // Find (d/dx)^l for l = 1..order (m = 0, n = 0) // Any (d/dx) term can be formed as // Tl00j = x * T(l-1)00(j+1) + (l-1) * T(l-2)00(j+1) // All intermediate terms are indexed as l*il + m*im + n*in + j; double current; double previous = work[1]; // Store the l=1 tensor T100 (d/dx) tensor[ti(1, 0, 0)] = x * previous; sb.append(format("%s = x * %s;\n", rlmn(1, 0, 0), term(0, 0, 0, 1))); // Starting the loop at l=2 avoids an if statement. for (int l = 2; l < o1; l++) { // Initial condition for the inner loop is formation of T100(l-1). // Starting the inner loop at a=2 avoid an if statement. // T100(l-1) = x * T000(l) current = x * work[l]; int iw = il + l - 1; work[iw] = current; sb.append(format("double %s = x * %s;\n", term(1, 0, 0, l - 1), term(0, 0, 0, l))); for (int a = 1; a < l - 1; a++) { // T200(l-2) = x * T100(l-1) + (2 - 1) * T000(l-1) // T300(l-3) = x * T200(l-2) + (3 - 1) * T100(l-2) // ... // T(l-1)001 = x * T(l-2)002 + (l - 2) * T(l-3)002 current = x * current + a * work[iw - il]; iw += il - 1; work[iw] = current; if (a > 1) { sb.append(format("double %s = x * %s + %d * %s;\n", term(a + 1, 0, 0, l - a - 1), term(a, 0, 0, l - a), a, term(a - 1, 0, 0, l - a))); } else { sb.append(format("double %s = x * %s + %s;\n", term(a + 1, 0, 0, l - a - 1), term(a, 0, 0, l - a), term(a - 1, 0, 0, l - a))); } } // Store the Tl00 tensor (d/dx)^l // Tl00 = x * [[ T(l-1)001 ]] + (l - 1) * T(l-2)001 tensor[ti(l, 0, 0)] = x * current + (l - 1) * previous; previous = current; if (l > 2) { sb.append(format("%s = x * %s + %d * %s;\n", rlmn(l, 0, 0), term(l - 1, 0, 0, 1), (l - 1), term(l - 2, 0, 0, 1))); } else { sb.append(format("%s = x * %s + %s;\n", rlmn(l, 0, 0), term(l - 1, 0, 0, 1), term(l - 2, 0, 0, 1), l, 0, 0)); } } // Find (d/dx)^l * (d/dy)^m for l+m = 1..order (m > 0, n = 0) // Any (d/dy) term can be formed as: // Tlm0j = y * Tl(m-1)00(j+1) + (m-1) * Tl(m-2)00(j+1) for (int l = 0; l < order; l++) { // Store the m=1 tensor (d/dx)^l *(d/dy) // Tl10 = y * Tl001 previous = work[l * il + 1]; tensor[ti(l, 1, 0)] = y * previous; sb.append(format("%s = y * %s;\n", rlmn(l, 1, 0), term(l, 0, 0, 1))); for (int m = 2; m + l < o1; m++) { // Tl10(m-1) = y * Tl00m; int iw = l * il + m; current = y * work[iw]; iw += im - 1; work[iw] = current; sb.append(format("double %s = y * %s;\n", term(l, 1, 0, m - 1), term(l, 0, 0, m))); for (int a = 1; a < m - 1; a++) { // Tl20(m-2) = Y * Tl10(m-1) + (2 - 1) * T100(m-1) // Tl30(m-3) = Y * Tl20(m-2) + (3 - 1) * Tl10(m-2) // ... // Tl(m-1)01 = Y * Tl(m-2)02 + (m - 2) * T(m-3)02 current = y * current + a * work[iw - im]; iw += im - 1; work[iw] = current; if (a > 1) { sb.append(format("double %s = y * %s + %d * %s;\n", term(l, a + 1, 0, m - a - 1), term(l, a, 0, m - a), a, term(l, a - 1, 0, m - a))); } else { sb.append(format("double %s = y * %s + %s;\n", term(l, a + 1, 0, m - a - 1), term(l, a, 0, m - a), term(l, a - 1, 0, m - a))); } } // Store the tensor (d/dx)^l * (d/dy)^m // Tlm0 = y * Tl(m-1)01 + (m - 1) * Tl(m-2)01 tensor[ti(l, m, 0)] = y * current + (m - 1) * previous; previous = current; if (m > 2) { sb.append(format("%s = y * %s + %d * %s;\n", rlmn(l, m, 0), term(l, m - 1, 0, 1), (m - 1), term(l, m - 2, 0, 1))); } else { sb.append(format("%s = y * %s + %s;\n", rlmn(l, m, 0), term(l, m - 1, 0, 1), term(l, m - 2, 0, 1))); } } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l+m+n = 1..order (n > 0) // Any (d/dz) term can be formed as: // Tlmnj = z * Tlm(n-1)(j+1) + (n-1) * Tlm(n-2)(j+1) for (int l = 0; l < order; l++) { for (int m = 0; m + l < order; m++) { // Store the n=1 tensor (d/dx)^l *(d/dy)^m * (d/dz) // Tlmn = z * Tlm01 final int lm = m + l; final int lilmim = l * il + m * im; previous = work[lilmim + 1]; tensor[ti(l, m, 1)] = z * previous; sb.append(format("%s = z * %s;\n", rlmn(l, m, 1), term(l, m, 0, 1), l, m, 1)); for (int n = 2; lm + n < o1; n++) { // Tlm1(n-1) = z * Tlm0n; int iw = lilmim + n; current = z * work[iw]; iw += in - 1; work[iw] = current; sb.append(format("double %s = z * %s;\n", term(l, m, 1, n - 1), term(l, m, 0, n))); final int n1 = n - 1; for (int a = 1; a < n1; a++) { // Tlm2(n-2) = z * Tlm1(n-1) + (2 - 1) * T1m0(n-1) // Tlm3(n-3) = z * Tlm2(n-2) + (3 - 1) * Tlm1(n-2) // ... // Tlm(n-1)1 = z * Tlm(n-2)2 + (n - 2) * Tlm(n-3)2 current = z * current + a * work[iw - in]; iw += in - 1; work[iw] = current; if (a > 1) { sb.append(format("double %s = z * %s + %d * %s;\n", term(l, m, a + 1, n - a - 1), term(l, m, a, n - a), a, term(l, m, a - 1, n - a))); } else { sb.append(format("double %s = z * %s + %s;\n", term(l, m, a + 1, n - a - 1), term(l, m, a, n - a), term(l, m, a - 1, n - a))); } } // Store the tensor (d/dx)^l * (d/dy)^m * (d/dz)^n // Tlmn = z * Tlm(n-1)1 + (n - 1) * Tlm(n-2)1 tensor[ti(l, m, n)] = z * current + n1 * previous; previous = current; if (n > 2) { sb.append(format("%s = z * %s + %d * %s;\n", rlmn(l, m, n), term(l, m, n - 1, 1), (n - 1), term(l, m, n - 2, 1), l, m, n)); } else { sb.append(format("%s = z * %s + %s;\n", rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1), l, m, n)); } } } } return sb.toString(); } /** * This function write Java code to collect elements of the Cartesian * multipole tensor for Quasi-Internal coordinate evaluation. * * The r[0] and r[1] (i.e. dx and dy) must be zero. * * Collecting all tensors scales better than O(order^4). * <p> * For a multipole expansion truncated at quadrupole order, for example, up * to order 5 is needed for energy gradients. The number of terms this * requires is binomial(5 + 3, 3) or 8! / (5! * 3!), which is 56. * <p> * The packing of the tensor elements for order = 1<br> tensor[0] = 1/|r| * <br> * tensor[1] = -x/|r|^3 <br> * tensor[2] = -y/|r|^3 <br> * tensor[3] = -z/|r|^3 <br> * <p> * * @param r double[] vector between two sites (assumes r[0] and r[1] = 0.0). * @param tensor double[] length must be at least binomial(order + 3, 3). * @return Java code to build the tensors. * * @since 1.0 */ public String codeTensorRecursionQI(final double r[], final double tensor[]) { setR_QI(r); assert (x == 0.0 && y == 0.0); source(work); StringBuilder sb = new StringBuilder(); tensor[0] = work[0]; if (work[0] > 0) { sb.append(format("%s = %s;\n", rlmn(0, 0, 0), term(0, 0, 0, 0))); } // Find (d/dx)^l for l = 1..order (m = 0, n = 0) // Any (d/dx) term can be formed as // Tl00j = x * T(l-1)00(j+1) + (l-1) * T(l-2)00(j+1) // All intermediate terms are indexed as l*il + m*im + n*in + j; // Store the l=1 tensor T100 (d/dx) // Starting the loop at l=2 avoids an if statement. double current; double previous = work[1]; tensor[ti(1, 0, 0)] = 0.0; for (int l = 2; l < o1; l++) { // Initial condition for the inner loop is formation of T100(l-1). // Starting the inner loop at a=1 avoids an if statement. // T100(l-1) = 0.0 * T000(l) current = 0.0; int iw = il + (l - 1); work[iw] = current; // sb.append(format("double %s = 0.0;\n", term(1, 0, 0, l - 1))); for (int a = 1; a < l - 1; a++) { // T200(l-2) = 0.0 * T100(l-1) + (2 - 1) * T000(l-1) // T300(l-3) = 0.0 * T200(l-2) + (3 - 1) * T100(l-2) // ... // T(l-1)001 = 0.0 * T(l-2)002 + (l - 2) * T(l-3)002 // iw = (a - 1) * il + (l - a) current = a * work[iw - il]; iw += il - 1; // iw = (a + 1) * il + (l - a - 1) work[iw] = current; if (current != 0) { if (a > 2) { sb.append(format("double %s = %d * %s;\n", term(a + 1, 0, 0, l - a - 1), a, term(a - 1, 0, 0, l - a))); } else { sb.append(format("double %s = %s;\n", term(a + 1, 0, 0, l - a - 1), term(a - 1, 0, 0, l - a))); } } } // Store the Tl00 tensor (d/dx)^l // Tl00 = 0.0 * T(l-1)001 + (l - 1) * T(l-2)001 int index = ti(l, 0, 0); tensor[index] = (l - 1) * previous; previous = current; if (tensor[index] != 0) { if (l > 2) { sb.append(format("%s = %d * %s;\n", rlmn(l, 0, 0), (l - 1), term(l - 2, 0, 0, 1))); } else { sb.append(format("%s = %s;\n", rlmn(l, 0, 0), term(l - 2, 0, 0, 1))); } } } // Find (d/dx)^l * (d/dy)^m for l+m = 1..order (m > 0, n = 0) // Any (d/dy) term can be formed as: // Tlm0j = y * Tl(m-1)00(j+1) + (m-1) * Tl(m-2)00(j+1) for (int l = 0; l < order; l++) { // Store the m=1 tensor (d/dx)^l *(d/dy) // Tl10 = y * Tl001 previous = work[l * il + 1]; tensor[ti(l, 1, 0)] = 0.0; for (int m = 2; m + l < o1; m++) { // Tl10(m-1) = y * Tl00m; int iw = l * il + m; current = 0.0; iw += im - 1; work[iw] = current; for (int a = 1; a < m - 1; a++) { // Tl20(m-2) = 0.0 * Tl10(m-1) + (2 - 1) * Tl00(m-1) // Tl30(m-3) = 0.0 * Tl20(m-2) + (3 - 1) * Tl10(m-2) // ... // Tl(m-1)01 = 0.0 * Tl(m-2)02 + (m - 2) * Tl(m-3)02 current = a * work[iw - im]; iw += im - 1; work[iw] = current; if (current != 0) { if (a > 1) { sb.append(format("double %s = %d * %s;\n", term(l, a + 1, 0, m - a - 1), a, term(l, a - 1, 0, m - a))); } else { sb.append(format("double %s = %s;\n", term(l, a + 1, 0, m - a - 1), term(l, a - 1, 0, m - a))); } } } // Store the tensor (d/dx)^l * (d/dy)^m // Tlm0 = y * Tl(m-1)01 + (m - 1) * Tl(m-2)01 int index = ti(l, m, 0); tensor[index] = (m - 1) * previous; previous = current; if (tensor[index] != 0) { if (m > 2) { sb.append(format("%s = %d * %s;\n", rlmn(l, m, 0), (m - 1), term(l, m - 2, 0, 1))); } else { sb.append(format("%s = %s;\n", rlmn(l, m, 0), term(l, m - 2, 0, 1))); } } } } // Find (d/dx)^l * (d/dy)^m * (d/dz)^n for l+m+n = 1..order (n > 0) // Any (d/dz) term can be formed as: // Tlmnj = z * Tlm(n-1)(j+1) + (n-1) * Tlm(n-2)(j+1) for (int l = 0; l < order; l++) { for (int m = 0; m + l < order; m++) { // Store the n=1 tensor (d/dx)^l *(d/dy)^m * (d/dz) // Tlmn = z * Tlm01 final int lm = m + l; final int lilmim = l * il + m * im; previous = work[lilmim + 1]; int index = ti(l, m, 1); tensor[index] = z * previous; if (tensor[index] != 0) { sb.append(format("%s = z * %s;\n", rlmn(l, m, 1), term(l, m, 0, 1))); } for (int n = 2; lm + n < o1; n++) { // Tlm1(n-1) = z * Tlm0n; int iw = lilmim + n; current = z * work[iw]; iw += in - 1; work[iw] = current; if (current != 0) { sb.append(format("double %s = z * %s;\n", term(l, m, 1, n - 1), term(l, m, 0, n))); } final int n1 = n - 1; for (int a = 1; a < n1; a++) { // Tlm2(n-2) = z * Tlm1(n-1) + (2 - 1) * T1m0(n-1) // Tlm3(n-3) = z * Tlm2(n-2) + (3 - 1) * Tlm1(n-2) // ... // Tlm(n-1)1 = z * Tlm(n-2)2 + (n - 2) * Tlm(n-3)2 current = z * current + a * work[iw - in]; iw += in - 1; work[iw] = current; if (current != 0) { if (a > 1) { sb.append(format("double %s = z * %s + %d * %s;\n", term(l, m, a + 1, n - a - 1), term(l, m, a, n - a), a, term(l, m, a - 1, n - a))); } else { sb.append(format("double %s = z * %s + %s;\n", term(l, m, a + 1, n - a - 1), term(l, m, a, n - a), term(l, m, a - 1, n - a))); } } } // Store the tensor (d/dx)^l * (d/dy)^m * (d/dz)^n // Tlmn = z * Tlm(n-1)1 + (n - 1) * Tlm(n-2)1 index = ti(l, m, n); tensor[index] = z * current + n1 * previous; previous = current; if (tensor[index] != 0) { if (n > 2) { sb.append(format("%s = z * %s + %d * %s;\n", rlmn(l, m, n), term(l, m, n - 1, 1), (n - 1), term(l, m, n - 2, 1))); } else { sb.append(format("%s = z * %s + %s;\n", rlmn(l, m, n), term(l, m, n - 1, 1), term(l, m, n - 2, 1))); } } } } } return sb.toString(); } /** * Contract multipole moments with their respective electrostatic potential * derivatives. * * @param T array of electrostatic potential and partial derivatives * @param l apply (d/dx)^l to the potential * @param m apply (d/dy)^l to the potential * @param n apply (d/dz)^l to the potential * @return the contracted interaction. */ public double contract(double T[], int l, int m, int n) { double total = 0.0; double total2 = 0.0; total += qi * T[ti(l, m, n)]; total += dxi * T[ti(l + 1, m, n)]; total += dyi * T[ti(l, m + 1, n)]; total += dzi * T[ti(l, m, n + 1)]; total += qxxi * T[ti(l + 2, m, n)]; total += qyyi * T[ti(l, m + 2, n)]; total += qzzi * T[ti(l, m, n + 2)]; total2 += qxyi * T[ti(l + 1, m + 1, n)]; total2 += qxzi * T[ti(l + 1, m, n + 1)]; total2 += qyzi * T[ti(l, m + 1, n + 1)]; return total + 2.0 * total2; } /** * Contract multipole moments with their respective electrostatic potential * derivatives. * * @param T array of electrostatic potential and partial derivatives * @param l apply (d/dx)^l to the potential * @param m apply (d/dy)^l to the potential * @param n apply (d/dz)^l to the potential * @param sb the code will be appended to the StringBuilfer. * @return the contracted interaction. */ public double codeContract(double T[], int l, int m, int n, StringBuilder sb) { double total = 0.0; String name = term(l, m, n); sb.append(format("double %s = 0.0;\n", name)); StringBuilder sb1 = new StringBuilder(); double term = qi * T[ti(l, m, n)]; if (term != 0) { total += term; sb1.append(format("%s += q000i * T[%s];\n", name, tlmn(l, m, n))); } term = dxi * T[ti(l + 1, m, n)]; if (term != 0) { total += term; sb1.append(format("%s += q100i * T[%s];\n", name, tlmn(l + 1, m, n))); } term = dyi * T[ti(l, m + 1, n)]; if (term != 0) { total += term; sb1.append(format("%s += q010i * T[%s];\n", name, tlmn(l, m + 1, n))); } term = dzi * T[ti(l, m, n + 1)]; if (term != 0) { total += term; sb1.append(format("%s += q001i * T[%s];\n", name, tlmn(l, m, n + 1))); } StringBuilder traceSB = new StringBuilder(); double trace = 0.0; term = qxxi * T[ti(l + 2, m, n)]; if (term != 0) { trace += term; // logger.info(format(" Qxx: %16.15f T: %16.15f Term: %16.15f", q200, T[ti(l + 2, m, n)], term)); traceSB.append(format("%s += q200i * T[%s];\n", name, tlmn(l + 2, m, n))); } term = qyyi * T[ti(l, m + 2, n)]; if (term != 0) { trace += term; // logger.info(format(" Qyy: %16.15f T: %16.15f Term: %16.15f", q020, T[ti(l, m + 2, n)], term)); traceSB.append(format("%s += q020i * T[%s];\n", name, tlmn(l, m + 2, n))); } term = qzzi * T[ti(l, m, n + 2)]; if (term != 0) { trace += term; // logger.info(format(" Qzz: %16.15f T: %16.15f Term: %16.15f", q002, T[ti(l, m, n + 2)], term)); traceSB.append(format("%s += q002i * T[%s];\n", name, tlmn(l, m, n + 2))); } total += trace; if (total != 0) { sb.append(sb1.toString()); if (trace != 0) { //logger.info(format(" Trace: %16.15f", trace)); sb.append(traceSB); } } StringBuilder sb2 = new StringBuilder(); double total2 = 0.0; term = qxyi * T[ti(l + 1, m + 1, n)]; if (term != 0) { total2 += term; sb2.append(format("%s2 += q110i * T[%s];\n", name, tlmn(l + 1, m + 1, n))); } term = qxzi * T[ti(l + 1, m, n + 1)]; if (term != 0) { total2 += term; sb2.append(format("%s2 += q101i * T[%s];\n", name, tlmn(l + 1, m, n + 1))); } term = qyzi * T[ti(l, m + 1, n + 1)]; if (term != 0) { total2 += term; sb2.append(format("%s2 += q011i * T[%s];\n", name, tlmn(l, m + 1, n + 1))); } if (total2 != 0.0) { sb.append(format("double %s2 = 0.0;\n", name)); sb.append(sb2); total += 2.0 * total2; sb.append(format("%s += 2.0 * %s2;\n", name, name)); } return total; } /** * Collect the field at R due to Q multipole moments at the origin. * * @param T Electrostatic potential and partial derivatives * @param l apply (d/dx)^l to the potential * @param m apply (d/dy)^l to the potential * @param n apply (d/dz)^l to the potential */ public void field(double T[], int l, int m, int n) { E000 = contract(T, l, m, n); E100 = contract(T, l + 1, m, n); E010 = contract(T, l, m + 1, n); E001 = contract(T, l, m, n + 1); E200 = contract(T, l + 2, m, n); E020 = contract(T, l, m + 2, n); E002 = contract(T, l, m, n + 2); E110 = contract(T, l + 1, n + 1, m); E101 = contract(T, l + 1, m, n + 1); E011 = contract(T, l, m + 1, n + 1); } /** * Collect the field at R due to Q multipole moments at the origin. * * @param T Electrostatic potential and partial derivatives * @param l apply (d/dx)^l to the potential * @param m apply (d/dy)^l to the potential * @param n apply (d/dz)^l to the potential * @param sb */ public void codeField(double T[], int l, int m, int n, StringBuilder sb) { E000 = codeContract(T, l, m, n, sb); if (E000 != 0) { sb.append(format("e000 = %s;\n", term(l, m, n))); } E100 = codeContract(T, l + 1, m, n, sb); if (E100 != 0) { sb.append(format("e100 = %s;\n", term(l + 1, m, n))); } E010 = codeContract(T, l, m + 1, n, sb); if (E100 != 0) { sb.append(format("e010 = %s;\n", term(l, m + 1, n))); } E001 = codeContract(T, l, m, n + 1, sb); if (E001 != 0) { sb.append(format("e001 = %s;\n", term(l, m, n + 1))); } E200 = codeContract(T, l + 2, m, n, sb); if (E200 != 0) { sb.append(format("e200 = %s;\n", term(l + 2, m, n))); } E020 = codeContract(T, l, m + 2, n, sb); if (E020 != 0) { sb.append(format("e020 = %s;\n", term(l, m + 2, n))); } E002 = codeContract(T, l, m, n + 2, sb); if (E002 != 0) { sb.append(format("e002 = %s;\n", term(l, m, n + 2))); } E110 = codeContract(T, l + 1, m + 1, n, sb); if (E110 != 0) { sb.append(format("e110 = %s;\n", term(l + 1, m + 1, n))); } E101 = codeContract(T, l + 1, m, n + 1, sb); if (E101 != 0) { sb.append(format("e101 = %s;\n", term(l + 1, m, n + 1))); } E011 = codeContract(T, l, m + 1, n + 1, sb); if (E011 != 0) { sb.append(format("e011 = %s;\n", term(l, m + 1, n + 1))); } } public double codeInteract5(double r[], double Qi[], double Qk[], double Fi[], double Fk[], double Ti[], double Tk[]) { double T[] = new double[tensorCount(5)]; recursion(r, T); setMultipoleI(Qi); setMultipoleK(Qk); StringBuilder sb = new StringBuilder("\n\npublic void E5(double T[]) {\n"); codeField(T, 0, 0, 0, sb); sb.append("}\n"); sb.append("\n\npublic void Ex5(double T[]) {\n"); codeField(T, 1, 0, 0, sb); sb.append("}\n"); sb.append("\n\npublic void Ey5(double T[]) {\n"); codeField(T, 0, 1, 0, sb); sb.append("}\n"); sb.append("\n\npublic void Ez5(double T[]) {\n"); codeField(T, 0, 0, 1, sb); sb.append("}\n"); logger.log(Level.INFO, sb.toString()); return 0.0; } public double codeInteractQI5(double r[], double Qi[], double Qk[], double Fi[], double Fk[], double Ti[], double Tk[]) { double T[] = new double[tensorCount(5)]; assert (r[0] == 0.0 && r[1] == 0.0); recursionQI(r, T); setMultipoleI(Qi); setMultipoleK(Qk); StringBuilder sb = new StringBuilder("\n\npublic void EQI5(double T[]) {\n"); codeField(T, 0, 0, 0, sb); sb.append("}\n"); sb.append("\n\npublic void ExQI5(double T[]) {\n"); codeField(T, 1, 0, 0, sb); sb.append("}\n"); sb.append("\n\npublic void EyQI5(double T[]) {\n"); codeField(T, 0, 1, 0, sb); sb.append("}\n"); sb.append("\n\npublic void EzQI5(double T[]) {\n"); codeField(T, 0, 0, 1, sb); sb.append("}\n"); logger.log(Level.INFO, sb.toString()); return 0.0; } /** * Hard coded computation of all Cartesian multipole tensors up to 4th * order, in the global frame, which is sufficient for quadrupole-induced * dipole forces. */ public void order4() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; R000 = term0000; R100 = x * term0001; double term1001 = x * term0002; R200 = x * term1001 + term0001; double term1002 = x * term0003; double term2001 = x * term1002 + term0002; R300 = x * term2001 + 2 * term1001; double term1003 = x * term0004; double term2002 = x * term1003 + term0003; double term3001 = x * term2002 + 2 * term1002; R400 = x * term3001 + 3 * term2001; R010 = y * term0001; double term0101 = y * term0002; R020 = y * term0101 + term0001; double term0102 = y * term0003; double term0201 = y * term0102 + term0002; R030 = y * term0201 + 2 * term0101; double term0103 = y * term0004; double term0202 = y * term0103 + term0003; double term0301 = y * term0202 + 2 * term0102; R040 = y * term0301 + 3 * term0201; R110 = y * term1001; double term1101 = y * term1002; R120 = y * term1101 + term1001; double term1102 = y * term1003; double term1201 = y * term1102 + term1002; R130 = y * term1201 + 2 * term1101; R210 = y * term2001; double term2101 = y * term2002; R220 = y * term2101 + term2001; R310 = y * term3001; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; R011 = z * term0101; double term0111 = z * term0102; R012 = z * term0111 + term0101; double term0112 = z * term0103; double term0121 = z * term0112 + term0102; R013 = z * term0121 + 2 * term0111; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; R031 = z * term0301; R101 = z * term1001; double term1011 = z * term1002; R102 = z * term1011 + term1001; double term1012 = z * term1003; double term1021 = z * term1012 + term1002; R103 = z * term1021 + 2 * term1011; R111 = z * term1101; double term1111 = z * term1102; R112 = z * term1111 + term1101; R121 = z * term1201; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; R211 = z * term2101; R301 = z * term3001; } /** * Hard coded computation of all Cartesian multipole tensors up to 4th * order, based on a quasi-internal frame, which is sufficient for * quadrupole-induced dipole forces. */ public void order4QI() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; R000 = term0000; R200 = term0001; double term2001 = term0002; double term2002 = term0003; R400 = 3 * term2001; R020 = term0001; double term0201 = term0002; double term0202 = term0003; R040 = 3 * term0201; R220 = term2001; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; } /** * Hard coded computation of all Cartesian multipole tensors up to 5th * order, in the global frame, which is sufficient for quadrupole-quadrupole * forces. */ public void order5() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; double term0005 = work[5]; R000 = term0000; R100 = x * term0001; double term1001 = x * term0002; R200 = x * term1001 + term0001; double term1002 = x * term0003; double term2001 = x * term1002 + term0002; R300 = x * term2001 + 2 * term1001; double term1003 = x * term0004; double term2002 = x * term1003 + term0003; double term3001 = x * term2002 + 2 * term1002; R400 = x * term3001 + 3 * term2001; double term1004 = x * term0005; double term2003 = x * term1004 + term0004; double term3002 = x * term2003 + 2 * term1003; double term4001 = x * term3002 + 3 * term2002; R500 = x * term4001 + 4 * term3001; R010 = y * term0001; double term0101 = y * term0002; R020 = y * term0101 + term0001; double term0102 = y * term0003; double term0201 = y * term0102 + term0002; R030 = y * term0201 + 2 * term0101; double term0103 = y * term0004; double term0202 = y * term0103 + term0003; double term0301 = y * term0202 + 2 * term0102; R040 = y * term0301 + 3 * term0201; double term0104 = y * term0005; double term0203 = y * term0104 + term0004; double term0302 = y * term0203 + 2 * term0103; double term0401 = y * term0302 + 3 * term0202; R050 = y * term0401 + 4 * term0301; R110 = y * term1001; double term1101 = y * term1002; R120 = y * term1101 + term1001; double term1102 = y * term1003; double term1201 = y * term1102 + term1002; R130 = y * term1201 + 2 * term1101; double term1103 = y * term1004; double term1202 = y * term1103 + term1003; double term1301 = y * term1202 + 2 * term1102; R140 = y * term1301 + 3 * term1201; R210 = y * term2001; double term2101 = y * term2002; R220 = y * term2101 + term2001; double term2102 = y * term2003; double term2201 = y * term2102 + term2002; R230 = y * term2201 + 2 * term2101; R310 = y * term3001; double term3101 = y * term3002; R320 = y * term3101 + term3001; R410 = y * term4001; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; double term0014 = z * term0005; double term0023 = z * term0014 + term0004; double term0032 = z * term0023 + 2 * term0013; double term0041 = z * term0032 + 3 * term0022; R005 = z * term0041 + 4 * term0031; R011 = z * term0101; double term0111 = z * term0102; R012 = z * term0111 + term0101; double term0112 = z * term0103; double term0121 = z * term0112 + term0102; R013 = z * term0121 + 2 * term0111; double term0113 = z * term0104; double term0122 = z * term0113 + term0103; double term0131 = z * term0122 + 2 * term0112; R014 = z * term0131 + 3 * term0121; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; double term0212 = z * term0203; double term0221 = z * term0212 + term0202; R023 = z * term0221 + 2 * term0211; R031 = z * term0301; double term0311 = z * term0302; R032 = z * term0311 + term0301; R041 = z * term0401; R101 = z * term1001; double term1011 = z * term1002; R102 = z * term1011 + term1001; double term1012 = z * term1003; double term1021 = z * term1012 + term1002; R103 = z * term1021 + 2 * term1011; double term1013 = z * term1004; double term1022 = z * term1013 + term1003; double term1031 = z * term1022 + 2 * term1012; R104 = z * term1031 + 3 * term1021; R111 = z * term1101; double term1111 = z * term1102; R112 = z * term1111 + term1101; double term1112 = z * term1103; double term1121 = z * term1112 + term1102; R113 = z * term1121 + 2 * term1111; R121 = z * term1201; double term1211 = z * term1202; R122 = z * term1211 + term1201; R131 = z * term1301; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; double term2012 = z * term2003; double term2021 = z * term2012 + term2002; R203 = z * term2021 + 2 * term2011; R211 = z * term2101; double term2111 = z * term2102; R212 = z * term2111 + term2101; R221 = z * term2201; R301 = z * term3001; double term3011 = z * term3002; R302 = z * term3011 + term3001; R311 = z * term3101; R401 = z * term4001; } /** * Hard coded computation of all Cartesian multipole tensors up to 5th * order, based on a quasi-internal frame, which is sufficient for * quadrupole-quadrupole forces. */ public void order5QI() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; double term0005 = work[5]; R000 = term0000; R200 = term0001; double term2001 = term0002; double term2002 = term0003; R400 = 3 * term2001; double term2003 = term0004; double term4001 = 3 * term2002; R020 = term0001; double term0201 = term0002; double term0202 = term0003; R040 = 3 * term0201; double term0203 = term0004; double term0401 = 3 * term0202; R220 = term2001; double term2201 = term2002; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; double term0014 = z * term0005; double term0023 = z * term0014 + term0004; double term0032 = z * term0023 + 2 * term0013; double term0041 = z * term0032 + 3 * term0022; R005 = z * term0041 + 4 * term0031; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; double term0212 = z * term0203; double term0221 = z * term0212 + term0202; R023 = z * term0221 + 2 * term0211; R041 = z * term0401; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; double term2012 = z * term2003; double term2021 = z * term2012 + term2002; R203 = z * term2021 + 2 * term2011; R221 = z * term2201; R401 = z * term4001; } /** * Hard coded computation of all Cartesian multipole tensors up to 5th * order, in the global frame, which is sufficient for quadrupole-quadrupole * forces and orthogonal space sampling. */ public void order6() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; double term0005 = work[5]; double term0006 = work[6]; R000 = term0000; R100 = x * term0001; double term1001 = x * term0002; R200 = x * term1001 + term0001; double term1002 = x * term0003; double term2001 = x * term1002 + term0002; R300 = x * term2001 + 2 * term1001; double term1003 = x * term0004; double term2002 = x * term1003 + term0003; double term3001 = x * term2002 + 2 * term1002; R400 = x * term3001 + 3 * term2001; double term1004 = x * term0005; double term2003 = x * term1004 + term0004; double term3002 = x * term2003 + 2 * term1003; double term4001 = x * term3002 + 3 * term2002; R500 = x * term4001 + 4 * term3001; double term1005 = x * term0006; double term2004 = x * term1005 + term0005; double term3003 = x * term2004 + 2 * term1004; double term4002 = x * term3003 + 3 * term2003; double term5001 = x * term4002 + 4 * term3002; R600 = x * term5001 + 5 * term4001; R010 = y * term0001; double term0101 = y * term0002; R020 = y * term0101 + term0001; double term0102 = y * term0003; double term0201 = y * term0102 + term0002; R030 = y * term0201 + 2 * term0101; double term0103 = y * term0004; double term0202 = y * term0103 + term0003; double term0301 = y * term0202 + 2 * term0102; R040 = y * term0301 + 3 * term0201; double term0104 = y * term0005; double term0203 = y * term0104 + term0004; double term0302 = y * term0203 + 2 * term0103; double term0401 = y * term0302 + 3 * term0202; R050 = y * term0401 + 4 * term0301; double term0105 = y * term0006; double term0204 = y * term0105 + term0005; double term0303 = y * term0204 + 2 * term0104; double term0402 = y * term0303 + 3 * term0203; double term0501 = y * term0402 + 4 * term0302; R060 = y * term0501 + 5 * term0401; R110 = y * term1001; double term1101 = y * term1002; R120 = y * term1101 + term1001; double term1102 = y * term1003; double term1201 = y * term1102 + term1002; R130 = y * term1201 + 2 * term1101; double term1103 = y * term1004; double term1202 = y * term1103 + term1003; double term1301 = y * term1202 + 2 * term1102; R140 = y * term1301 + 3 * term1201; double term1104 = y * term1005; double term1203 = y * term1104 + term1004; double term1302 = y * term1203 + 2 * term1103; double term1401 = y * term1302 + 3 * term1202; R150 = y * term1401 + 4 * term1301; R210 = y * term2001; double term2101 = y * term2002; R220 = y * term2101 + term2001; double term2102 = y * term2003; double term2201 = y * term2102 + term2002; R230 = y * term2201 + 2 * term2101; double term2103 = y * term2004; double term2202 = y * term2103 + term2003; double term2301 = y * term2202 + 2 * term2102; R240 = y * term2301 + 3 * term2201; R310 = y * term3001; double term3101 = y * term3002; R320 = y * term3101 + term3001; double term3102 = y * term3003; double term3201 = y * term3102 + term3002; R330 = y * term3201 + 2 * term3101; R410 = y * term4001; double term4101 = y * term4002; R420 = y * term4101 + term4001; R510 = y * term5001; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; double term0014 = z * term0005; double term0023 = z * term0014 + term0004; double term0032 = z * term0023 + 2 * term0013; double term0041 = z * term0032 + 3 * term0022; R005 = z * term0041 + 4 * term0031; double term0015 = z * term0006; double term0024 = z * term0015 + term0005; double term0033 = z * term0024 + 2 * term0014; double term0042 = z * term0033 + 3 * term0023; double term0051 = z * term0042 + 4 * term0032; R006 = z * term0051 + 5 * term0041; R011 = z * term0101; double term0111 = z * term0102; R012 = z * term0111 + term0101; double term0112 = z * term0103; double term0121 = z * term0112 + term0102; R013 = z * term0121 + 2 * term0111; double term0113 = z * term0104; double term0122 = z * term0113 + term0103; double term0131 = z * term0122 + 2 * term0112; R014 = z * term0131 + 3 * term0121; double term0114 = z * term0105; double term0123 = z * term0114 + term0104; double term0132 = z * term0123 + 2 * term0113; double term0141 = z * term0132 + 3 * term0122; R015 = z * term0141 + 4 * term0131; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; double term0212 = z * term0203; double term0221 = z * term0212 + term0202; R023 = z * term0221 + 2 * term0211; double term0213 = z * term0204; double term0222 = z * term0213 + term0203; double term0231 = z * term0222 + 2 * term0212; R024 = z * term0231 + 3 * term0221; R031 = z * term0301; double term0311 = z * term0302; R032 = z * term0311 + term0301; double term0312 = z * term0303; double term0321 = z * term0312 + term0302; R033 = z * term0321 + 2 * term0311; R041 = z * term0401; double term0411 = z * term0402; R042 = z * term0411 + term0401; R051 = z * term0501; R101 = z * term1001; double term1011 = z * term1002; R102 = z * term1011 + term1001; double term1012 = z * term1003; double term1021 = z * term1012 + term1002; R103 = z * term1021 + 2 * term1011; double term1013 = z * term1004; double term1022 = z * term1013 + term1003; double term1031 = z * term1022 + 2 * term1012; R104 = z * term1031 + 3 * term1021; double term1014 = z * term1005; double term1023 = z * term1014 + term1004; double term1032 = z * term1023 + 2 * term1013; double term1041 = z * term1032 + 3 * term1022; R105 = z * term1041 + 4 * term1031; R111 = z * term1101; double term1111 = z * term1102; R112 = z * term1111 + term1101; double term1112 = z * term1103; double term1121 = z * term1112 + term1102; R113 = z * term1121 + 2 * term1111; double term1113 = z * term1104; double term1122 = z * term1113 + term1103; double term1131 = z * term1122 + 2 * term1112; R114 = z * term1131 + 3 * term1121; R121 = z * term1201; double term1211 = z * term1202; R122 = z * term1211 + term1201; double term1212 = z * term1203; double term1221 = z * term1212 + term1202; R123 = z * term1221 + 2 * term1211; R131 = z * term1301; double term1311 = z * term1302; R132 = z * term1311 + term1301; R141 = z * term1401; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; double term2012 = z * term2003; double term2021 = z * term2012 + term2002; R203 = z * term2021 + 2 * term2011; double term2013 = z * term2004; double term2022 = z * term2013 + term2003; double term2031 = z * term2022 + 2 * term2012; R204 = z * term2031 + 3 * term2021; R211 = z * term2101; double term2111 = z * term2102; R212 = z * term2111 + term2101; double term2112 = z * term2103; double term2121 = z * term2112 + term2102; R213 = z * term2121 + 2 * term2111; R221 = z * term2201; double term2211 = z * term2202; R222 = z * term2211 + term2201; R231 = z * term2301; R301 = z * term3001; double term3011 = z * term3002; R302 = z * term3011 + term3001; double term3012 = z * term3003; double term3021 = z * term3012 + term3002; R303 = z * term3021 + 2 * term3011; R311 = z * term3101; double term3111 = z * term3102; R312 = z * term3111 + term3101; R321 = z * term3201; R401 = z * term4001; double term4011 = z * term4002; R402 = z * term4011 + term4001; R411 = z * term4101; R501 = z * term5001; } /** * Hard coded computation of all Cartesian multipole tensors up to 6th * order, based on a quasi-internal frame, which is sufficient for * quadrupole-quadrupole forces and orthogonal space sampling. */ public void order6QI() { source(work); double term0000 = work[0]; double term0001 = work[1]; double term0002 = work[2]; double term0003 = work[3]; double term0004 = work[4]; double term0005 = work[5]; double term0006 = work[6]; R000 = term0000; R200 = term0001; double term2001 = term0002; double term2002 = term0003; R400 = 3 * term2001; double term2003 = term0004; double term4001 = 3 * term2002; double term2004 = term0005; double term4002 = 3 * term2003; R600 = 5 * term4001; R020 = term0001; double term0201 = term0002; double term0202 = term0003; R040 = 3 * term0201; double term0203 = term0004; double term0401 = 3 * term0202; double term0204 = term0005; double term0402 = 3 * term0203; R060 = 5 * term0401; R220 = term2001; double term2201 = term2002; double term2202 = term2003; R240 = 3 * term2201; R420 = term4001; R001 = z * term0001; double term0011 = z * term0002; R002 = z * term0011 + term0001; double term0012 = z * term0003; double term0021 = z * term0012 + term0002; R003 = z * term0021 + 2 * term0011; double term0013 = z * term0004; double term0022 = z * term0013 + term0003; double term0031 = z * term0022 + 2 * term0012; R004 = z * term0031 + 3 * term0021; double term0014 = z * term0005; double term0023 = z * term0014 + term0004; double term0032 = z * term0023 + 2 * term0013; double term0041 = z * term0032 + 3 * term0022; R005 = z * term0041 + 4 * term0031; double term0015 = z * term0006; double term0024 = z * term0015 + term0005; double term0033 = z * term0024 + 2 * term0014; double term0042 = z * term0033 + 3 * term0023; double term0051 = z * term0042 + 4 * term0032; R006 = z * term0051 + 5 * term0041; R021 = z * term0201; double term0211 = z * term0202; R022 = z * term0211 + term0201; double term0212 = z * term0203; double term0221 = z * term0212 + term0202; R023 = z * term0221 + 2 * term0211; double term0213 = z * term0204; double term0222 = z * term0213 + term0203; double term0231 = z * term0222 + 2 * term0212; R024 = z * term0231 + 3 * term0221; R041 = z * term0401; double term0411 = z * term0402; R042 = z * term0411 + term0401; R201 = z * term2001; double term2011 = z * term2002; R202 = z * term2011 + term2001; double term2012 = z * term2003; double term2021 = z * term2012 + term2002; R203 = z * term2021 + 2 * term2011; double term2013 = z * term2004; double term2022 = z * term2013 + term2003; double term2031 = z * term2022 + 2 * term2012; R204 = z * term2031 + 3 * term2021; R221 = z * term2201; double term2211 = z * term2202; R222 = z * term2211 + term2201; R401 = z * term4001; double term4011 = z * term4002; R402 = z * term4011 + term4001; } /** * Re-use R, Qi, and Qk from a previous call. * * @param Fi Output force on i. * @param Ti Output torque on i. * @param Tk Output torque on k. * * @return the energy. */ private double multipoleEnergyGlobal(double Fi[], double Ti[], double Tk[]) { multipoleIField(); double energy = dotMultipoleK(); // Torques multipoleKTorque(Tk); multipoleKField(); multipoleITorque(Ti); // Forces multipoleIdX(); Fi[0] = -dotMultipoleK(); multipoleIdY(); Fi[1] = -dotMultipoleK(); multipoleIdZ(); Fi[2] = -dotMultipoleK(); dEdF = -Fi[2]; /** * if (order > 5) { // multipoleIdZ2(); // d2EdF2 = dotMultipoleK(); } */ return energy; } /** * Re-use Qi and Qk from a previous call. * * @param Fi Output force on i. * @param Ti Output torque on i. * @param Tk Output torque on k. * * @return the energy. */ public double multipoleEnergyQI(double Fi[], double Ti[], double Tk[]) { // Compute the potential due to site I at site K. multipoleIFieldQI(); // Dot the potential, field, field gradient with multipole K. double energy = dotMultipoleK(); // Compute the torque on site K due to the field from site I. multipoleKTorque(Tk); // Compute the field at site I due to site K. multipoleKFieldQI(); // Compute the torque on site I due to the field from site K. multipoleITorque(Ti); // Compute the force on site I F = {-dE/dx, -dE/dy, -dE/dz}. multipoleIdXQI(); Fi[0] = -dotMultipoleK(); multipoleIdYQI(); Fi[1] = -dotMultipoleK(); multipoleIdZQI(); Fi[2] = -dotMultipoleK(); if (bufferCoordinates == COORDINATES.QI) { dEdF = -Fi[2]; if (order > 5) { multipoleIdZ2QI(); d2EdF2 = dotMultipoleK(); } } // Rotate the force and torques from the QI frame into the Global frame. qiToGlobal(Fi, Ti, Tk); if (bufferCoordinates == COORDINATES.GLOBAL) { // dEdL = dEdF = -Fi[2] dEdF = -Fi[2]; if (order > 5) { multipoleIdZ2QI(); d2EdF2 = dotMultipoleK(); } } return energy; } public void setBufferCoordinates(COORDINATES coord) { this.bufferCoordinates = coord; } public double[] getR() { return new double[] { x, y, z, R, R * R }; } public double getdEdF() { return dEdF; } public double getd2EdF2() { return d2EdF2; } private double polarizationEnergyGlobal(double scaleField, double scaleEnergy, double scaleMutual, double Fi[], double Ti[], double Tk[]) { // Find the potential, field, etc at k due to the induced dipole i. inducedIField(); // Energy of multipole k in the field of induced dipole i. double energy = scaleEnergy * dotMultipoleK(); /** * Get the induced-induced portion of the force. */ Fi[0] = -0.5 * scaleMutual * (pxk * E200 + pyk * E110 + pzk * E101); Fi[1] = -0.5 * scaleMutual * (pxk * E110 + pyk * E020 + pzk * E011); Fi[2] = -0.5 * scaleMutual * (pxk * E101 + pyk * E011 + pzk * E002); // Find the potential, field, etc at i due to the induced dipole k. inducedKField(); // Energy of multipole i in the field of induced dipole k. energy += scaleEnergy * dotMultipoleI(); /** * Get the induced-induced portion of the force. */ Fi[0] += 0.5 * scaleMutual * (pxi * E200 + pyi * E110 + pzi * E101); Fi[1] += 0.5 * scaleMutual * (pxi * E110 + pyi * E020 + pzi * E011); Fi[2] += 0.5 * scaleMutual * (pxi * E101 + pyi * E011 + pzi * E002); /** * Apply scale factors directly to induced dipole components for * efficiency and convenience in computing remaining force terms and * torques. */ scaleInduced(scaleField, scaleEnergy); // Find the potential, field, etc at k due to (ind + indCR) at i. inducedIFieldCR(); // Torque on multipole k. multipoleKTorque(Tk); // Find the potential, field, etc at i due to (ind + indCR) at k. inducedKFieldCR(); // Torque on multipole i. multipoleITorque(Ti); // Forces inducedIdX(); Fi[0] -= dotMultipoleK(); inducedIdY(); Fi[1] -= dotMultipoleK(); inducedIdZ(); Fi[2] -= dotMultipoleK(); inducedKdX(); Fi[0] -= dotMultipoleI(); inducedKdY(); Fi[1] -= dotMultipoleI(); inducedKdZ(); Fi[2] -= dotMultipoleI(); return energy; } public double polarizationEnergyQI(double scaleField, double scaleEnergy, double scaleMutual, double Fi[], double Ti[], double Tk[]) { // Find the potential, field, etc at k due to the induced dipole i. inducedIFieldQI(); // Energy of multipole k in the field of induced dipole i. double energy = scaleEnergy * dotMultipoleK(); /** * Get the induced-induced portion of the force. */ Fi[0] = -0.5 * scaleMutual * (pxk * E200 + pyk * E110 + pzk * E101); Fi[1] = -0.5 * scaleMutual * (pxk * E110 + pyk * E020 + pzk * E011); Fi[2] = -0.5 * scaleMutual * (pxk * E101 + pyk * E011 + pzk * E002); // Find the potential, field, etc at i due to the induced dipole k. inducedKFieldQI(); // Energy of multipole i in the field of induced dipole k. energy += scaleEnergy * dotMultipoleI(); /** * Get the induced-induced portion of the force. */ Fi[0] += 0.5 * scaleMutual * (pxi * E200 + pyi * E110 + pzi * E101); Fi[1] += 0.5 * scaleMutual * (pxi * E110 + pyi * E020 + pzi * E011); Fi[2] += 0.5 * scaleMutual * (pxi * E101 + pyi * E011 + pzi * E002); /** * Apply scale factors directly to induced dipole components for * efficiency and convenience in computing remaining force terms and * torques. */ scaleInduced(scaleField, scaleEnergy); // Find the potential, field, etc at k due to (ind + indCR) at i. inducedIFieldCRQI(); // Torque on multipole k. multipoleKTorque(Tk); // Find the potential, field, etc at i due to (ind + indCR) at k. inducedKFieldCRQI(); // Torque on multipole i. multipoleITorque(Ti); // Forces inducedIdXQI(); Fi[0] -= dotMultipoleK(); inducedIdYQI(); Fi[1] -= dotMultipoleK(); inducedIdZQI(); Fi[2] -= dotMultipoleK(); inducedKdXQI(); Fi[0] -= dotMultipoleI(); inducedKdYQI(); Fi[1] -= dotMultipoleI(); inducedKdZQI(); Fi[2] -= dotMultipoleI(); // Rotate the force and torques from the QI frame into the Global frame. qiToGlobal(Fi, Ti, Tk); return energy; } private void multipoleIField() { double term000 = qi * R000; term000 -= dxi * R100; term000 -= dyi * R010; term000 -= dzi * R001; term000 += qxxi * R200; term000 += qyyi * R020; term000 += qzzi * R002; term000 += qxyi * R110; term000 += qxzi * R101; term000 += qyzi * R011; E000 = term000; double term100 = qi * R100; term100 -= dxi * R200; term100 -= dyi * R110; term100 -= dzi * R101; term100 += qxxi * R300; term100 += qyyi * R120; term100 += qzzi * R102; term100 += qxyi * R210; term100 += qxzi * R201; term100 += qyzi * R111; E100 = term100; double term010 = qi * R010; term010 -= dxi * R110; term010 -= dyi * R020; term010 -= dzi * R011; term010 += qxxi * R210; term010 += qyyi * R030; term010 += qzzi * R012; term010 += qxyi * R120; term010 += qxzi * R111; term010 += qyzi * R021; E010 = term010; double term001 = qi * R001; term001 -= dxi * R101; term001 -= dyi * R011; term001 -= dzi * R002; term001 += qxxi * R201; term001 += qyyi * R021; term001 += qzzi * R003; term001 += qxyi * R111; term001 += qxzi * R102; term001 += qyzi * R012; E001 = term001; double term200 = qi * R200; term200 -= dxi * R300; term200 -= dyi * R210; term200 -= dzi * R201; term200 += qxxi * R400; term200 += qyyi * R220; term200 += qzzi * R202; term200 += qxyi * R310; term200 += qxzi * R301; term200 += qyzi * R211; E200 = term200; double term020 = qi * R020; term020 -= dxi * R120; term020 -= dyi * R030; term020 -= dzi * R021; term020 += qxxi * R220; term020 += qyyi * R040; term020 += qzzi * R022; term020 += qxyi * R130; term020 += qxzi * R121; term020 += qyzi * R031; E020 = term020; double term002 = qi * R002; term002 -= dxi * R102; term002 -= dyi * R012; term002 -= dzi * R003; term002 += qxxi * R202; term002 += qyyi * R022; term002 += qzzi * R004; term002 += qxyi * R112; term002 += qxzi * R103; term002 += qyzi * R013; E002 = term002; double term110 = qi * R110; term110 -= dxi * R210; term110 -= dyi * R120; term110 -= dzi * R111; term110 += qxxi * R310; term110 += qyyi * R130; term110 += qzzi * R112; term110 += qxyi * R220; term110 += qxzi * R211; term110 += qyzi * R121; E110 = term110; double term101 = qi * R101; term101 -= dxi * R201; term101 -= dyi * R111; term101 -= dzi * R102; term101 += qxxi * R301; term101 += qyyi * R121; term101 += qzzi * R103; term101 += qxyi * R211; term101 += qxzi * R202; term101 += qyzi * R112; E101 = term101; double term011 = qi * R011; term011 -= dxi * R111; term011 -= dyi * R021; term011 -= dzi * R012; term011 += qxxi * R211; term011 += qyyi * R031; term011 += qzzi * R013; term011 += qxyi * R121; term011 += qxzi * R112; term011 += qyzi * R022; E011 = term011; } private void multipoleKField() { double term000 = 0.0; term000 += qk * R000; term000 += dxk * R100; term000 += dyk * R010; term000 += dzk * R001; term000 += qxxk * R200; term000 += qyyk * R020; term000 += qzzk * R002; term000 += qxyk * R110; term000 += qxzk * R101; term000 += qyzk * R011; E000 = term000; double term100 = 0.0; term100 += qk * R100; term100 += dxk * R200; term100 += dyk * R110; term100 += dzk * R101; term100 += qxxk * R300; term100 += qyyk * R120; term100 += qzzk * R102; term100 += qxyk * R210; term100 += qxzk * R201; term100 += qyzk * R111; E100 = term100; double term010 = 0.0; term010 += qk * R010; term010 += dxk * R110; term010 += dyk * R020; term010 += dzk * R011; term010 += qxxk * R210; term010 += qyyk * R030; term010 += qzzk * R012; term010 += qxyk * R120; term010 += qxzk * R111; term010 += qyzk * R021; E010 = term010; double term001 = 0.0; term001 += qk * R001; term001 += dxk * R101; term001 += dyk * R011; term001 += dzk * R002; term001 += qxxk * R201; term001 += qyyk * R021; term001 += qzzk * R003; term001 += qxyk * R111; term001 += qxzk * R102; term001 += qyzk * R012; E001 = term001; double term200 = 0.0; term200 += qk * R200; term200 += dxk * R300; term200 += dyk * R210; term200 += dzk * R201; term200 += qxxk * R400; term200 += qyyk * R220; term200 += qzzk * R202; term200 += qxyk * R310; term200 += qxzk * R301; term200 += qyzk * R211; E200 = term200; double term020 = 0.0; term020 += qk * R020; term020 += dxk * R120; term020 += dyk * R030; term020 += dzk * R021; term020 += qxxk * R220; term020 += qyyk * R040; term020 += qzzk * R022; term020 += qxyk * R130; term020 += qxzk * R121; term020 += qyzk * R031; E020 = term020; double term002 = 0.0; term002 += qk * R002; term002 += dxk * R102; term002 += dyk * R012; term002 += dzk * R003; term002 += qxxk * R202; term002 += qyyk * R022; term002 += qzzk * R004; term002 += qxyk * R112; term002 += qxzk * R103; term002 += qyzk * R013; E002 = term002; double term110 = 0.0; term110 += qk * R110; term110 += dxk * R210; term110 += dyk * R120; term110 += dzk * R111; term110 += qxxk * R310; term110 += qyyk * R130; term110 += qzzk * R112; term110 += qxyk * R220; term110 += qxzk * R211; term110 += qyzk * R121; E110 = term110; double term101 = 0.0; term101 += qk * R101; term101 += dxk * R201; term101 += dyk * R111; term101 += dzk * R102; term101 += qxxk * R301; term101 += qyyk * R121; term101 += qzzk * R103; term101 += qxyk * R211; term101 += qxzk * R202; term101 += qyzk * R112; E101 = term101; double term011 = 0.0; term011 += qk * R011; term011 += dxk * R111; term011 += dyk * R021; term011 += dzk * R012; term011 += qxxk * R211; term011 += qyyk * R031; term011 += qzzk * R013; term011 += qxyk * R121; term011 += qxzk * R112; term011 += qyzk * R022; E011 = term011; } private void multipoleIdX() { double term100 = 0.0; term100 += qi * R100; term100 -= dxi * R200; term100 -= dyi * R110; term100 -= dzi * R101; term100 += qxxi * R300; term100 += qyyi * R120; term100 += qzzi * R102; term100 += qxyi * R210; term100 += qxzi * R201; term100 += qyzi * R111; E000 = term100; double term200 = 0.0; term200 += qi * R200; term200 -= dxi * R300; term200 -= dyi * R210; term200 -= dzi * R201; term200 += qxxi * R400; term200 += qyyi * R220; term200 += qzzi * R202; term200 += qxyi * R310; term200 += qxzi * R301; term200 += qyzi * R211; E100 = term200; double term110 = 0.0; term110 += qi * R110; term110 -= dxi * R210; term110 -= dyi * R120; term110 -= dzi * R111; term110 += qxxi * R310; term110 += qyyi * R130; term110 += qzzi * R112; term110 += qxyi * R220; term110 += qxzi * R211; term110 += qyzi * R121; E010 = term110; double term101 = 0.0; term101 += qi * R101; term101 -= dxi * R201; term101 -= dyi * R111; term101 -= dzi * R102; term101 += qxxi * R301; term101 += qyyi * R121; term101 += qzzi * R103; term101 += qxyi * R211; term101 += qxzi * R202; term101 += qyzi * R112; E001 = term101; double term300 = 0.0; term300 += qi * R300; term300 -= dxi * R400; term300 -= dyi * R310; term300 -= dzi * R301; term300 += qxxi * R500; term300 += qyyi * R320; term300 += qzzi * R302; term300 += qxyi * R410; term300 += qxzi * R401; term300 += qyzi * R311; E200 = term300; double term120 = 0.0; term120 += qi * R120; term120 -= dxi * R220; term120 -= dyi * R130; term120 -= dzi * R121; term120 += qxxi * R320; term120 += qyyi * R140; term120 += qzzi * R122; term120 += qxyi * R230; term120 += qxzi * R221; term120 += qyzi * R131; E020 = term120; double term102 = 0.0; term102 += qi * R102; term102 -= dxi * R202; term102 -= dyi * R112; term102 -= dzi * R103; term102 += qxxi * R302; term102 += qyyi * R122; term102 += qzzi * R104; term102 += qxyi * R212; term102 += qxzi * R203; term102 += qyzi * R113; E002 = term102; double term210 = 0.0; term210 += qi * R210; term210 -= dxi * R310; term210 -= dyi * R220; term210 -= dzi * R211; term210 += qxxi * R410; term210 += qyyi * R230; term210 += qzzi * R212; term210 += qxyi * R320; term210 += qxzi * R311; term210 += qyzi * R221; E110 = term210; double term201 = 0.0; term201 += qi * R201; term201 -= dxi * R301; term201 -= dyi * R211; term201 -= dzi * R202; term201 += qxxi * R401; term201 += qyyi * R221; term201 += qzzi * R203; term201 += qxyi * R311; term201 += qxzi * R302; term201 += qyzi * R212; E101 = term201; double term111 = 0.0; term111 += qi * R111; term111 -= dxi * R211; term111 -= dyi * R121; term111 -= dzi * R112; term111 += qxxi * R311; term111 += qyyi * R131; term111 += qzzi * R113; term111 += qxyi * R221; term111 += qxzi * R212; term111 += qyzi * R122; E011 = term111; } private void multipoleIdY() { double term010 = 0.0; term010 += qi * R010; term010 -= dxi * R110; term010 -= dyi * R020; term010 -= dzi * R011; term010 += qxxi * R210; term010 += qyyi * R030; term010 += qzzi * R012; term010 += qxyi * R120; term010 += qxzi * R111; term010 += qyzi * R021; E000 = term010; double term110 = 0.0; term110 += qi * R110; term110 -= dxi * R210; term110 -= dyi * R120; term110 -= dzi * R111; term110 += qxxi * R310; term110 += qyyi * R130; term110 += qzzi * R112; term110 += qxyi * R220; term110 += qxzi * R211; term110 += qyzi * R121; E100 = term110; double term020 = 0.0; term020 += qi * R020; term020 -= dxi * R120; term020 -= dyi * R030; term020 -= dzi * R021; term020 += qxxi * R220; term020 += qyyi * R040; term020 += qzzi * R022; term020 += qxyi * R130; term020 += qxzi * R121; term020 += qyzi * R031; E010 = term020; double term011 = 0.0; term011 += qi * R011; term011 -= dxi * R111; term011 -= dyi * R021; term011 -= dzi * R012; term011 += qxxi * R211; term011 += qyyi * R031; term011 += qzzi * R013; term011 += qxyi * R121; term011 += qxzi * R112; term011 += qyzi * R022; E001 = term011; double term210 = 0.0; term210 += qi * R210; term210 -= dxi * R310; term210 -= dyi * R220; term210 -= dzi * R211; term210 += qxxi * R410; term210 += qyyi * R230; term210 += qzzi * R212; term210 += qxyi * R320; term210 += qxzi * R311; term210 += qyzi * R221; E200 = term210; double term030 = 0.0; term030 += qi * R030; term030 -= dxi * R130; term030 -= dyi * R040; term030 -= dzi * R031; term030 += qxxi * R230; term030 += qyyi * R050; term030 += qzzi * R032; term030 += qxyi * R140; term030 += qxzi * R131; term030 += qyzi * R041; E020 = term030; double term012 = 0.0; term012 += qi * R012; term012 -= dxi * R112; term012 -= dyi * R022; term012 -= dzi * R013; term012 += qxxi * R212; term012 += qyyi * R032; term012 += qzzi * R014; term012 += qxyi * R122; term012 += qxzi * R113; term012 += qyzi * R023; E002 = term012; double term120 = 0.0; term120 += qi * R120; term120 -= dxi * R220; term120 -= dyi * R130; term120 -= dzi * R121; term120 += qxxi * R320; term120 += qyyi * R140; term120 += qzzi * R122; term120 += qxyi * R230; term120 += qxzi * R221; term120 += qyzi * R131; E110 = term120; double term111 = 0.0; term111 += qi * R111; term111 -= dxi * R211; term111 -= dyi * R121; term111 -= dzi * R112; term111 += qxxi * R311; term111 += qyyi * R131; term111 += qzzi * R113; term111 += qxyi * R221; term111 += qxzi * R212; term111 += qyzi * R122; E101 = term111; double term021 = 0.0; term021 += qi * R021; term021 -= dxi * R121; term021 -= dyi * R031; term021 -= dzi * R022; term021 += qxxi * R221; term021 += qyyi * R041; term021 += qzzi * R023; term021 += qxyi * R131; term021 += qxzi * R122; term021 += qyzi * R032; E011 = term021; } private void multipoleIdZ() { double term001 = 0.0; term001 += qi * R001; term001 -= dxi * R101; term001 -= dyi * R011; term001 -= dzi * R002; term001 += qxxi * R201; term001 += qyyi * R021; term001 += qzzi * R003; term001 += qxyi * R111; term001 += qxzi * R102; term001 += qyzi * R012; E000 = term001; double term101 = 0.0; term101 += qi * R101; term101 -= dxi * R201; term101 -= dyi * R111; term101 -= dzi * R102; term101 += qxxi * R301; term101 += qyyi * R121; term101 += qzzi * R103; term101 += qxyi * R211; term101 += qxzi * R202; term101 += qyzi * R112; E100 = term101; double term011 = 0.0; term011 += qi * R011; term011 -= dxi * R111; term011 -= dyi * R021; term011 -= dzi * R012; term011 += qxxi * R211; term011 += qyyi * R031; term011 += qzzi * R013; term011 += qxyi * R121; term011 += qxzi * R112; term011 += qyzi * R022; E010 = term011; double term002 = 0.0; term002 += qi * R002; term002 -= dxi * R102; term002 -= dyi * R012; term002 -= dzi * R003; term002 += qxxi * R202; term002 += qyyi * R022; term002 += qzzi * R004; term002 += qxyi * R112; term002 += qxzi * R103; term002 += qyzi * R013; E001 = term002; double term201 = 0.0; term201 += qi * R201; term201 -= dxi * R301; term201 -= dyi * R211; term201 -= dzi * R202; term201 += qxxi * R401; term201 += qyyi * R221; term201 += qzzi * R203; term201 += qxyi * R311; term201 += qxzi * R302; term201 += qyzi * R212; E200 = term201; double term021 = 0.0; term021 += qi * R021; term021 -= dxi * R121; term021 -= dyi * R031; term021 -= dzi * R022; term021 += qxxi * R221; term021 += qyyi * R041; term021 += qzzi * R023; term021 += qxyi * R131; term021 += qxzi * R122; term021 += qyzi * R032; E020 = term021; double term003 = 0.0; term003 += qi * R003; term003 -= dxi * R103; term003 -= dyi * R013; term003 -= dzi * R004; term003 += qxxi * R203; term003 += qyyi * R023; term003 += qzzi * R005; term003 += qxyi * R113; term003 += qxzi * R104; term003 += qyzi * R014; E002 = term003; double term111 = 0.0; term111 += qi * R111; term111 -= dxi * R211; term111 -= dyi * R121; term111 -= dzi * R112; term111 += qxxi * R311; term111 += qyyi * R131; term111 += qzzi * R113; term111 += qxyi * R221; term111 += qxzi * R212; term111 += qyzi * R122; E110 = term111; double term102 = 0.0; term102 += qi * R102; term102 -= dxi * R202; term102 -= dyi * R112; term102 -= dzi * R103; term102 += qxxi * R302; term102 += qyyi * R122; term102 += qzzi * R104; term102 += qxyi * R212; term102 += qxzi * R203; term102 += qyzi * R113; E101 = term102; double term012 = 0.0; term012 += qi * R012; term012 -= dxi * R112; term012 -= dyi * R022; term012 -= dzi * R013; term012 += qxxi * R212; term012 += qyyi * R032; term012 += qzzi * R014; term012 += qxyi * R122; term012 += qxzi * R113; term012 += qyzi * R023; E011 = term012; } private void multipoleIFieldQI() { double term000 = qi * R000; term000 -= dzi * R001; term000 += qxxi * R200; term000 += qyyi * R020; term000 += qzzi * R002; E000 = term000; double term100 = -dxi * R200; term100 += qxzi * R201; E100 = term100; double term010 = -dyi * R020; term010 += qyzi * R021; E010 = term010; double term001 = qi * R001; term001 -= dzi * R002; term001 += qxxi * R201; term001 += qyyi * R021; term001 += qzzi * R003; E001 = term001; double term200 = qi * R200; term200 -= dzi * R201; term200 += qxxi * R400; term200 += qyyi * R220; term200 += qzzi * R202; E200 = term200; double term020 = qi * R020; term020 -= dzi * R021; term020 += qxxi * R220; term020 += qyyi * R040; term020 += qzzi * R022; E020 = term020; double term002 = qi * R002; term002 -= dzi * R003; term002 += qxxi * R202; term002 += qyyi * R022; term002 += qzzi * R004; E002 = term002; double term110 = qxyi * R220; E110 = term110; double term101 = -dxi * R201; term101 += qxzi * R202; E101 = term101; double term011 = -dyi * R021; term011 += qyzi * R022; E011 = term011; } private void multipoleKFieldQI() { double term000 = qk * R000; term000 += dzk * R001; term000 += qxxk * R200; term000 += qyyk * R020; term000 += qzzk * R002; E000 = term000; double term100 = dxk * R200; term100 += qxzk * R201; E100 = term100; double term010 = dyk * R020; term010 += qyzk * R021; E010 = term010; double term001 = qk * R001; term001 += dzk * R002; term001 += qxxk * R201; term001 += qyyk * R021; term001 += qzzk * R003; E001 = term001; double term200 = qk * R200; term200 += dzk * R201; term200 += qxxk * R400; term200 += qyyk * R220; term200 += qzzk * R202; E200 = term200; double term020 = qk * R020; term020 += dzk * R021; term020 += qxxk * R220; term020 += qyyk * R040; term020 += qzzk * R022; E020 = term020; double term002 = qk * R002; term002 += dzk * R003; term002 += qxxk * R202; term002 += qyyk * R022; term002 += qzzk * R004; E002 = term002; double term110 = qxyk * R220; E110 = term110; double term101 = dxk * R201; term101 += qxzk * R202; E101 = term101; double term011 = dyk * R021; term011 += qyzk * R022; E011 = term011; } private void multipoleIdXQI() { double term100 = -dxi * R200; term100 += qxzi * R201; E000 = term100; double term200 = qi * R200; term200 -= dzi * R201; term200 += qxxi * R400; term200 += qyyi * R220; term200 += qzzi * R202; E100 = term200; double term110 = qxyi * R220; E010 = term110; double term101 = -dxi * R201; term101 += qxzi * R202; E001 = term101; double term300 = -dxi * R400; term300 += qxzi * R401; E200 = term300; double term120 = -dxi * R220; term120 += qxzi * R221; E020 = term120; double term102 = -dxi * R202; term102 += qxzi * R203; E002 = term102; double term210 = -dyi * R220; term210 += qyzi * R221; E110 = term210; double term201 = qi * R201; term201 -= dzi * R202; term201 += qxxi * R401; term201 += qyyi * R221; term201 += qzzi * R203; E101 = term201; double term111 = qxyi * R221; E011 = term111; } private void multipoleIdYQI() { double term010 = -dyi * R020; term010 += qyzi * R021; E000 = term010; double term110 = qxyi * R220; E100 = term110; double term020 = qi * R020; term020 -= dzi * R021; term020 += qxxi * R220; term020 += qyyi * R040; term020 += qzzi * R022; E010 = term020; double term011 = -dyi * R021; term011 += qyzi * R022; E001 = term011; double term210 = -dyi * R220; term210 += qyzi * R221; E200 = term210; double term030 = -dyi * R040; term030 += qyzi * R041; E020 = term030; double term012 = -dyi * R022; term012 += qyzi * R023; E002 = term012; double term120 = -dxi * R220; term120 += qxzi * R221; E110 = term120; double term111 = qxyi * R221; E101 = term111; double term021 = qi * R021; term021 -= dzi * R022; term021 += qxxi * R221; term021 += qyyi * R041; term021 += qzzi * R023; E011 = term021; } private void multipoleIdZQI() { double term001 = qi * R001; term001 -= dzi * R002; term001 += qxxi * R201; term001 += qyyi * R021; term001 += qzzi * R003; E000 = term001; double term101 = -dxi * R201; term101 += qxzi * R202; E100 = term101; double term011 = -dyi * R021; term011 += qyzi * R022; E010 = term011; double term002 = qi * R002; term002 -= dzi * R003; term002 += qxxi * R202; term002 += qyyi * R022; term002 += qzzi * R004; E001 = term002; double term201 = qi * R201; term201 -= dzi * R202; term201 += qxxi * R401; term201 += qyyi * R221; term201 += qzzi * R203; E200 = term201; double term021 = qi * R021; term021 -= dzi * R022; term021 += qxxi * R221; term021 += qyyi * R041; term021 += qzzi * R023; E020 = term021; double term003 = qi * R003; term003 -= dzi * R004; term003 += qxxi * R203; term003 += qyyi * R023; term003 += qzzi * R005; E002 = term003; double term111 = qxyi * R221; E110 = term111; double term102 = -dxi * R202; term102 += qxzi * R203; E101 = term102; double term012 = -dyi * R022; term012 += qyzi * R023; E011 = term012; } private void multipoleIdZ2QI() { if (order < 6) { logger.severe("Use higher order tensor for lambda derivatives."); } double term001 = qi * R002; term001 -= dzi * R003; term001 += qxxi * R202; term001 += qyyi * R022; term001 += qzzi * R004; E000 = term001; double term101 = -dxi * R202; term101 += qxzi * R203; E100 = term101; double term011 = -dyi * R022; term011 += qyzi * R023; E010 = term011; double term002 = qi * R003; term002 -= dzi * R004; term002 += qxxi * R203; term002 += qyyi * R023; term002 += qzzi * R005; E001 = term002; double term201 = qi * R202; term201 -= dzi * R203; term201 += qxxi * R402; term201 += qyyi * R222; term201 += qzzi * R204; E200 = term201; double term021 = qi * R022; term021 -= dzi * R023; term021 += qxxi * R222; term021 += qyyi * R042; term021 += qzzi * R024; E020 = term021; double term003 = qi * R004; term003 -= dzi * R005; term003 += qxxi * R204; term003 += qyyi * R024; term003 += qzzi * R006; E002 = term003; double term111 = qxyi * R222; E110 = term111; double term102 = -dxi * R203; term102 += qxzi * R204; E101 = term102; double term012 = -dyi * R023; term012 += qyzi * R024; E011 = term012; } private void inducedIField() { double term000 = -uxi * R100; term000 -= uyi * R010; term000 -= uzi * R001; E000 = term000; double term100 = -uxi * R200; term100 -= uyi * R110; term100 -= uzi * R101; E100 = term100; double term010 = -uxi * R110; term010 -= uyi * R020; term010 -= uzi * R011; E010 = term010; double term001 = -uxi * R101; term001 -= uyi * R011; term001 -= uzi * R002; E001 = term001; double term200 = -uxi * R300; term200 -= uyi * R210; term200 -= uzi * R201; E200 = term200; double term020 = -uxi * R120; term020 -= uyi * R030; term020 -= uzi * R021; E020 = term020; double term002 = -uxi * R102; term002 -= uyi * R012; term002 -= uzi * R003; E002 = term002; double term110 = -uxi * R210; term110 -= uyi * R120; term110 -= uzi * R111; E110 = term110; double term101 = -uxi * R201; term101 -= uyi * R111; term101 -= uzi * R102; E101 = term101; double term011 = -uxi * R111; term011 -= uyi * R021; term011 -= uzi * R012; E011 = term011; } private void inducedKField() { double term000 = uxk * R100; term000 += uyk * R010; term000 += uzk * R001; E000 = term000; double term100 = uxk * R200; term100 += uyk * R110; term100 += uzk * R101; E100 = term100; double term010 = uxk * R110; term010 += uyk * R020; term010 += uzk * R011; E010 = term010; double term001 = uxk * R101; term001 += uyk * R011; term001 += uzk * R002; E001 = term001; double term200 = uxk * R300; term200 += uyk * R210; term200 += uzk * R201; E200 = term200; double term020 = uxk * R120; term020 += uyk * R030; term020 += uzk * R021; E020 = term020; double term002 = uxk * R102; term002 += uyk * R012; term002 += uzk * R003; E002 = term002; double term110 = uxk * R210; term110 += uyk * R120; term110 += uzk * R111; E110 = term110; double term101 = uxk * R201; term101 += uyk * R111; term101 += uzk * R102; E101 = term101; double term011 = uxk * R111; term011 += uyk * R021; term011 += uzk * R012; E011 = term011; } private void inducedIFieldCR() { double term000 = -sxi * R100; term000 -= syi * R010; term000 -= szi * R001; E000 = term000; double term100 = -sxi * R200; term100 -= syi * R110; term100 -= szi * R101; E100 = term100; double term010 = -sxi * R110; term010 -= syi * R020; term010 -= szi * R011; E010 = term010; double term001 = -sxi * R101; term001 -= syi * R011; term001 -= szi * R002; E001 = term001; double term200 = -sxi * R300; term200 -= syi * R210; term200 -= szi * R201; E200 = term200; double term020 = -sxi * R120; term020 -= syi * R030; term020 -= szi * R021; E020 = term020; double term002 = -sxi * R102; term002 -= syi * R012; term002 -= szi * R003; E002 = term002; double term110 = -sxi * R210; term110 -= syi * R120; term110 -= szi * R111; E110 = term110; double term101 = -sxi * R201; term101 -= syi * R111; term101 -= szi * R102; E101 = term101; double term011 = -sxi * R111; term011 -= syi * R021; term011 -= szi * R012; E011 = term011; } private void inducedKFieldCR() { double term000 = sxk * R100; term000 += syk * R010; term000 += szk * R001; E000 = term000; double term100 = sxk * R200; term100 += syk * R110; term100 += szk * R101; E100 = term100; double term010 = sxk * R110; term010 += syk * R020; term010 += szk * R011; E010 = term010; double term001 = sxk * R101; term001 += syk * R011; term001 += szk * R002; E001 = term001; double term200 = sxk * R300; term200 += syk * R210; term200 += szk * R201; E200 = term200; double term020 = sxk * R120; term020 += syk * R030; term020 += szk * R021; E020 = term020; double term002 = sxk * R102; term002 += syk * R012; term002 += szk * R003; E002 = term002; double term110 = sxk * R210; term110 += syk * R120; term110 += szk * R111; E110 = term110; double term101 = sxk * R201; term101 += syk * R111; term101 += szk * R102; E101 = term101; double term011 = sxk * R111; term011 += syk * R021; term011 += szk * R012; E011 = term011; } private void inducedIdX() { double term100 = 0.0; term100 -= sxi * R200; term100 -= syi * R110; term100 -= szi * R101; E000 = term100; double term200 = 0.0; term200 -= sxi * R300; term200 -= syi * R210; term200 -= szi * R201; E100 = term200; double term110 = 0.0; term110 -= sxi * R210; term110 -= syi * R120; term110 -= szi * R111; E010 = term110; double term101 = 0.0; term101 -= sxi * R201; term101 -= syi * R111; term101 -= szi * R102; E001 = term101; double term300 = 0.0; term300 -= sxi * R400; term300 -= syi * R310; term300 -= szi * R301; E200 = term300; double term120 = 0.0; term120 -= sxi * R220; term120 -= syi * R130; term120 -= szi * R121; E020 = term120; double term102 = 0.0; term102 -= sxi * R202; term102 -= syi * R112; term102 -= szi * R103; E002 = term102; double term210 = 0.0; term210 -= sxi * R310; term210 -= syi * R220; term210 -= szi * R211; E110 = term210; double term201 = 0.0; term201 -= sxi * R301; term201 -= syi * R211; term201 -= szi * R202; E101 = term201; double term111 = 0.0; term111 -= sxi * R211; term111 -= syi * R121; term111 -= szi * R112; E011 = term111; } private void inducedIdY() { double term010 = 0.0; term010 -= sxi * R110; term010 -= syi * R020; term010 -= szi * R011; E000 = term010; double term110 = 0.0; term110 -= sxi * R210; term110 -= syi * R120; term110 -= szi * R111; E100 = term110; double term020 = 0.0; term020 -= sxi * R120; term020 -= syi * R030; term020 -= szi * R021; E010 = term020; double term011 = 0.0; term011 -= sxi * R111; term011 -= syi * R021; term011 -= szi * R012; E001 = term011; double term210 = 0.0; term210 -= sxi * R310; term210 -= syi * R220; term210 -= szi * R211; E200 = term210; double term030 = 0.0; term030 -= sxi * R130; term030 -= syi * R040; term030 -= szi * R031; E020 = term030; double term012 = 0.0; term012 -= sxi * R112; term012 -= syi * R022; term012 -= szi * R013; E002 = term012; double term120 = 0.0; term120 -= sxi * R220; term120 -= syi * R130; term120 -= szi * R121; E110 = term120; double term111 = 0.0; term111 -= sxi * R211; term111 -= syi * R121; term111 -= szi * R112; E101 = term111; double term021 = 0.0; term021 -= sxi * R121; term021 -= syi * R031; term021 -= szi * R022; E011 = term021; } private void inducedIdZ() { double term001 = 0.0; term001 -= sxi * R101; term001 -= syi * R011; term001 -= szi * R002; E000 = term001; double term101 = 0.0; term101 -= sxi * R201; term101 -= syi * R111; term101 -= szi * R102; E100 = term101; double term011 = 0.0; term011 -= sxi * R111; term011 -= syi * R021; term011 -= szi * R012; E010 = term011; double term002 = 0.0; term002 -= sxi * R102; term002 -= syi * R012; term002 -= szi * R003; E001 = term002; double term201 = 0.0; term201 -= sxi * R301; term201 -= syi * R211; term201 -= szi * R202; E200 = term201; double term021 = 0.0; term021 -= sxi * R121; term021 -= syi * R031; term021 -= szi * R022; E020 = term021; double term003 = 0.0; term003 -= sxi * R103; term003 -= syi * R013; term003 -= szi * R004; E002 = term003; double term111 = 0.0; term111 -= sxi * R211; term111 -= syi * R121; term111 -= szi * R112; E110 = term111; double term102 = 0.0; term102 -= sxi * R202; term102 -= syi * R112; term102 -= szi * R103; E101 = term102; double term012 = 0.0; term012 -= sxi * R112; term012 -= syi * R022; term012 -= szi * R013; E011 = term012; } private void inducedKdX() { double term100 = 0.0; term100 += sxk * R200; term100 += syk * R110; term100 += szk * R101; E000 = term100; double term200 = 0.0; term200 += sxk * R300; term200 += syk * R210; term200 += szk * R201; E100 = term200; double term110 = 0.0; term110 += sxk * R210; term110 += syk * R120; term110 += szk * R111; E010 = term110; double term101 = 0.0; term101 += sxk * R201; term101 += syk * R111; term101 += szk * R102; E001 = term101; double term300 = 0.0; term300 += sxk * R400; term300 += syk * R310; term300 += szk * R301; E200 = term300; double term120 = 0.0; term120 += sxk * R220; term120 += syk * R130; term120 += szk * R121; E020 = term120; double term102 = 0.0; term102 += sxk * R202; term102 += syk * R112; term102 += szk * R103; E002 = term102; double term210 = 0.0; term210 += sxk * R310; term210 += syk * R220; term210 += szk * R211; E110 = term210; double term201 = 0.0; term201 += sxk * R301; term201 += syk * R211; term201 += szk * R202; E101 = term201; double term111 = 0.0; term111 += sxk * R211; term111 += syk * R121; term111 += szk * R112; E011 = term111; } private void inducedKdY() { double term010 = 0.0; term010 += sxk * R110; term010 += syk * R020; term010 += szk * R011; E000 = term010; double term110 = 0.0; term110 += sxk * R210; term110 += syk * R120; term110 += szk * R111; E100 = term110; double term020 = 0.0; term020 += sxk * R120; term020 += syk * R030; term020 += szk * R021; E010 = term020; double term011 = 0.0; term011 += sxk * R111; term011 += syk * R021; term011 += szk * R012; E001 = term011; double term210 = 0.0; term210 += sxk * R310; term210 += syk * R220; term210 += szk * R211; E200 = term210; double term030 = 0.0; term030 += sxk * R130; term030 += syk * R040; term030 += szk * R031; E020 = term030; double term012 = 0.0; term012 += sxk * R112; term012 += syk * R022; term012 += szk * R013; E002 = term012; double term120 = 0.0; term120 += sxk * R220; term120 += syk * R130; term120 += szk * R121; E110 = term120; double term111 = 0.0; term111 += sxk * R211; term111 += syk * R121; term111 += szk * R112; E101 = term111; double term021 = 0.0; term021 += sxk * R121; term021 += syk * R031; term021 += szk * R022; E011 = term021; } private void inducedKdZ() { double term001 = 0.0; term001 += sxk * R101; term001 += syk * R011; term001 += szk * R002; E000 = term001; double term101 = 0.0; term101 += sxk * R201; term101 += syk * R111; term101 += szk * R102; E100 = term101; double term011 = 0.0; term011 += sxk * R111; term011 += syk * R021; term011 += szk * R012; E010 = term011; double term002 = 0.0; term002 += sxk * R102; term002 += syk * R012; term002 += szk * R003; E001 = term002; double term201 = 0.0; term201 += sxk * R301; term201 += syk * R211; term201 += szk * R202; E200 = term201; double term021 = 0.0; term021 += sxk * R121; term021 += syk * R031; term021 += szk * R022; E020 = term021; double term003 = 0.0; term003 += sxk * R103; term003 += syk * R013; term003 += szk * R004; E002 = term003; double term111 = 0.0; term111 += sxk * R211; term111 += syk * R121; term111 += szk * R112; E110 = term111; double term102 = 0.0; term102 += sxk * R202; term102 += syk * R112; term102 += szk * R103; E101 = term102; double term012 = 0.0; term012 += sxk * R112; term012 += syk * R022; term012 += szk * R013; E011 = term012; } private void inducedIFieldQI() { E000 = -uzi * R001; E100 = -uxi * R200; E010 = -uyi * R020; E001 = -uzi * R002; E200 = -uzi * R201; E020 = -uzi * R021; E002 = -uzi * R003; E110 = 0.0; E101 = -uxi * R201; E011 = -uyi * R021; } private void inducedKFieldQI() { E000 = uzk * R001; E100 = uxk * R200; E010 = uyk * R020; E001 = uzk * R002; E200 = uzk * R201; E020 = uzk * R021; E002 = uzk * R003; E110 = 0.0; E101 = uxk * R201; E011 = uyk * R021; } private void inducedIFieldCRQI() { E000 = -szi * R001; E100 = -sxi * R200; E010 = -syi * R020; E001 = -szi * R002; E200 = -szi * R201; E020 = -szi * R021; E002 = -szi * R003; E110 = 0.0; E101 = -sxi * R201; E011 = -syi * R021; } private void inducedKFieldCRQI() { E000 = szk * R001; E100 = sxk * R200; E010 = syk * R020; E001 = szk * R002; E200 = szk * R201; E020 = szk * R021; E002 = szk * R003; E110 = 0.0; E101 = sxk * R201; E011 = syk * R021; } private void inducedIdXQI() { E000 = -sxi * R200; E100 = -szi * R201; E010 = 0.0; E001 = -sxi * R201; E200 = -sxi * R400; E020 = -sxi * R220; E002 = -sxi * R202; E110 = -syi * R220; E101 = -szi * R202; E011 = 0.0; } private void inducedIdYQI() { E000 = -syi * R020; E100 = 0.0; E010 = -szi * R021; E001 = -syi * R021; E200 = -syi * R220; E020 = -syi * R040; E002 = -syi * R022; E110 = -sxi * R220; E101 = 0.0; E011 = -szi * R022; } private void inducedIdZQI() { E000 = -szi * R002; E100 = -sxi * R201; E010 = -syi * R021; E001 = -szi * R003; E200 = -szi * R202; E020 = -szi * R022; E002 = -szi * R004; E110 = 0.0; E101 = -sxi * R202; E011 = -syi * R022; } private void inducedKdXQI() { E000 = sxk * R200; E100 = szk * R201; E010 = 0.0; E001 = sxk * R201; E200 = sxk * R400; E020 = sxk * R220; E002 = sxk * R202; E110 = syk * R220; E101 = szk * R202; E011 = 0.0; } private void inducedKdYQI() { E000 = syk * R020; E100 = 0.0; E010 = szk * R021; E001 = syk * R021; E200 = syk * R220; E020 = syk * R040; E002 = syk * R022; E110 = sxk * R220; E101 = 0.0; E011 = szk * R022; } private void inducedKdZQI() { E000 = szk * R002; E100 = sxk * R201; E010 = syk * R021; E001 = szk * R003; E200 = szk * R202; E020 = szk * R022; E002 = szk * R004; E110 = 0.0; E101 = sxk * R202; E011 = syk * R022; } private void multipoleITorque(double torque[]) { // Torque on dipole moments due to the field. double dx = dyi * E001 - dzi * E010; double dy = dzi * E100 - dxi * E001; double dz = dxi * E010 - dyi * E100; // Torque on quadrupole moments due to the gradient of the field. double qx = qxyi * E101 + 2.0 * qyyi * E011 + qyzi * E002 - (qxzi * E110 + qyzi * E020 + 2.0 * qzzi * E011); double qy = qxzi * E200 + qyzi * E110 + 2.0 * qzzi * E101 - (2.0 * qxxi * E101 + qxyi * E011 + qxzi * E002); double qz = 2.0 * qxxi * E110 + qxyi * E020 + qxzi * E011 - (qxyi * E200 + 2.0 * qyyi * E110 + qyzi * E101); torque[0] = dx - qx; torque[1] = dy - qy; torque[2] = dz - qz; } private void multipoleKTorque(double torque[]) { // Torque on dipole moments due to the field. double dx = dyk * E001 - dzk * E010; double dy = dzk * E100 - dxk * E001; double dz = dxk * E010 - dyk * E100; // Torque on quadrupole moments due to the gradkent of the fkeld. double qx = qxyk * E101 + 2.0 * qyyk * E011 + qyzk * E002 - (qxzk * E110 + qyzk * E020 + 2.0 * qzzk * E011); double qy = qxzk * E200 + qyzk * E110 + 2.0 * qzzk * E101 - (2.0 * qxxk * E101 + qxyk * E011 + qxzk * E002); double qz = 2.0 * qxxk * E110 + qxyk * E020 + qxzk * E011 - (qxyk * E200 + 2.0 * qyyk * E110 + qyzk * E101); torque[0] = -(dx + qx); torque[1] = -(dy + qy); torque[2] = -(dz + qz); } private double dotMultipoleK() { double total = qk * E000; total += dxk * E100; total += dyk * E010; total += dzk * E001; total += qxxk * E200; total += qyyk * E020; total += qzzk * E002; total += qxyk * E110; total += qxzk * E101; total += qyzk * E011; return total; } private double dotMultipoleI() { double total = qi * E000; total -= dxi * E100; total -= dyi * E010; total -= dzi * E001; total += qxxi * E200; total += qyyi * E020; total += qzzi * E002; total += qxyi * E110; total += qxzi * E101; total += qyzi * E011; return total; } private void validate(boolean print) { double[] qiVals = new double[] { qi, dxi, dyi, dzi, qxxi, qyyi, qzzi, qxyi, qxzi, qyzi }; double[] qkVals = new double[] { qk, dxk, dyk, dzk, qxxk, qyyk, qzzk, qxyk, qxzk, qyzk }; double[] exxxVals = new double[] { E000, E100, E010, E001, E200, E020, E002, E110, E101, E011 }; double[] rxxxVals = new double[] { R000, R100, R010, R001, R200, R020, R002, R110, R101, R011, R300, R120, R102, R210, R201, R111, R210, R030, R012, R120, R111, R021, R201, R021, R003, R111, R102, R012, R400, R220, R202, R310, R301, R211, R220, R040, R022, R130, R121, R031, R202, R022, R004, R112, R103, R013 }; double[][] all = new double[][] { qiVals, qkVals, exxxVals, rxxxVals }; if (!print) { for (int i = 0; i < all.length; i++) { for (int j = 0; j < all[i].length; j++) { if (Double.isNaN(all[i][j])) { logger.warning(format("MT::validate(): NaN @ (%d,%d)", i, j)); } if (Double.isInfinite(all[i][j])) { logger.warning(format("MT::validate(): Inf @ (%d,%d)", i, j)); } } } } else { logger.info(format("MT::ALL_VALS: %s", formArr(all))); } } public static void main(String args[]) { if (args == null || args.length < 4) { logger.info(" Usage: java ffx.numerics.MultipoleTensor order dx dy dz"); } int order = Integer.parseInt(args[0]); double dx = Double.parseDouble(args[1]); double dy = Double.parseDouble(args[2]); double dz = Double.parseDouble(args[3]); double r[] = { dx, dy, dz }; double n2 = 710643; double cycles = 10; logger.info(format(" 6th Order Tensor Count: %d", tensorCount(6))); MultipoleTensor multipoleTensor = new MultipoleTensor(OPERATOR.SCREENED_COULOMB, COORDINATES.GLOBAL, order, 1e-6); double[] Fi = new double[3]; double[] Ti = new double[3]; double[] Tk = new double[3]; double[] Qi = new double[] { 0.11, 0.21, 0.31, 0.41, -0.51, -0.61, 1.12, 0.71, 0.81, 0.91 }; double[] Qk = new double[] { 0.11, 0.21, 0.31, 0.41, -0.51, -0.61, 1.12, 0.70, 0.81, 0.91 }; for (int j = 0; j < cycles; j++) { long timeGlobalT = -System.nanoTime(); for (int i = 0; i < n2; i++) { r[0] = Math.random(); r[1] = Math.random(); r[2] = Math.random(); multipoleTensor.setR(r); multipoleTensor.order6(); } timeGlobalT += System.nanoTime(); long timeGlobal = -System.nanoTime(); for (int i = 0; i < n2; i++) { r[0] = Math.random(); r[1] = Math.random(); r[2] = Math.random(); multipoleTensor.setR(r); multipoleTensor.setMultipoles(Qi, Qk); double e = multipoleTensor.multipoleEnergy(Fi, Ti, Tk); if (Double.isNaN(e) || Double.isInfinite(e)) { multipoleTensor.nanWarning(e, r, Qi, Qk, Fi, Ti, Tk); } } timeGlobal += System.nanoTime(); multipoleTensor.setCoordinateSystem(COORDINATES.QI); long timeQIT = -System.nanoTime(); for (int i = 0; i < n2; i++) { r[0] = Math.random(); r[1] = Math.random(); r[2] = Math.random(); multipoleTensor.setR_QI(r); multipoleTensor.order6QI(); } timeQIT += System.nanoTime(); long timeQI = -System.nanoTime(); for (int i = 0; i < n2; i++) { r[0] = 0.0; r[1] = 0.0; r[2] = Math.random(); multipoleTensor.setR_QI(r); multipoleTensor.setMultipoles(Qi, Qk); double e = multipoleTensor.multipoleEnergy(Fi, Ti, Tk); if (Double.isNaN(e) || Double.isInfinite(e)) { multipoleTensor.nanWarning(e, r, Qi, Qk, Fi, Ti, Tk); } } timeQI += System.nanoTime(); logger.info(String.format("\n Global Frame: %6.4f %6.4f\n QI: %6.4f %6.4f\n", timeGlobalT * 1.0e-9, timeGlobal * 1.0e-9, timeQIT * 1.0e-9, timeQI * 1.0e-9)); } /** * double tensors[] = new double[MultipoleTensor.tensorCount(order)]; * String string = multipoleTensor.codeTensorRecursion(r, tensors); * logger.info(" Java Code:\n" + string); string = * multipoleTensor.codeTensorRecursionQI(r, tensors); logger.info(" Java * Code:\n" + string); */ } /** * Convenience method for writing out intermediate terms in the recursion. * * @param l number of d/dx partial derivatives. * @param m number of d/dx partial derivatives. * @param n number of d/dx partial derivatives. * @param j the jth intermediate term. * @return a String of the form <code>termlmnj</code>. */ private static String term(int l, int m, int n, int j) { return String.format("term%d%d%d%d", l, m, n, j); } /** * Convenience method for writing out intermediate terms in the recursion. * * @param l number of d/dx partial derivatives. * @param m number of d/dx partial derivatives. * @param n number of d/dx partial derivatives. * @param j the jth intermediate term. * @return a String of the form <code>termlmnj</code>. */ private static String term(int l, int m, int n) { return String.format("term%d%d%d", l, m, n); } /** * Convenience method for writing out tensor indeces. * * @param l number of d/dx partial derivatives. * @param m number of d/dx partial derivatives. * @param n number of d/dx partial derivatives. * @return a String of the form <code>tlmn</code>. */ private static String tlmn(int l, int m, int n) { return String.format("t%d%d%d", l, m, n); } /** * Convenience method for writing out tensor indeces. * * @param l number of d/dx partial derivatives. * @param m number of d/dx partial derivatives. * @param n number of d/dx partial derivatives. * @return a String of the form <code>Rlmn</code>. */ private static String rlmn(int l, int m, int n) { return String.format("R%d%d%d", l, m, n); } public void getTensor(double T[]) { if (T == null || order < 0) { return; } // l + m + n = 0 (1) T[t000] = R000; if (order < 1) { return; } // l + m + n = 1 (3) 4 T[t100] = R100; T[t010] = R010; T[t001] = R001; if (order < 2) { return; } // l + m + n = 2 (6) 10 T[t200] = R200; T[t020] = R020; T[t002] = R002; T[t110] = R110; T[t101] = R101; T[t011] = R011; if (order < 3) { return; } // l + m + n = 3 (10) 20 T[t300] = R300; T[t030] = R030; T[t003] = R003; T[t210] = R210; T[t201] = R201; T[t120] = R120; T[t021] = R021; T[t102] = R102; T[t012] = R012; T[t111] = R111; if (order < 4) { return; } // l + m + n = 4 (15) 35 T[t400] = R400; T[t040] = R040; T[t004] = R004; T[t310] = R310; T[t301] = R301; T[t130] = R130; T[t031] = R031; T[t103] = R103; T[t013] = R013; T[t220] = R220; T[t202] = R202; T[t022] = R022; T[t211] = R211; T[t121] = R121; T[t112] = R112; if (order < 5) { return; } // l + m + n = 5 (21) 56 T[t500] = R500; T[t050] = R050; T[t005] = R005; T[t410] = R410; T[t401] = R401; T[t140] = R140; T[t041] = R041; T[t104] = R104; T[t014] = R014; T[t320] = R320; T[t302] = R302; T[t230] = R230; T[t032] = R032; T[t203] = R203; T[t023] = R023; T[t311] = R311; T[t131] = R131; T[t113] = R113; T[t221] = R221; T[t212] = R212; T[t122] = R122; } public void setTensor(double T[]) { if (T == null || order < 0) { return; } // l + m + n = 0 (1) R000 = T[t000]; if (order < 1) { return; } // l + m + n = 1 (3) 4 R100 = T[t100]; R010 = T[t010]; R001 = T[t001]; if (order < 2) { return; } // l + m + n = 2 (6) 10 R200 = T[t200]; R020 = T[t020]; R002 = T[t002]; R110 = T[t110]; R101 = T[t101]; R011 = T[t011]; if (order < 3) { return; } // l + m + n = 3 (10) 20 R300 = T[t300]; R030 = T[t030]; R003 = T[t003]; R210 = T[t210]; R201 = T[t201]; R120 = T[t120]; R021 = T[t021]; R102 = T[t102]; R012 = T[t012]; R111 = T[t111]; if (order < 4) { return; } // l + m + n = 4 (15) 35 R400 = T[t400]; R040 = T[t040]; R004 = T[t004]; R310 = T[t310]; R301 = T[t301]; R130 = T[t130]; R031 = T[t031]; R103 = T[t103]; R013 = T[t013]; R220 = T[t220]; R202 = T[t202]; R022 = T[t022]; R211 = T[t211]; R121 = T[t121]; R112 = T[t112]; if (order < 5) { return; } // l + m + n = 5 (21) 56 R500 = T[t500]; R050 = T[t050]; R005 = T[t005]; R410 = T[t410]; R401 = T[t401]; R140 = T[t140]; R041 = T[t041]; R104 = T[t104]; R014 = T[t014]; R320 = T[t320]; R302 = T[t302]; R230 = T[t230]; R032 = T[t032]; R203 = T[t203]; R023 = T[t023]; R311 = T[t311]; R131 = T[t131]; R113 = T[t113]; R221 = T[t221]; R212 = T[t212]; R122 = T[t122]; } private void setQIRotationMatrix(double dx, double dy, double dz) { double zAxis[] = { dx, dy, dz }; double xAxis[] = { dx + 1.0, dy, dz }; norm(zAxis, zAxis); ir02 = zAxis[0]; ir12 = zAxis[1]; ir22 = zAxis[2]; double dot = dot(xAxis, zAxis); scalar(zAxis, dot, zAxis); diff(xAxis, zAxis, xAxis); norm(xAxis, xAxis); ir00 = xAxis[0]; ir10 = xAxis[1]; ir20 = xAxis[2]; ir01 = ir20 * ir12 - ir10 * ir22; ir11 = ir00 * ir22 - ir20 * ir02; ir21 = ir10 * ir02 - ir00 * ir12; // Set the forward elements as the transpose of the inverse matrix. r00 = ir00; r11 = ir11; r22 = ir22; r01 = ir10; r02 = ir20; r10 = ir01; r12 = ir21; r20 = ir02; r21 = ir12; } private static final double ONE_THIRD = 1.0 / 3.0; private static final double TWO_THIRD = 2.0 / 3.0; public void setMultipoles(double Qi[], double Qk[]) { switch (coordinates) { case GLOBAL: default: setMultipoleI(Qi); setMultipoleK(Qk); break; case QI: multipoleItoQI(Qi); multipoleKtoQI(Qk); break; } } public void setMultipolesQI(double Qi[], double Qk[]) { multipoleItoQI(Qi); multipoleKtoQI(Qk); } public void setDipoles(double Ui[], double UiCR[], double Uk[], double UkCR[]) { switch (coordinates) { case GLOBAL: default: setDipoleI(Ui, UiCR); setDipoleK(Uk, UkCR); break; case QI: dipoleItoQI(Ui, UiCR); dipoleKtoQI(Uk, UkCR); break; } } public void setDipolesQI(double Ui[], double UiCR[], double Uk[], double UkCR[]) { dipoleItoQI(Ui, UiCR); dipoleKtoQI(Uk, UkCR); } private void setDipoleI(double Ui[], double UiCR[]) { uxi = Ui[0]; uyi = Ui[1]; uzi = Ui[2]; pxi = UiCR[0]; pyi = UiCR[1]; pzi = UiCR[2]; } private void setDipoleK(double Uk[], double UkCR[]) { uxk = Uk[0]; uyk = Uk[1]; uzk = Uk[2]; pxk = UkCR[0]; pyk = UkCR[1]; pzk = UkCR[2]; } private void scaleInduced(double scaleField, double scaleEnergy) { uxi *= scaleEnergy; uyi *= scaleEnergy; uzi *= scaleEnergy; pxi *= scaleField; pyi *= scaleField; pzi *= scaleField; sxi = 0.5 * (uxi + pxi); syi = 0.5 * (uyi + pyi); szi = 0.5 * (uzi + pzi); uxk *= scaleEnergy; uyk *= scaleEnergy; uzk *= scaleEnergy; pxk *= scaleField; pyk *= scaleField; pzk *= scaleField; sxk = 0.5 * (uxk + pxk); syk = 0.5 * (uyk + pyk); szk = 0.5 * (uzk + pzk); } private void setMultipoleI(double Qi[]) { qi = Qi[0]; dxi = Qi[1]; dyi = Qi[2]; dzi = Qi[3]; qxxi = Qi[4] * ONE_THIRD; qyyi = Qi[5] * ONE_THIRD; qzzi = Qi[6] * ONE_THIRD; qxyi = Qi[7] * TWO_THIRD; qxzi = Qi[8] * TWO_THIRD; qyzi = Qi[9] * TWO_THIRD; } private void setMultipoleK(double Qk[]) { qk = Qk[0]; dxk = Qk[1]; dyk = Qk[2]; dzk = Qk[3]; qxxk = Qk[4] * ONE_THIRD; qyyk = Qk[5] * ONE_THIRD; qzzk = Qk[6] * ONE_THIRD; qxyk = Qk[7] * TWO_THIRD; qxzk = Qk[8] * TWO_THIRD; qyzk = Qk[9] * TWO_THIRD; } private void dipoleItoQI(double Ui[], double UiCR[]) { double dx = Ui[0]; double dy = Ui[1]; double dz = Ui[2]; uxi = r00 * dx + r01 * dy + r02 * dz; uyi = r10 * dx + r11 * dy + r12 * dz; uzi = r20 * dx + r21 * dy + r22 * dz; dx = UiCR[0]; dy = UiCR[1]; dz = UiCR[2]; pxi = r00 * dx + r01 * dy + r02 * dz; pyi = r10 * dx + r11 * dy + r12 * dz; pzi = r20 * dx + r21 * dy + r22 * dz; } private void dipoleKtoQI(double Uk[], double UkCR[]) { double dx = Uk[0]; double dy = Uk[1]; double dz = Uk[2]; uxk = r00 * dx + r01 * dy + r02 * dz; uyk = r10 * dx + r11 * dy + r12 * dz; uzk = r20 * dx + r21 * dy + r22 * dz; dx = UkCR[0]; dy = UkCR[1]; dz = UkCR[2]; pxk = r00 * dx + r01 * dy + r02 * dz; pyk = r10 * dx + r11 * dy + r12 * dz; pzk = r20 * dx + r21 * dy + r22 * dz; } private void multipoleItoQI(double Qi[]) { qi = Qi[0]; double dx = Qi[1]; double dy = Qi[2]; double dz = Qi[3]; dxi = r00 * dx + r01 * dy + r02 * dz; dyi = r10 * dx + r11 * dy + r12 * dz; dzi = r20 * dx + r21 * dy + r22 * dz; double qxx = Qi[4] * ONE_THIRD; double qyy = Qi[5] * ONE_THIRD; double qzz = Qi[6] * ONE_THIRD; double qxy = Qi[7] * ONE_THIRD; double qxz = Qi[8] * ONE_THIRD; double qyz = Qi[9] * ONE_THIRD; // i=0, j=0 // qij r0k * r00 * qkx + r01 * qky + r02 * qkz qxxi = r00 * (r00 * qxx + r01 * qxy + r02 * qxz) + r01 * (r00 * qxy + r01 * qyy + r02 * qyz) + r02 * (r00 * qxz + r01 * qyz + r02 * qzz); // i=0, j=1 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qxyi = r00 * (r10 * qxx + r11 * qxy + r12 * qxz) + r01 * (r10 * qxy + r11 * qyy + r12 * qyz) + r02 * (r10 * qxz + r11 * qyz + r12 * qzz); qxyi *= 2.0; // i=0, j=2 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qxzi = r00 * (r20 * qxx + r21 * qxy + r22 * qxz) + r01 * (r20 * qxy + r21 * qyy + r22 * qyz) + r02 * (r20 * qxz + r21 * qyz + r22 * qzz); qxzi *= 2.0; // i=1, j=1 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qyyi = r10 * (r10 * qxx + r11 * qxy + r12 * qxz) + r11 * (r10 * qxy + r11 * qyy + r12 * qyz) + r12 * (r10 * qxz + r11 * qyz + r12 * qzz); // i=1, j=2 // qij r1k * r20 * qkx + r21 * qky + r22 * qkz qyzi = r10 * (r20 * qxx + r21 * qxy + r22 * qxz) + r11 * (r20 * qxy + r21 * qyy + r22 * qyz) + r12 * (r20 * qxz + r21 * qyz + r22 * qzz); qyzi *= 2.0; // i=2, j=2 // qij r2k * r20 * qkx + r21 * qky + r22 * qkz qzzi = r20 * (r20 * qxx + r21 * qxy + r22 * qxz) + r21 * (r20 * qxy + r21 * qyy + r22 * qyz) + r22 * (r20 * qxz + r21 * qyz + r22 * qzz); } private void multipoleKtoQI(double Qk[]) { qk = Qk[0]; double dx = Qk[1]; double dy = Qk[2]; double dz = Qk[3]; dxk = r00 * dx + r01 * dy + r02 * dz; dyk = r10 * dx + r11 * dy + r12 * dz; dzk = r20 * dx + r21 * dy + r22 * dz; double qxx = Qk[4] * ONE_THIRD; double qyy = Qk[5] * ONE_THIRD; double qzz = Qk[6] * ONE_THIRD; double qxy = Qk[7] * ONE_THIRD; double qxz = Qk[8] * ONE_THIRD; double qyz = Qk[9] * ONE_THIRD; // i=0, j=0 // qij r0k * r00 * qkx + r01 * qky + r02 * qkz qxxk = r00 * (r00 * qxx + r01 * qxy + r02 * qxz) + r01 * (r00 * qxy + r01 * qyy + r02 * qyz) + r02 * (r00 * qxz + r01 * qyz + r02 * qzz); // i=0, j=1 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qxyk = r00 * (r10 * qxx + r11 * qxy + r12 * qxz) + r01 * (r10 * qxy + r11 * qyy + r12 * qyz) + r02 * (r10 * qxz + r11 * qyz + r12 * qzz); qxyk *= 2.0; // i=0, j=2 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qxzk = r00 * (r20 * qxx + r21 * qxy + r22 * qxz) + r01 * (r20 * qxy + r21 * qyy + r22 * qyz) + r02 * (r20 * qxz + r21 * qyz + r22 * qzz); qxzk *= 2.0; // i=1, j=1 // qij rik * rj0 * qkx + rj1 * qky + rj2 * qkz qyyk = r10 * (r10 * qxx + r11 * qxy + r12 * qxz) + r11 * (r10 * qxy + r11 * qyy + r12 * qyz) + r12 * (r10 * qxz + r11 * qyz + r12 * qzz); // i=1, j=2 // qij r1k * r20 * qkx + r21 * qky + r22 * qkz qyzk = r10 * (r20 * qxx + r21 * qxy + r22 * qxz) + r11 * (r20 * qxy + r21 * qyy + r22 * qyz) + r12 * (r20 * qxz + r21 * qyz + r22 * qzz); qyzk *= 2.0; // i=2, j=2 // qij r2k * r20 * qkx + r21 * qky + r22 * qkz qzzk = r20 * (r20 * qxx + r21 * qxy + r22 * qxz) + r21 * (r20 * qxy + r21 * qyy + r22 * qyz) + r22 * (r20 * qxz + r21 * qyz + r22 * qzz); } private void qiToGlobal(double v1[], double v2[], double v3[]) { double vx = v1[0]; double vy = v1[1]; double vz = v1[2]; v1[0] = ir00 * vx + ir01 * vy + ir02 * vz; v1[1] = ir10 * vx + ir11 * vy + ir12 * vz; v1[2] = ir20 * vx + ir21 * vy + ir22 * vz; vx = v2[0]; vy = v2[1]; vz = v2[2]; v2[0] = ir00 * vx + ir01 * vy + ir02 * vz; v2[1] = ir10 * vx + ir11 * vy + ir12 * vz; v2[2] = ir20 * vx + ir21 * vy + ir22 * vz; vx = v3[0]; vy = v3[1]; vz = v3[2]; v3[0] = ir00 * vx + ir01 * vy + ir02 * vz; v3[1] = ir10 * vx + ir11 * vy + ir12 * vz; v3[2] = ir20 * vx + ir21 * vy + ir22 * vz; } private void nanWarning(double energy, double[] r, double[] Qi, double[] Qk, double[] Fi, double[] Ti, double[] Tk) { StringBuilder sb = new StringBuilder(); if (Double.isInfinite(energy)) { sb.append(format("DotK infinite: \n")); } else if (Double.isNaN(energy)) { sb.append(format("DotK was NaN: \n")); } sb.append(format(" r: %s\n Qi: %s\n Qk: %s\n Fi: %s\n Ti: %s\n Tk: %s\n", formArr(r), formArr(Qi), formArr(Qk), formArr(Fi), formArr(Ti), formArr(Tk))); double total = qk * E000; double total2 = qxyi * E110; // sb.append(format("DotK components:" // + "\n (1) %.4f %.4f %.4f %.4f %.4f\n (2) %.4f %.4f %.4f %.4f %.4f" // + "\n (3) %.4f %.4f %.4f %.4f %.4f\n (4) %.4f %.4f %.4f %.4f %.4f" // + "\n (5) %.4f %.4f %.4f", // E000, E100, E010, E001, E200, // E020, E002, E110, E101, E011, // qi, dxi, dyi, dzi, qxxi, // qyyi, qzzi, qxyi, qxzi, qyzi, // qk, dxk, dyk, dzk, qxxk, // qyyk, qzzk, qxyi, qxzk, qyzk, // total, total2, total + 2.0*total2)); double[] Exxx = new double[] { E000, E100, E010, E001, E200, E020, E002, E110, E010, E001 }; double[] mpoleI = new double[] { qi, dxi, dyi, dzi, qxxi, qyyi, qzzi, qxyi, qxzi, qyzi }; double[] mpoleK = new double[] { qk, dxk, dyk, dzk, qxxk, qyyk, qzzk, qxyk, qxzk, qyzk }; sb.append(format("DotK components:\n Exxx: %s\n mpoleI: %s\n mpoleK: %s", formArr(Exxx), formArr(mpoleI), formArr(mpoleK))); (new ArithmeticException()).printStackTrace(); logger.warning(sb.toString()); } // Helper method for logging distance and multipole arrays. private static String formArr(double[] x) { StringBuilder sb = new StringBuilder(); sb.append("["); for (int i = 0; i < x.length; i++) { sb.append(String.format("%.4f", x[i])); if (i + 1 < x.length) { sb.append(", "); } } sb.append("]"); return sb.toString(); } // Helper method for logging distance and multipole arrays. private static String formArr(double[][] x) { StringBuilder sb = new StringBuilder(); for (int i = 0; i < x.length; i++) { sb.append("\n["); if (i == 0) { sb.append("["); } for (int j = 0; j < x[i].length; j++) { sb.append(String.format("%.4f", x[i][j])); if (j + 1 < x[i].length) { sb.append(", "); } } sb.append("]"); if (i + 1 < x.length) { sb.append("; "); } } sb.append("]\n"); return sb.toString(); } // Rotation Matrix from Global to QI. private double r00, r01, r02; private double r10, r11, r12; private double r20, r21, r22; // Rotation Matrix from QI to Global. private double ir00, ir01, ir02; private double ir10, ir11, ir12; private double ir20, ir21, ir22; // Multipole components for atom i. private double qi; private double dxi; private double dyi; private double dzi; private double qxxi; private double qyyi; private double qzzi; private double qxyi; private double qxzi; private double qyzi; // Induced dipole components for atom i. private double uxi; private double uyi; private double uzi; // Induced dipole CR components for atom i. private double pxi; private double pyi; private double pzi; // Induced dipole + Induced dipole CR components for atom i. private double sxi; private double syi; private double szi; // Multipole components for atom k. private double qk; private double dxk; private double dyk; private double dzk; private double qxxk; private double qyyk; private double qzzk; private double qxyk; private double qxzk; private double qyzk; // Induced dipole components for atom k. private double uxk; private double uyk; private double uzk; // Induced dipole CR components for atom k. private double pxk; private double pyk; private double pzk; // Induced dipole + Induced dipole CR components for atom k. private double sxk; private double syk; private double szk; // Components of the potential, field and field gradient. private double E000; // Potential private double E100; // X Component of the Field private double E010; // Y Component of the Field private double E001; // Z Component of the Field private double E200; // XX Component of the Field Gradient private double E020; // YY Component of the Field Gradient private double E002; // ZZ Component of the Field Gradient private double E110; // XY Component of the Field Gradient private double E101; // XZ Component of the Field Gradient private double E011; // YZ Component of the Field Gradient // Cartesian tensor elements (for 1/R, erfc(Beta*R)/R or Thole damping. // l + m + n = 0 (1) private double R000; // l + m + n = 1 (3) 4 private double R100; private double R010; private double R001; // l + m + n = 2 (6) 10 private double R200; private double R020; private double R002; private double R110; private double R101; private double R011; // l + m + n = 3 (10) 20 private double R300; private double R030; private double R003; private double R210; private double R201; private double R120; private double R021; private double R102; private double R012; private double R111; // l + m + n = 4 (15) 35 private double R400; private double R040; private double R004; private double R310; private double R301; private double R130; private double R031; private double R103; private double R013; private double R220; private double R202; private double R022; private double R211; private double R121; private double R112; // l + m + n = 5 (21) 56 private double R500; private double R050; private double R005; private double R410; private double R401; private double R140; private double R041; private double R104; private double R014; private double R320; private double R302; private double R230; private double R032; private double R203; private double R023; private double R311; private double R131; private double R113; private double R221; private double R212; private double R122; // l + m + n = 6 (28) 84 private double R006; private double R402; private double R042; private double R204; private double R024; private double R222; private double R600; private double R060; private double R510; private double R501; private double R150; private double R051; private double R105; private double R015; private double R420; private double R240; private double R411; private double R141; private double R114; private double R330; private double R303; private double R033; private double R321; private double R231; private double R213; private double R312; private double R132; private double R123; // l + m + n = 0 (1) public final int t000; // l + m + n = 1 (3) 4 public final int t100; public final int t010; public final int t001; // l + m + n = 2 (6) 10 public final int t200; public final int t020; public final int t002; public final int t110; public final int t101; public final int t011; // l + m + n = 3 (10) 20 public final int t300; public final int t030; public final int t003; public final int t210; public final int t201; public final int t120; public final int t021; public final int t102; public final int t012; public final int t111; // l + m + n = 4 (15) 35 public final int t400; public final int t040; public final int t004; public final int t310; public final int t301; public final int t130; public final int t031; public final int t103; public final int t013; public final int t220; public final int t202; public final int t022; public final int t211; public final int t121; public final int t112; // l + m + n = 5 (21) 56 public final int t500; public final int t050; public final int t005; public final int t410; public final int t401; public final int t140; public final int t041; public final int t104; public final int t014; public final int t320; public final int t302; public final int t230; public final int t032; public final int t203; public final int t023; public final int t311; public final int t131; public final int t113; public final int t221; public final int t212; public final int t122; // l + m + n = 6 (28) 84 public final int t600; public final int t060; public final int t006; public final int t510; public final int t501; public final int t150; public final int t051; public final int t105; public final int t015; public final int t420; public final int t402; public final int t240; public final int t042; public final int t204; public final int t024; public final int t411; public final int t141; public final int t114; public final int t330; public final int t303; public final int t033; public final int t321; public final int t231; public final int t213; public final int t312; public final int t132; public final int t123; public final int t222; }