1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.nonbonded.pme;
39
40 import static ffx.numerics.math.ScalarMath.binomial;
41 import static java.lang.String.format;
42 import static java.util.Arrays.fill;
43
44 import ffx.potential.parameters.ForceField;
45 import java.util.logging.Level;
46 import java.util.logging.Logger;
47 import org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
48 import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
49 import org.apache.commons.math3.optimization.general.LevenbergMarquardtOptimizer;
50
51 @SuppressWarnings("deprecation")
52 public class SCFPredictorParameters {
53
54 private static final Logger logger = Logger.getLogger(SCFPredictorParameters.class.getName());
55
56
57 public int predictorOrder;
58
59 public int predictorStartIndex;
60
61 public int predictorCount;
62
63 public double[][][][] predictorInducedDipole;
64
65 public double[][][][] predictorInducedDipoleCR;
66
67 private LeastSquaresPredictor leastSquaresPredictor;
68 public LevenbergMarquardtOptimizer leastSquaresOptimizer;
69
70 public final SCFPredictor scfPredictor;
71 private final int nAtoms;
72
73 public SCFPredictorParameters(SCFPredictor scfPredictor, int nAtoms) {
74 this.nAtoms = nAtoms;
75 this.scfPredictor = scfPredictor;
76 }
77
78
79 public void aspcPredictor(LambdaMode lambdaMode,
80 double[][][] inducedDipole, double[][][] inducedDipoleCR) {
81
82 if (predictorCount < 6) {
83 return;
84 }
85
86 int mode;
87 switch (lambdaMode) {
88 case CONDENSED_NO_LIGAND:
89 mode = 1;
90 break;
91 case VAPOR:
92 mode = 2;
93 break;
94 default:
95 mode = 0;
96 }
97
98 final double[] aspc = {
99 22.0 / 7.0, -55.0 / 14.0, 55.0 / 21.0, -22.0 / 21.0, 5.0 / 21.0, -1.0 / 42.0
100 };
101
102
103 int index = predictorStartIndex;
104
105
106 for (int k = 0; k < 6; k++) {
107
108
109 double c = aspc[k];
110 for (int i = 0; i < nAtoms; i++) {
111 for (int j = 0; j < 3; j++) {
112 inducedDipole[0][i][j] += c * predictorInducedDipole[mode][index][i][j];
113 inducedDipoleCR[0][i][j] += c * predictorInducedDipoleCR[mode][index][i][j];
114 }
115 }
116 index++;
117 if (index >= predictorOrder) {
118 index = 0;
119 }
120 }
121 }
122
123 public void init(ForceField forceField) {
124 predictorCount = 0;
125 int defaultOrder = 6;
126 predictorOrder = forceField.getInteger("SCF_PREDICTOR_ORDER", defaultOrder);
127 if (scfPredictor == SCFPredictor.LS) {
128 leastSquaresPredictor = new LeastSquaresPredictor();
129 double eps = 1.0e-4;
130 leastSquaresOptimizer =
131 new LevenbergMarquardtOptimizer(
132 new org.apache.commons.math3.optimization.SimpleVectorValueChecker(eps, eps));
133 } else if (scfPredictor == SCFPredictor.ASPC) {
134 predictorOrder = 6;
135 }
136 predictorStartIndex = 0;
137 }
138
139
140
141
142
143 public void leastSquaresPredictor(LambdaMode lambdaMode,
144 double[][][] inducedDipole, double[][][] inducedDipoleCR) {
145 if (predictorCount < 2) {
146 return;
147 }
148 try {
149
150
151
152
153
154 leastSquaresPredictor.lambdaMode = lambdaMode;
155 leastSquaresPredictor.updateJacobianAndTarget();
156 int maxEvals = 100;
157 fill(leastSquaresPredictor.initialSolution, 0.0);
158 leastSquaresPredictor.initialSolution[0] = 1.0;
159 org.apache.commons.math3.optimization.PointVectorValuePair optimum =
160 leastSquaresOptimizer.optimize(
161 maxEvals,
162 leastSquaresPredictor,
163 leastSquaresPredictor.calculateTarget(),
164 leastSquaresPredictor.weights,
165 leastSquaresPredictor.initialSolution);
166 double[] optimalValues = optimum.getPoint();
167 if (logger.isLoggable(Level.FINEST)) {
168 logger.finest(format("\n LS RMS: %10.6f", leastSquaresOptimizer.getRMS()));
169 logger.finest(format(" LS Iterations: %10d", leastSquaresOptimizer.getEvaluations()));
170 logger.finest(
171 format(" Jacobian Evals: %10d", leastSquaresOptimizer.getJacobianEvaluations()));
172 logger.finest(format(" Chi Square: %10.6f", leastSquaresOptimizer.getChiSquare()));
173 logger.finest(" LS Coefficients");
174 for (int i = 0; i < predictorOrder - 1; i++) {
175 logger.finest(format(" %2d %10.6f", i + 1, optimalValues[i]));
176 }
177 }
178
179 int mode;
180 switch (lambdaMode) {
181 case CONDENSED_NO_LIGAND:
182 mode = 1;
183 break;
184 case VAPOR:
185 mode = 2;
186 break;
187 default:
188 mode = 0;
189 }
190
191
192 int index = predictorStartIndex;
193
194
195
196 for (int k = 0; k < predictorOrder - 1; k++) {
197
198
199 double c = optimalValues[k];
200 for (int i = 0; i < nAtoms; i++) {
201 for (int j = 0; j < 3; j++) {
202 inducedDipole[0][i][j] += c * predictorInducedDipole[mode][index][i][j];
203 inducedDipoleCR[0][i][j] += c * predictorInducedDipoleCR[mode][index][i][j];
204 }
205 }
206 index++;
207 if (index >= predictorOrder) {
208 index = 0;
209 }
210 }
211 } catch (Exception e) {
212 logger.log(Level.WARNING, " Exception computing predictor coefficients", e);
213 }
214 }
215
216
217 public void polynomialPredictor(LambdaMode lambdaMode,
218 double[][][] inducedDipole, double[][][] inducedDipoleCR) {
219
220 if (predictorCount == 0) {
221 return;
222 }
223
224 int mode;
225 switch (lambdaMode) {
226 case CONDENSED_NO_LIGAND:
227 mode = 1;
228 break;
229 case VAPOR:
230 mode = 2;
231 break;
232 default:
233 mode = 0;
234 }
235
236
237 int n = predictorOrder;
238 if (predictorCount < predictorOrder) {
239 n = predictorCount;
240 }
241
242
243 int index = predictorStartIndex;
244
245
246 double sign = -1.0;
247
248
249 for (int k = 0; k < n; k++) {
250
251
252 sign *= -1.0;
253 double c = sign * binomial(n, k);
254 for (int i = 0; i < nAtoms; i++) {
255 for (int j = 0; j < 3; j++) {
256 inducedDipole[0][i][j] += c * predictorInducedDipole[mode][index][i][j];
257 inducedDipoleCR[0][i][j] += c * predictorInducedDipoleCR[mode][index][i][j];
258 }
259 }
260 index++;
261 if (index >= predictorOrder) {
262 index = 0;
263 }
264 }
265 }
266
267
268 public void saveMutualInducedDipoles(LambdaMode lambdaMode,
269 double[][][] inducedDipole, double[][][] inducedDipoleCR,
270 double[][] directDipole, double[][] directDipoleCR) {
271
272 int mode;
273 switch (lambdaMode) {
274 case CONDENSED_NO_LIGAND:
275 mode = 1;
276 break;
277 case VAPOR:
278 mode = 2;
279 break;
280 default:
281 mode = 0;
282 }
283
284
285 predictorStartIndex--;
286 if (predictorStartIndex < 0) {
287 predictorStartIndex = predictorOrder - 1;
288 }
289
290 if (predictorCount < predictorOrder) {
291 predictorCount++;
292 }
293
294 for (int i = 0; i < nAtoms; i++) {
295 for (int j = 0; j < 3; j++) {
296 predictorInducedDipole[mode][predictorStartIndex][i][j] =
297 inducedDipole[0][i][j] - directDipole[i][j];
298 predictorInducedDipoleCR[mode][predictorStartIndex][i][j] =
299 inducedDipoleCR[0][i][j] - directDipoleCR[i][j];
300 }
301 }
302 }
303
304 private class LeastSquaresPredictor implements DifferentiableMultivariateVectorFunction {
305
306 double[] weights;
307 double[] target;
308 double[] values;
309 double[][] jacobian;
310 double[] initialSolution;
311 private final MultivariateMatrixFunction multivariateMatrixFunction = this::jacobian;
312
313 public LambdaMode lambdaMode = null;
314
315 LeastSquaresPredictor() {
316 weights = new double[2 * nAtoms * 3];
317 target = new double[2 * nAtoms * 3];
318 values = new double[2 * nAtoms * 3];
319 jacobian = new double[2 * nAtoms * 3][predictorOrder - 1];
320 initialSolution = new double[predictorOrder - 1];
321 fill(weights, 1.0);
322 initialSolution[0] = 1.0;
323 }
324
325 @Override
326 public MultivariateMatrixFunction jacobian() {
327 return multivariateMatrixFunction;
328 }
329
330 @Override
331 public double[] value(double[] variables) {
332 int mode;
333 switch (lambdaMode) {
334 case CONDENSED_NO_LIGAND:
335 mode = 1;
336 break;
337 case VAPOR:
338 mode = 2;
339 break;
340 default:
341 mode = 0;
342 }
343
344 for (int i = 0; i < nAtoms; i++) {
345 int index = 6 * i;
346 values[index] = 0;
347 values[index + 1] = 0;
348 values[index + 2] = 0;
349 values[index + 3] = 0;
350 values[index + 4] = 0;
351 values[index + 5] = 0;
352 int pi = predictorStartIndex + 1;
353 if (pi >= predictorOrder) {
354 pi = 0;
355 }
356 for (int j = 0; j < predictorOrder - 1; j++) {
357 values[index] += variables[j] * predictorInducedDipole[mode][pi][i][0];
358 values[index + 1] += variables[j] * predictorInducedDipole[mode][pi][i][1];
359 values[index + 2] += variables[j] * predictorInducedDipole[mode][pi][i][2];
360 values[index + 3] += variables[j] * predictorInducedDipoleCR[mode][pi][i][0];
361 values[index + 4] += variables[j] * predictorInducedDipoleCR[mode][pi][i][1];
362 values[index + 5] += variables[j] * predictorInducedDipoleCR[mode][pi][i][2];
363 pi++;
364 if (pi >= predictorOrder) {
365 pi = 0;
366 }
367 }
368 }
369 return values;
370 }
371
372 double[] calculateTarget() {
373 return target;
374 }
375
376 void updateJacobianAndTarget() {
377 int mode;
378 switch (lambdaMode) {
379 case CONDENSED_NO_LIGAND:
380 mode = 1;
381 break;
382 case VAPOR:
383 mode = 2;
384 break;
385 default:
386 mode = 0;
387 }
388
389
390 int index = 0;
391 for (int i = 0; i < nAtoms; i++) {
392 target[index++] = predictorInducedDipole[mode][predictorStartIndex][i][0];
393 target[index++] = predictorInducedDipole[mode][predictorStartIndex][i][1];
394 target[index++] = predictorInducedDipole[mode][predictorStartIndex][i][2];
395 target[index++] = predictorInducedDipoleCR[mode][predictorStartIndex][i][0];
396 target[index++] = predictorInducedDipoleCR[mode][predictorStartIndex][i][1];
397 target[index++] = predictorInducedDipoleCR[mode][predictorStartIndex][i][2];
398 }
399
400
401 index = predictorStartIndex + 1;
402 if (index >= predictorOrder) {
403 index = 0;
404 }
405 for (int j = 0; j < predictorOrder - 1; j++) {
406 int ji = 0;
407 for (int i = 0; i < nAtoms; i++) {
408 jacobian[ji++][j] = predictorInducedDipole[mode][index][i][0];
409 jacobian[ji++][j] = predictorInducedDipole[mode][index][i][1];
410 jacobian[ji++][j] = predictorInducedDipole[mode][index][i][2];
411 jacobian[ji++][j] = predictorInducedDipoleCR[mode][index][i][0];
412 jacobian[ji++][j] = predictorInducedDipoleCR[mode][index][i][1];
413 jacobian[ji++][j] = predictorInducedDipoleCR[mode][index][i][2];
414 }
415 index++;
416 if (index >= predictorOrder) {
417 index = 0;
418 }
419 }
420 }
421
422 private double[][] jacobian(double[] variables) {
423 return jacobian;
424 }
425 }
426 }