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.algorithms.optimize.manybody;
39
40 import static java.lang.String.format;
41 import static java.util.Arrays.fill;
42 import static org.apache.commons.math3.util.FastMath.min;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44
45 import edu.rit.pj.ParallelTeam;
46 import ffx.algorithms.AlgorithmListener;
47 import ffx.algorithms.optimize.RotamerOptimization;
48 import ffx.crystal.Crystal;
49 import ffx.crystal.SymOp;
50 import ffx.potential.MolecularAssembly;
51 import ffx.potential.bonded.Atom;
52 import ffx.potential.bonded.MultiResidue;
53 import ffx.potential.bonded.Residue;
54 import ffx.potential.bonded.ResidueState;
55 import ffx.potential.bonded.Rotamer;
56 import ffx.potential.bonded.RotamerLibrary;
57 import ffx.potential.nonbonded.NeighborList;
58
59 import java.util.Arrays;
60 import java.util.List;
61 import java.util.logging.Level;
62 import java.util.logging.Logger;
63
64 import org.apache.commons.math3.util.CombinatoricsUtils;
65 import org.apache.commons.math3.util.FastMath;
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82 public class DistanceMatrix {
83
84 private static final Logger logger = Logger.getLogger(DistanceMatrix.class.getName());
85
86
87
88
89 private final MolecularAssembly molecularAssembly;
90
91
92
93 private final AlgorithmListener algorithmListener;
94
95
96
97
98 private final Residue[] allResiduesArray;
99
100
101
102
103 private final List<Residue> allResiduesList;
104
105
106
107 private final int numResidues;
108
109
110
111 private final RotamerOptimization.DistanceMethod distanceMethod;
112
113
114
115 private final double distance;
116
117
118
119 private final double twoBodyCutoffDist;
120
121
122
123 private final double threeBodyCutoffDist;
124
125
126
127
128 private final boolean lazyMatrix;
129
130
131
132
133
134
135 private double[][][][] distanceMatrix;
136
137 public DistanceMatrix(MolecularAssembly molecularAssembly, AlgorithmListener algorithmListener,
138 Residue[] allResiduesArray, List<Residue> allResiduesList,
139 RotamerOptimization.DistanceMethod distanceMethod, double distance, double twoBodyCutoffDist,
140 double threeBodyCutoffDist, boolean lazyMatrix) {
141 this.molecularAssembly = molecularAssembly;
142 this.algorithmListener = algorithmListener;
143 this.allResiduesArray = allResiduesArray;
144 this.allResiduesList = allResiduesList;
145 this.numResidues = allResiduesArray.length;
146 this.distanceMethod = distanceMethod;
147 this.distance = distance;
148 this.twoBodyCutoffDist = twoBodyCutoffDist;
149 this.threeBodyCutoffDist = threeBodyCutoffDist;
150 this.lazyMatrix = lazyMatrix;
151
152 distanceMatrix();
153 }
154
155
156
157
158
159
160
161 private static boolean areFiniteAndNotMax(double... values) {
162 return Arrays.stream(values)
163 .allMatch((double val) -> Double.isFinite(val) && val < Double.MAX_VALUE);
164 }
165
166
167
168
169
170
171
172
173
174
175
176 public double checkDistMatrix(int i, int ri, int j, int rj) {
177 if (i > j) {
178 int temp = i;
179 i = j;
180 j = temp;
181 temp = ri;
182 ri = rj;
183 rj = temp;
184 }
185 double dist = distanceMatrix[i][ri][j][rj];
186 if (dist < 0) {
187 dist = evaluateDistance(i, ri, j, rj);
188 distanceMatrix[i][ri][j][rj] = dist;
189 }
190 return dist;
191 }
192
193
194
195
196
197
198
199
200
201
202 public boolean checkPairDistThreshold(int i, int ri, int j, int rj) {
203 if (twoBodyCutoffDist <= 0 || !Double.isFinite(twoBodyCutoffDist)) {
204 return false;
205 }
206 return (get2BodyDistance(i, ri, j, rj) > twoBodyCutoffDist);
207 }
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223 public boolean checkQuadDistThreshold(int i, int ri, int j, int rj, int k, int rk, int l, int rl) {
224 if (checkTriDistThreshold(i, ri, j, rj, k, rk) || checkTriDistThreshold(i, ri, j, rj, l, rl)
225 || checkTriDistThreshold(i, ri, k, rk, l, rl) || checkTriDistThreshold(j, rj, k, rk, l,
226 rl)) {
227 return true;
228 }
229
230 if (threeBodyCutoffDist <= 0 || !Double.isFinite(threeBodyCutoffDist)) {
231 return false;
232 }
233 return get4BodyDistance(i, ri, j, rj, k, rk, l, rl) > threeBodyCutoffDist;
234 }
235
236
237
238
239
240
241
242
243
244
245
246
247
248 public boolean checkTriDistThreshold(int i, int ri, int j, int rj, int k, int rk) {
249 if (checkPairDistThreshold(i, ri, j, rj) || checkPairDistThreshold(i, ri, k, rk)
250 || checkPairDistThreshold(j, rj, k, rk)) {
251 return true;
252 }
253 if (threeBodyCutoffDist <= 0 || !Double.isFinite(threeBodyCutoffDist)) {
254 return false;
255 }
256 return (get3BodyDistance(i, ri, j, rj, k, rk) > threeBodyCutoffDist);
257 }
258
259
260
261
262
263
264
265
266
267
268
269
270
271 public double get2BodyDistance(int i, int ri, int j, int rj) {
272 if (i > j) {
273 int temp = i;
274 i = j;
275 j = temp;
276 temp = ri;
277 ri = rj;
278 rj = temp;
279 }
280
281 if (distanceMethod == RotamerOptimization.DistanceMethod.ROTAMER) {
282 return checkDistMatrix(i, ri, j, rj);
283 } else {
284 return getResidueDistance(i, ri, j, rj);
285 }
286 }
287
288
289
290
291
292
293
294
295
296
297
298
299
300 public double get3BodyResidueDistance(int i, int ri, int j, int rj, int k, int rk) {
301 double ij = getResidueDistance(i, ri, j, rj);
302 double ik = getResidueDistance(i, ri, k, rk);
303 double jk = getResidueDistance(j, rj, k, rk);
304 if (areFiniteAndNotMax(ij, ik, jk)) {
305 return sqrt((ij * ij + ik * ik + jk * jk) / 3.0);
306 } else {
307 return Double.MAX_VALUE;
308 }
309 }
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325 public double get4BodyResidueDistance(int i, int ri, int j, int rj, int k, int rk, int l, int rl) {
326 double ij = getResidueDistance(i, ri, j, rj);
327 double ik = getResidueDistance(i, ri, k, rk);
328 double il = getResidueDistance(i, ri, l, rl);
329 double jk = getResidueDistance(j, rj, k, rk);
330 double jl = getResidueDistance(j, rj, l, rl);
331 double kl = getResidueDistance(k, rk, l, rl);
332 if (areFiniteAndNotMax(ij, ik, il, jk, jl, kl)) {
333 return sqrt((ij * ij + ik * ik + il * il + jk * jk + jl * jl + kl * kl) / 6.0);
334 } else {
335 return Double.MAX_VALUE;
336 }
337 }
338
339
340
341
342
343
344
345 public double getRawNBodyDistance(int... resRot) {
346 int nRes = resRot.length;
347 if (nRes % 2 != 0) {
348 throw new IllegalArgumentException(" Must have an even number of arguments; res-rot pairs!");
349 }
350 nRes /= 2;
351 if (nRes < 2) {
352 throw new IllegalArgumentException(" Must have >= 4 arguments; at least 2 res-rot pairs!");
353 }
354
355 int numDists = (int) CombinatoricsUtils.binomialCoefficient(nRes, 2);
356 double mult = 1.0 / numDists;
357 double totDist2 = 0.0;
358
359 for (int i = 0; i < nRes - 1; i++) {
360 int i2 = 2 * i;
361 for (int j = i + 1; j < nRes; j++) {
362 int j2 = 2 * j;
363 double rawDist = checkDistMatrix(resRot[i2], resRot[i2 + 1], resRot[j2], resRot[j2 + 1]);
364 if (!Double.isFinite(rawDist) || rawDist == Double.MAX_VALUE) {
365 return Double.MAX_VALUE;
366 }
367 totDist2 += (rawDist * mult);
368 }
369 }
370 return FastMath.sqrt(totDist2);
371 }
372
373
374
375
376
377
378
379
380
381
382
383 public double getResidueDistance(int i, int ri, int j, int rj) {
384 if (i > j) {
385 int temp = i;
386 i = j;
387 j = temp;
388 temp = ri;
389 ri = rj;
390 rj = temp;
391 }
392
393 double minDist = Double.MAX_VALUE;
394 final int lenri = distanceMatrix[i].length;
395 final int lenrj = distanceMatrix[i][ri][j].length;
396 for (int roti = 0; roti < lenri; roti++) {
397 for (int rotj = 0; rotj < lenrj; rotj++) {
398 minDist = Math.min(minDist, checkDistMatrix(i, roti, j, rotj));
399 }
400 }
401 return minDist;
402 }
403
404
405
406
407
408
409
410
411
412 public double interResidueDistance(double[][] resi, double[][] resj, SymOp symOp) {
413 double dist = Double.MAX_VALUE;
414 Crystal crystal = molecularAssembly.getCrystal();
415 for (double[] xi : resi) {
416 for (double[] xj : resj) {
417 if (symOp != null) {
418 crystal.applySymOp(xj, xj, symOp);
419 }
420 double r = crystal.image(xi[0] - xj[0], xi[1] - xj[1], xi[2] - xj[2]);
421 if (r < dist) {
422 dist = r;
423 }
424 }
425 }
426 return sqrt(dist);
427 }
428
429 private void distanceMatrix() {
430
431 distanceMatrix = new double[numResidues - 1][][][];
432 long numDistances = 0L;
433 for (int i = 0; i < (numResidues - 1); i++) {
434 Residue residuei = allResiduesArray[i];
435 int lengthRi;
436 try {
437 lengthRi = residuei.getRotamers().length;
438 } catch (IndexOutOfBoundsException ex) {
439 logger.warning(format(" Residue i %s has null rotamers.", residuei.toFormattedString(false, true)));
440 continue;
441 }
442 distanceMatrix[i] = new double[lengthRi][][];
443 for (int ri = 0; ri < lengthRi; ri++) {
444 distanceMatrix[i][ri] = new double[numResidues][];
445 for (int j = (i + 1); j < numResidues; j++) {
446 Residue residuej = allResiduesArray[j];
447 int lengthRj;
448 try {
449 lengthRj = residuej.getRotamers().length;
450 } catch (IndexOutOfBoundsException ex) {
451 logger.warning(format(" Residue j %s has null rotamers.", residuej.toFormattedString(false, true)));
452 continue;
453 }
454 distanceMatrix[i][ri][j] = new double[lengthRj];
455 numDistances += lengthRj;
456 if (!lazyMatrix) {
457 fill(distanceMatrix[i][ri][j], Double.MAX_VALUE);
458 } else {
459 fill(distanceMatrix[i][ri][j], -1.0);
460 }
461 }
462 }
463 }
464
465 logger.info(format(" Number of pairwise distances: %d", numDistances));
466
467 if (!lazyMatrix) {
468 ResidueState[] orig = ResidueState.storeAllCoordinates(allResiduesList);
469 int nMultiRes = 0;
470
471
472
473 Atom[] atoms = new Atom[numResidues];
474 for (int i = 0; i < numResidues; i++) {
475 Residue residuei = allResiduesArray[i];
476 atoms[i] = residuei.getReferenceAtom();
477 if (residuei instanceof MultiResidue) {
478 ++nMultiRes;
479 }
480 }
481
482
483
484
485
486
487
488 ParallelTeam parallelTeam = getParallelTeam(nMultiRes);
489 Crystal crystal = molecularAssembly.getCrystal();
490 int nSymm = crystal.spaceGroup.getNumberOfSymOps();
491 logger.info("\n Computing Residue Distance Matrix");
492
493 double nlistCutoff = Math.max(Math.max(distance, twoBodyCutoffDist), threeBodyCutoffDist);
494
495
496
497 double conservativeBuffer = 25.0;
498 nlistCutoff += conservativeBuffer;
499
500
501
502
503
504
505 if (!crystal.aperiodic()) {
506 double sphere = min(min(crystal.interfacialRadiusA, crystal.interfacialRadiusB),
507 crystal.interfacialRadiusC);
508 if (nlistCutoff > sphere) {
509 nlistCutoff = sphere;
510 }
511 }
512
513 NeighborList neighborList = new NeighborList(null, crystal, atoms, nlistCutoff, 0.0,
514 parallelTeam);
515
516
517 double[][] xyz = new double[nSymm][3 * numResidues];
518 double[] in = new double[3];
519 double[] out = new double[3];
520 for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
521 SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
522 for (int i = 0; i < numResidues; i++) {
523 Atom atom = atoms[i];
524 in[0] = atom.getX();
525 in[1] = atom.getY();
526 in[2] = atom.getZ();
527 crystal.applySymOp(in, out, symOp);
528 int iX = i * 3;
529 int iY = iX + 1;
530 int iZ = iX + 2;
531 xyz[iSymOp][iX] = out[0];
532 xyz[iSymOp][iY] = out[1];
533 xyz[iSymOp][iZ] = out[2];
534 }
535 }
536
537
538 int[][][] lists = new int[nSymm][numResidues][];
539 boolean[] use = new boolean[numResidues];
540 fill(use, true);
541 boolean forceRebuild = true;
542 boolean printLists = false;
543 long neighborTime = -System.nanoTime();
544 neighborList.buildList(xyz, lists, use, forceRebuild, printLists);
545
546 neighborTime += System.nanoTime();
547 logger.info(format(" Built residue neighbor list: %8.3f sec", neighborTime * 1.0e-9));
548
549 DistanceRegion distanceRegion = new DistanceRegion(parallelTeam.getThreadCount(), numResidues,
550 crystal, lists, neighborList.getPairwiseSchedule());
551
552 long parallelTime = -System.nanoTime();
553 try {
554 distanceRegion.init(this, molecularAssembly, allResiduesArray, algorithmListener,
555 distanceMatrix);
556 parallelTeam.execute(distanceRegion);
557 } catch (Exception e) {
558 String message = " Exception compting residue distance matrix.";
559 logger.log(Level.SEVERE, message, e);
560 }
561 parallelTime += System.nanoTime();
562 logger.info(format(" Pairwise distance matrix: %8.3f sec", parallelTime * 1.0e-9));
563
564 ResidueState.revertAllCoordinates(allResiduesList, orig);
565 try {
566 parallelTeam.shutdown();
567 } catch (Exception ex) {
568 logger.warning(format(" Exception shutting down parallel team for the distance matrix: %s", ex));
569 }
570 }
571 }
572
573 private ParallelTeam getParallelTeam(int nMultiRes) {
574 int nThreads;
575 if (molecularAssembly.getPotentialEnergy().getParallelTeam() != null) {
576 nThreads = (nMultiRes > 1) ? 1
577 : molecularAssembly.getPotentialEnergy().getParallelTeam().getThreadCount();
578 } else {
579
580 nThreads = 16;
581 }
582 ParallelTeam parallelTeam = new ParallelTeam(nThreads);
583 return parallelTeam;
584 }
585
586
587
588
589
590
591
592
593
594
595
596
597
598 private double get3BodyDistance(int i, int ri, int j, int rj, int k, int rk) {
599 double ij = get2BodyDistance(i, ri, j, rj);
600 double ik = get2BodyDistance(i, ri, k, rk);
601 double jk = get2BodyDistance(j, rj, k, rk);
602 if (areFiniteAndNotMax(ij, ik, jk)) {
603 return sqrt((ij * ij + ik * ik + jk * jk) / 3.0);
604 } else {
605 return Double.MAX_VALUE;
606 }
607 }
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623 private double get4BodyDistance(int i, int ri, int j, int rj, int k, int rk, int l, int rl) {
624 double ij = get2BodyDistance(i, ri, j, rj);
625 double ik = get2BodyDistance(i, ri, k, rk);
626 double il = get2BodyDistance(i, ri, l, rl);
627 double jk = get2BodyDistance(j, rj, k, rk);
628 double jl = get2BodyDistance(j, rj, l, rl);
629 double kl = get2BodyDistance(k, rk, l, rl);
630 if (areFiniteAndNotMax(ij, ik, il, jk, jl, kl)) {
631 return sqrt((ij * ij + ik * ik + il * il + jk * jk + jl * jl + kl * kl) / 6.0);
632 } else {
633 return Double.MAX_VALUE;
634 }
635 }
636
637
638
639
640
641
642
643
644
645
646
647 private double evaluateDistance(int i, int ri, int j, int rj) {
648 Residue resi = allResiduesArray[i];
649 Rotamer[] rotamersI = resi.getRotamers();
650 Rotamer roti = rotamersI[ri];
651 double[][] xi;
652 if (roti.equals(resi.getRotamer())) {
653 xi = resi.storeCoordinateArray();
654 } else {
655 ResidueState origI = resi.storeState();
656 RotamerLibrary.applyRotamer(resi, roti);
657 xi = resi.storeCoordinateArray();
658 resi.revertState(origI);
659 }
660
661 Residue resj = allResiduesArray[j];
662 Rotamer[] rotamersJ = resj.getRotamers();
663 Rotamer rotj = rotamersJ[rj];
664 double[][] xj;
665 if (rotj.equals(resj.getRotamer())) {
666 xj = resj.storeCoordinateArray();
667 } else {
668 ResidueState origJ = resj.storeState();
669 RotamerLibrary.applyRotamer(resj, rotj);
670 xj = resj.storeCoordinateArray();
671 resj.revertState(origJ);
672 }
673
674 Crystal crystal = molecularAssembly.getCrystal();
675 int nSymm = crystal.spaceGroup.getNumberOfSymOps();
676 double minDist = Double.MAX_VALUE;
677 for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
678 SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
679 double dist = interResidueDistance(xi, xj, symOp);
680 minDist = Math.min(dist, minDist);
681 }
682 return minDist;
683 }
684 }