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.bonded;
39
40 import static ffx.numerics.math.DoubleMath.X;
41 import static ffx.numerics.math.DoubleMath.angle;
42 import static ffx.numerics.math.DoubleMath.length;
43 import static ffx.numerics.math.DoubleMath.normalize;
44 import static ffx.numerics.math.DoubleMath.scale;
45 import static ffx.numerics.math.DoubleMath.sub;
46 import static java.lang.String.format;
47 import static org.apache.commons.math3.util.FastMath.max;
48 import static org.apache.commons.math3.util.FastMath.min;
49
50 import ffx.numerics.atomic.AtomicDoubleArray3D;
51 import ffx.numerics.math.Double3;
52 import ffx.numerics.math.DoubleMath;
53 import ffx.potential.bonded.RendererCache.ViewModel;
54 import ffx.potential.parameters.AtomType;
55 import ffx.potential.parameters.BondType;
56 import ffx.potential.parameters.ForceField;
57
58 import java.io.Serial;
59 import java.util.ArrayList;
60 import java.util.List;
61 import java.util.logging.Logger;
62
63 import org.jogamp.java3d.BranchGroup;
64 import org.jogamp.java3d.Geometry;
65 import org.jogamp.java3d.LineArray;
66 import org.jogamp.java3d.Shape3D;
67 import org.jogamp.java3d.Transform3D;
68 import org.jogamp.java3d.TransformGroup;
69 import org.jogamp.vecmath.AxisAngle4d;
70 import org.jogamp.vecmath.Vector3d;
71
72
73
74
75
76
77
78 @SuppressWarnings("CloneableImplementsClone")
79 public class Bond extends BondedTerm {
80
81 @Serial
82 private static final long serialVersionUID = 1L;
83
84
85
86
87
88 static final float BUFF = 0.7f;
89
90 private static final Logger logger = Logger.getLogger(Bond.class.getName());
91
92 private static final float[] a0col = {
93 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f
94 };
95 private static final float[] f4a = {0.0f, 0.0f, 0.0f, 0.9f};
96 private static final float[] f4b = {0.0f, 0.0f, 0.0f, 0.9f};
97 private static final float[] f16 = {
98 0.0f, 0.0f, 0.0f, 0.9f, 0.0f, 0.0f, 0.0f, 0.9f, 0.0f, 0.0f, 0.0f, 0.9f, 0.0f, 0.0f, 0.0f, 0.9f
99 };
100 private static final double[] a13d = new double[3];
101 private static final double[] a23d = new double[3];
102 private static final double[] mid = new double[3];
103 private static final double[] diff3d = new double[3];
104 private static final double[] sum3d = new double[3];
105 private static final double[] coord = new double[12];
106 private static final double[] y = {0.0d, 1.0d, 0.0d};
107 private static final AxisAngle4d axisAngle = new AxisAngle4d();
108 private static final double[] bcross = new double[4];
109 private static final double[] cstart = new double[3];
110 private static final Vector3d pos3d = new Vector3d();
111
112
113
114 private final ArrayList<Bond> formsAngleWith = new ArrayList<>();
115
116
117
118 public BondType bondType = null;
119
120
121
122 private double rigidScale = 1.0;
123
124 private RendererCache.ViewModel viewModel = RendererCache.ViewModel.INVISIBLE;
125 private BranchGroup branchGroup;
126 private TransformGroup cy1tg, cy2tg;
127 private Transform3D cy1t3d, cy2t3d;
128 private Shape3D cy1, cy2;
129 private Vector3d scale;
130 private int detail = 3;
131 private LineArray la;
132 private int lineIndex;
133 private boolean wireVisible = true;
134
135
136
137
138
139
140 public Bond(String n) {
141 super(n);
142 }
143
144
145
146
147
148
149
150 public Bond(Atom a1, Atom a2) {
151 atoms = new Atom[2];
152 int i1 = a1.getIndex();
153 int i2 = a2.getIndex();
154 if (i1 < i2) {
155 atoms[0] = a1;
156 atoms[1] = a2;
157 } else {
158 atoms[0] = a2;
159 atoms[1] = a1;
160 }
161 setID_Key(false);
162 viewModel = RendererCache.ViewModel.WIREFRAME;
163 a1.setBond(this);
164 a2.setBond(this);
165 }
166
167
168
169
170
171
172
173
174 public static void logNoBondType(Atom a1, Atom a2, ForceField forceField) {
175 AtomType atomType1 = a1.getAtomType();
176 AtomType atomType2 = a2.getAtomType();
177 int[] c = {atomType1.atomClass, atomType2.atomClass};
178 String key = BondType.sortKey(c);
179 StringBuilder sb = new StringBuilder(
180 format(" No BondType for key: %s\n %s -> %s\n %s -> %s", key,
181 a1, atomType1, a2, atomType2));
182 int c1 = atomType1.atomClass;
183 int c2 = atomType2.atomClass;
184 List<AtomType> types1 = forceField.getSimilarAtomTypes(atomType1);
185 List<AtomType> types2 = forceField.getSimilarAtomTypes(atomType2);
186 List<BondType> bondTypes = new ArrayList<>();
187 boolean match = false;
188 for (AtomType type1 : types1) {
189 for (AtomType type2 : types2) {
190
191 if ((type1.atomClass != c1) && (type1.atomClass != c2) &&
192 (type2.atomClass != c1) && (type2.atomClass != c2)) {
193 continue;
194 }
195 BondType bondType = forceField.getBondType(type1, type2);
196 if (bondType != null && !bondTypes.contains(bondType)) {
197 if (!match) {
198 match = true;
199 sb.append("\n Similar Bond Types:");
200 }
201 bondTypes.add(bondType);
202 sb.append(format("\n %s", bondType));
203 }
204 }
205 }
206
207 logger.severe(sb.toString());
208 }
209
210
211
212
213 @Override
214 public int compareTo(BondedTerm b) {
215 if (b == null) {
216 throw new NullPointerException();
217 }
218 if (b == this) {
219 return 0;
220 }
221 if (!b.getClass().isInstance(this)) {
222 return super.compareTo(b);
223 }
224 int this0 = atoms[0].getIndex();
225 int a0 = b.atoms[0].getIndex();
226 if (this0 < a0) {
227 return -1;
228 }
229 if (this0 > a0) {
230 return 1;
231 }
232 int this1 = atoms[1].getIndex();
233 int a1 = b.atoms[1].getIndex();
234
235 return Integer.compare(this1, a1);
236 }
237
238
239
240
241
242
243 @Override
244 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
245 value = 0.0;
246 energy = 0.0;
247
248 if (!getUse()) {
249 return energy;
250 }
251 var atomA = atoms[0];
252 var atomB = atoms[1];
253 var va = atomA.getXYZ();
254 var vb = atomB.getXYZ();
255 var vab = va.sub(vb);
256 value = vab.length();
257 var prefactor = bondType.bondUnit * rigidScale * bondType.forceConstant;
258
259 var dv = value - bondType.distance;
260 if (bondType.bondFunction.hasFlatBottom()) {
261 if (dv > 0) {
262 dv = max(0, dv - bondType.flatBottomRadius);
263 } else if (dv < 0) {
264 dv = min(0, dv + bondType.flatBottomRadius);
265 }
266 }
267 var dv2 = dv * dv;
268
269 energy = prefactor * dv2 * (1.0 + bondType.cubic * dv + bondType.quartic * dv2);
270 if (gradient) {
271
272 var dedr = 2.0 * prefactor * dv * (1.0 + 1.5 * bondType.cubic * dv + 2.0 * bondType.quartic * dv2);
273 var de = 0.0;
274 if (value > 0.0) {
275 de = dedr / value;
276 }
277 var ga = vab.scale(de);
278 var ia = atomA.getIndex() - 1;
279 var ib = atomB.getIndex() - 1;
280 grad.add(threadID, ia, ga);
281 grad.sub(threadID, ib, ga);
282 }
283 value = dv;
284 return energy;
285 }
286
287
288
289
290
291
292
293
294
295 public Atom get1_2(Atom a) {
296 if (a == atoms[0]) {
297 return atoms[1];
298 }
299 if (a == atoms[1]) {
300 return atoms[0];
301 }
302 return null;
303 }
304
305
306
307
308
309
310 public double getCurrentDistance() {
311 double[] x1 = new double[3];
312 x1 = atoms[0].getXYZ(x1);
313 double[] x2 = new double[3];
314 x2 = atoms[1].getXYZ(x2);
315 return DoubleMath.dist(x1, x2);
316 }
317
318
319
320
321 public void log() {
322 logger.info(
323 format(
324 " %-8s %6d-%s %6d-%s %6.4f %6.4f %10.4f",
325 "Bond",
326 atoms[0].getIndex(),
327 atoms[0].getAtomType().name,
328 atoms[1].getIndex(),
329 atoms[1].getAtomType().name,
330 bondType.distance,
331 value,
332 energy));
333 }
334
335
336
337
338 @Override
339 public void removeFromParent() {
340 super.removeFromParent();
341 cy1 = null;
342 cy2 = null;
343 cy1tg = null;
344 cy2tg = null;
345 if (cy1t3d != null) {
346 RendererCache.poolTransform3D(cy1t3d);
347 RendererCache.poolTransform3D(cy2t3d);
348 cy1t3d = null;
349 cy2t3d = null;
350 }
351 if (branchGroup != null) {
352 branchGroup.detach();
353 branchGroup.setUserData(null);
354 RendererCache.poolDoubleCylinder(branchGroup);
355 branchGroup = null;
356 }
357 }
358
359
360
361
362
363
364 public void setBondType(BondType bondType) {
365 this.bondType = bondType;
366 }
367
368
369
370
371
372
373 public BondType getBondType() {
374 return bondType;
375 }
376
377
378
379
380
381
382 public void setColor(Atom a) {
383 if (viewModel != ViewModel.INVISIBLE
384 && viewModel != ViewModel.WIREFRAME
385 && branchGroup != null) {
386 if (a == atoms[0]) {
387 cy1.setAppearance(a.getAtomAppearance());
388 } else if (a == atoms[1]) {
389 cy2.setAppearance(a.getAtomAppearance());
390 }
391 }
392 setWireVisible(wireVisible);
393 }
394
395
396
397
398
399
400 public void setRigidScale(double rigidScale) {
401 this.rigidScale = rigidScale;
402 }
403
404
405
406
407
408
409 @Override
410 public void setView(RendererCache.ViewModel newViewModel, List<BranchGroup> newShapes) {
411 switch (newViewModel) {
412 case WIREFRAME -> {
413 viewModel = ViewModel.WIREFRAME;
414 setWireVisible(true);
415 setCylinderVisible(false, newShapes);
416 }
417 case SPACEFILL, INVISIBLE, RMIN -> {
418 viewModel = ViewModel.INVISIBLE;
419 setWireVisible(false);
420 setCylinderVisible(false, newShapes);
421 }
422 case RESTRICT -> {
423 if (!atoms[0].isSelected() || !atoms[1].isSelected()) {
424 viewModel = ViewModel.INVISIBLE;
425 setWireVisible(false);
426 setCylinderVisible(false, newShapes);
427 }
428 }
429 case BALLANDSTICK, TUBE -> {
430 viewModel = newViewModel;
431
432 double rad;
433 double len = getValue() / 2.0d;
434 if (viewModel == ViewModel.BALLANDSTICK) {
435 rad = 0.1d * RendererCache.radius;
436 } else {
437 rad = 0.2d * RendererCache.radius;
438 }
439 if (scale == null) {
440 scale = new Vector3d();
441 }
442 scale.set(rad, len, rad);
443 setWireVisible(false);
444 setCylinderVisible(true, newShapes);
445 }
446 case DETAIL -> {
447 int res = RendererCache.detail;
448 if (res != detail) {
449 detail = res;
450 if (branchGroup != null) {
451 Geometry geom1 = RendererCache.getCylinderGeom(0, detail);
452 Geometry geom2 = RendererCache.getCylinderGeom(1, detail);
453 Geometry geom3 = RendererCache.getCylinderGeom(2, detail);
454 cy1.removeAllGeometries();
455 cy2.removeAllGeometries();
456 cy1.addGeometry(geom1);
457 cy1.addGeometry(geom2);
458 cy1.addGeometry(geom3);
459 cy2.addGeometry(geom1);
460 cy2.addGeometry(geom2);
461 cy2.addGeometry(geom3);
462 }
463 }
464 if (scale == null) {
465 scale = new Vector3d();
466 }
467 double newRadius;
468 if (viewModel == ViewModel.BALLANDSTICK) {
469 newRadius = 0.1d * RendererCache.radius;
470 } else if (viewModel == ViewModel.TUBE) {
471 newRadius = 0.2d * RendererCache.radius;
472 } else {
473 break;
474 }
475 if (newRadius != scale.x) {
476 scale.x = newRadius;
477 scale.y = newRadius;
478 if (branchGroup != null) {
479 setView(viewModel, newShapes);
480 }
481 }
482 }
483 case SHOWHYDROGEN -> {
484 if (atoms[0].getAtomicNumber() == 1 || atoms[1].getAtomicNumber() == 1) {
485 setView(viewModel, newShapes);
486 }
487 }
488 case HIDEHYDROGEN -> {
489 if (atoms[0].getAtomicNumber() == 1 || atoms[1].getAtomicNumber() == 1) {
490 viewModel = ViewModel.INVISIBLE;
491 setWireVisible(false);
492 setCylinderVisible(false, newShapes);
493 }
494 }
495 case FILL, POINTS, LINES -> {
496 if (branchGroup != null && viewModel != ViewModel.INVISIBLE) {
497 cy1.setAppearance(atoms[0].getAtomAppearance());
498 cy2.setAppearance(atoms[1].getAtomAppearance());
499 }
500 }
501 }
502 }
503
504
505
506
507
508
509
510 public void setWire(LineArray l, int i) {
511 la = l;
512 lineIndex = i;
513 }
514
515
516
517
518
519
520 @Override
521 public void update() {
522
523 atoms[0].getXYZ(a13d);
524 atoms[1].getXYZ(a23d);
525 sub(a13d, a23d, diff3d);
526 double d = length(diff3d);
527 setValue(d);
528 DoubleMath.add(a13d, a23d, sum3d);
529 scale(sum3d, 0.5d, mid);
530
531 if (la != null) {
532 for (int i = 0; i < 3; i++) {
533 coord[i] = a13d[i];
534 coord[3 + i] = mid[i];
535 coord[6 + i] = mid[i];
536 coord[9 + i] = a23d[i];
537 }
538 la.setCoordinates(lineIndex, coord);
539 }
540
541 if (branchGroup != null) {
542 normalize(diff3d, diff3d);
543 scale.y = d / 2.0d;
544 setBondTransform3d(cy1t3d, mid, diff3d, d, true);
545 scale(diff3d, -1.0d, diff3d);
546 setBondTransform3d(cy2t3d, mid, diff3d, d, false);
547 cy1tg.setTransform(cy1t3d);
548 cy2tg.setTransform(cy2t3d);
549 }
550 }
551
552
553
554
555
556
557
558 boolean formsAngleWith(Bond b) {
559 for (Bond bond : formsAngleWith) {
560 if (b == bond) {
561 return true;
562 }
563 }
564 return false;
565 }
566
567
568
569
570
571
572
573
574 Atom getCommonAtom(Bond b) {
575 if (b == this || b == null) {
576 return null;
577 }
578 if (b.atoms[0] == atoms[0]) {
579 return atoms[0];
580 }
581 if (b.atoms[0] == atoms[1]) {
582 return atoms[1];
583 }
584 if (b.atoms[1] == atoms[0]) {
585 return atoms[0];
586 }
587 if (b.atoms[1] == atoms[1]) {
588 return atoms[1];
589 }
590 return null;
591 }
592
593
594
595
596
597
598
599
600 Atom getOtherAtom(Bond b) {
601 if (b == this || b == null) {
602 return null;
603 }
604 if (b.atoms[0] == atoms[0]) {
605 return atoms[1];
606 }
607 if (b.atoms[0] == atoms[1]) {
608 return atoms[0];
609 }
610 if (b.atoms[1] == atoms[0]) {
611 return atoms[1];
612 }
613 if (b.atoms[1] == atoms[1]) {
614 return atoms[0];
615 }
616 return null;
617 }
618
619
620
621
622
623
624 boolean sameGroup() {
625 return atoms[0].getParent() == atoms[1].getParent();
626 }
627
628
629
630
631
632
633 void setAngleWith(Bond b) {
634 formsAngleWith.add(b);
635 }
636
637
638
639
640
641
642 private void initJ3D(List<BranchGroup> newShapes) {
643 detail = RendererCache.detail;
644 branchGroup = RendererCache.doubleCylinderFactory(atoms[0], atoms[1], detail);
645 cy1tg = (TransformGroup) branchGroup.getChild(0);
646 cy2tg = (TransformGroup) branchGroup.getChild(1);
647 cy1 = (Shape3D) cy1tg.getChild(0);
648 cy2 = (Shape3D) cy2tg.getChild(0);
649 newShapes.add(branchGroup);
650 cy1t3d = RendererCache.transform3DFactory();
651 cy2t3d = RendererCache.transform3DFactory();
652 update();
653 }
654
655
656
657
658
659
660
661
662
663
664 private void setBondTransform3d(
665 Transform3D t3d, double[] pos, double[] orient, double len, boolean newRot) {
666
667 if (newRot) {
668 double angle = angle(orient, y);
669 X(y, orient, bcross);
670 bcross[3] = angle - Math.PI;
671 axisAngle.set(bcross);
672 }
673
674
675 scale(orient, len / 4.0d, cstart);
676 DoubleMath.add(cstart, pos, cstart);
677 pos3d.set(cstart);
678 t3d.setTranslation(pos3d);
679 t3d.setRotation(axisAngle);
680 t3d.setScale(scale);
681 }
682
683
684
685
686
687
688
689 private void setCylinderVisible(boolean visible, List<BranchGroup> newShapes) {
690 if (!visible) {
691
692 if (branchGroup != null) {
693 cy1.setPickable(false);
694 cy1.setAppearance(RendererCache.nullAp);
695 cy2.setPickable(false);
696 cy2.setAppearance(RendererCache.nullAp);
697
698 }
699 } else if (branchGroup == null) {
700
701 initJ3D(newShapes);
702 } else {
703
704 cy1t3d.setScale(scale);
705 cy1tg.setTransform(cy1t3d);
706 cy2t3d.setScale(scale);
707 cy2tg.setTransform(cy2t3d);
708 cy1.setAppearance(atoms[0].getAtomAppearance());
709 cy2.setAppearance(atoms[1].getAtomAppearance());
710 }
711 }
712
713
714
715
716
717
718 private void setWireVisible(boolean visible) {
719 if (!visible) {
720 wireVisible = false;
721 la.setColors(lineIndex, a0col);
722 } else {
723 wireVisible = true;
724 float[] cols = f16;
725 float[] col1 = f4a;
726 float[] col2 = f4b;
727 atoms[0].getAtomColor().get(col1);
728 atoms[1].getAtomColor().get(col2);
729 for (int i = 0; i < 3; i++) {
730 cols[i] = col1[i];
731 cols[4 + i] = col1[i];
732 cols[8 + i] = col2[i];
733 cols[12 + i] = col2[i];
734 }
735 la.setColors(lineIndex, cols);
736 }
737 }
738 }