1 // ****************************************************************************** 2 // 3 // Title: Force Field X. 4 // Description: Force Field X - Software for Molecular Biophysics. 5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024. 6 // 7 // This file is part of Force Field X. 8 // 9 // Force Field X is free software; you can redistribute it and/or modify it 10 // under the terms of the GNU General Public License version 3 as published by 11 // the Free Software Foundation. 12 // 13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT 14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 16 // details. 17 // 18 // You should have received a copy of the GNU General Public License along with 19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple 20 // Place, Suite 330, Boston, MA 02111-1307 USA 21 // 22 // Linking this library statically or dynamically with other modules is making a 23 // combined work based on this library. Thus, the terms and conditions of the 24 // GNU General Public License cover the whole combination. 25 // 26 // As a special exception, the copyright holders of this library give you 27 // permission to link this library with independent modules to produce an 28 // executable, regardless of the license terms of these independent modules, and 29 // to copy and distribute the resulting executable under terms of your choice, 30 // provided that you also meet, for each linked independent module, the terms 31 // and conditions of the license of that module. An independent module is a 32 // module which is not derived from or based on this library. If you modify this 33 // library, you may extend this exception to your version of the library, but 34 // you are not obligated to do so. If you do not wish to do so, delete this 35 // exception statement from your version. 36 // 37 // ****************************************************************************** 38 package ffx.numerics.multipole; 39 40 import static java.lang.Math.fma; 41 42 /** 43 * The PolarizableMultipole class defines a polarizable multipole. 44 * 45 * @author Michael J. Schnieders 46 * @since 1.0 47 */ 48 public class PolarizableMultipole { 49 50 private static final double oneThird = 1.0 / 3.0; 51 private static final double twoThirds = 2.0 / 3.0; 52 53 /** 54 * Partial charge. 55 */ 56 protected double q; 57 /** 58 * Dipole x-component. 59 */ 60 protected double dx; 61 /** 62 * Dipole y-component. 63 */ 64 protected double dy; 65 /** 66 * Dipole z-component. 67 */ 68 protected double dz; 69 /** 70 * Quadrupole xx-component multiplied by 1/3. 71 */ 72 protected double qxx; 73 /** 74 * Quadrupole yy-component multiplied by 1/3. 75 */ 76 protected double qyy; 77 /** 78 * Quadrupole zz-component multiplied by 1/3. 79 */ 80 protected double qzz; 81 /** 82 * Quadrupole xy-component multiplied by 2/3. 83 */ 84 protected double qxy; 85 /** 86 * Quadrupole xz-component multiplied by 2/3. 87 */ 88 protected double qxz; 89 /** 90 * Quadrupole xz-component multiplied by 2/3. 91 */ 92 protected double qyz; 93 94 /** 95 * Array of permanent multipole moments. 96 */ 97 private final double[] M = new double[10]; 98 99 /** 100 * Induced dipole x-component. 101 */ 102 protected double ux; 103 /** 104 * Induced dipole y-component. 105 */ 106 protected double uy; 107 /** 108 * Induced dipole z-component. 109 */ 110 protected double uz; 111 /** 112 * Induced dipole chain rule x-component. 113 */ 114 protected double px; 115 /** 116 * Induced dipole chain rule y-component. 117 */ 118 protected double py; 119 /** 120 * Induced dipole chain rule z-component. 121 */ 122 protected double pz; 123 /** 124 * Averaged induced dipole + induced dipole chain-rule x-component: sx = 0.5 * (ux + px). 125 */ 126 protected double sx; 127 /** 128 * Averaged induced dipole + induced dipole chain-rule y-component: sy = 0.5 * (uy + py). 129 */ 130 protected double sy; 131 /** 132 * Averaged induced dipole + induced dipole chain-rule z-component: sz = 0.5 * (uz + pz). 133 */ 134 protected double sz; 135 136 /** 137 * PolarizableMultipole constructor with zero moments. 138 */ 139 public PolarizableMultipole() { 140 } 141 142 /** 143 * PolarizableMultipole constructor. 144 * 145 * @param Q Multipole Q[q, dx, dy, dz, qxx, qyy, qzz, qxy, qxz, qyz] 146 * @param u Induced dipole u[ux, uy, uz] 147 * @param uCR Induced dipole chain-rule uCR[ux, uy, uz] 148 */ 149 public PolarizableMultipole(double[] Q, double[] u, double[] uCR) { 150 setPermanentMultipole(Q); 151 setInducedDipole(u, uCR); 152 } 153 154 /** 155 * Set the permanent multipole. 156 * <p> 157 * Note that the quadrupole trace components are multiplied by 1/3 and the 158 * off-diagonal components are multiplied by 2/3. 159 * 160 * @param Q Multipole Q[q, dx, dy, dz, qxx, qyy, qzz, qxy, qxz, qyz] 161 * @param u Induced dipole u[ux, uy, uz] 162 * @param uCR Induced dipole chain-rule uCR[ux, uy, uz] 163 */ 164 public void set(double[] Q, double[] u, double[] uCR) { 165 setPermanentMultipole(Q); 166 setInducedDipole(u, uCR); 167 } 168 169 /** 170 * Set the permanent multipole. 171 * <p> 172 * Note that the quadrupole trace components are multiplied by 1/3 and the 173 * off-diagonal components are multiplied by 2/3. 174 * 175 * @param Q Multipole Q[q, dx, dy, dz, qxx, qyy, qzz, qxy, qxz, qyz] 176 */ 177 public final void setPermanentMultipole(double[] Q) { 178 q = Q[0]; 179 dx = Q[1]; 180 dy = Q[2]; 181 dz = Q[3]; 182 qxx = Q[4] * oneThird; 183 qyy = Q[5] * oneThird; 184 qzz = Q[6] * oneThird; 185 qxy = Q[7] * twoThirds; 186 qxz = Q[8] * twoThirds; 187 qyz = Q[9] * twoThirds; 188 189 M[0] = q; 190 M[1] = dx; 191 M[2] = dy; 192 M[3] = dz; 193 M[4] = qxx; 194 M[5] = qyy; 195 M[6] = qzz; 196 M[7] = qxy; 197 M[8] = qxz; 198 M[9] = qyz; 199 } 200 201 /** 202 * Set the induced dipole. 203 * 204 * @param u Induced dipole u[ux, uy, uz] 205 * @param uCR Induced dipole chain-rule uCR[ux, uy, uz] 206 */ 207 public final void setInducedDipole(double[] u, double[] uCR) { 208 ux = u[0]; 209 uy = u[1]; 210 uz = u[2]; 211 px = uCR[0]; 212 py = uCR[1]; 213 pz = uCR[2]; 214 sx = 0.5 * (ux + px); 215 sy = 0.5 * (uy + py); 216 sz = 0.5 * (uz + pz); 217 } 218 219 /** 220 * Clear the induced dipoles. 221 */ 222 public void clearInducedDipoles() { 223 ux = 0.0; 224 uy = 0.0; 225 uz = 0.0; 226 px = 0.0; 227 py = 0.0; 228 pz = 0.0; 229 sx = 0.0; 230 sy = 0.0; 231 sz = 0.0; 232 } 233 234 /** 235 * Compute the scaled and averaged induced dipole. 236 * 237 * @param scaleInduction Induction mask scale factor. 238 * @param scaleEnergy Energy mask scale factor. 239 */ 240 public final void applyMasks(double scaleInduction, double scaleEnergy) { 241 // [Ux, Uy, Uz] resulted from induction masking rules, and we now apply the energy mask. 242 // [Px, Py, Pz] resulted from energy masking rules, and we now apply the induction mask. 243 sx = 0.5 * (ux * scaleEnergy + px * scaleInduction); 244 sy = 0.5 * (uy * scaleEnergy + py * scaleInduction); 245 sz = 0.5 * (uz * scaleEnergy + pz * scaleInduction); 246 } 247 248 /** 249 * Contract this multipole with the potential and its derivatives. 250 * 251 * @return The permanent multipole energy. 252 */ 253 protected final double multipoleEnergy(double[] field) { 254 double total = 0.0; 255 for (int i = 0; i < 10; i++) { 256 total = fma(M[i], field[i], total); 257 } 258 return total; 259 } 260 }