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.length;
41 import static ffx.numerics.math.DoubleMath.normalize;
42 import static ffx.numerics.math.DoubleMath.scale;
43 import static ffx.numerics.math.DoubleMath.sub;
44 import static ffx.potential.bonded.Bond.logNoBondType;
45 import static ffx.potential.bonded.NamingUtils.nameAcetylCap;
46 import static java.lang.String.format;
47 import static java.lang.System.arraycopy;
48 import static org.apache.commons.math3.util.FastMath.abs;
49 import static org.apache.commons.math3.util.FastMath.cos;
50 import static org.apache.commons.math3.util.FastMath.max;
51 import static org.apache.commons.math3.util.FastMath.sin;
52 import static org.apache.commons.math3.util.FastMath.sqrt;
53 import static org.apache.commons.math3.util.FastMath.toRadians;
54
55 import ffx.numerics.math.DoubleMath;
56 import ffx.potential.MolecularAssembly;
57 import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
58 import ffx.potential.bonded.AminoAcidUtils.SideChainType;
59 import ffx.potential.parameters.AtomType;
60 import ffx.potential.parameters.BioType;
61 import ffx.potential.parameters.BondType;
62 import ffx.potential.parameters.ForceField;
63 import ffx.utilities.StringUtils;
64
65 import java.io.Serial;
66 import java.util.ArrayList;
67 import java.util.Arrays;
68 import java.util.Comparator;
69 import java.util.List;
70 import java.util.Optional;
71 import java.util.logging.Level;
72 import java.util.logging.Logger;
73 import java.util.stream.Collectors;
74
75
76
77
78
79
80
81 public class BondedUtils {
82
83 private static final Logger logger = Logger.getLogger(BondedUtils.class.getName());
84 private static final double eps = 0.0000001d;
85
86
87
88
89
90
91
92
93 public static boolean atomAttachedToAtom(Atom a1, Atom a2) {
94 assert a1 != a2;
95 return a1.getBonds().stream().anyMatch((Bond b) -> b.get1_2(a1) == a2);
96 }
97
98
99
100
101
102
103
104
105
106
107 public static Bond buildBond(Atom a1, Atom a2, ForceField forceField, List<Bond> bondList) {
108 Bond bond = new Bond(a1, a2);
109 BondType bondType = forceField.getBondType(a1.getAtomType(), a2.getAtomType());
110 if (bondType == null) {
111 logNoBondType(a1, a2, forceField);
112 } else {
113 bond.setBondType(bondType);
114 }
115 if (bondList != null) {
116 bondList.add(bond);
117 }
118 return bond;
119 }
120
121
122
123
124
125
126
127
128
129
130
131
132
133 public static Atom buildHeavy(MSGroup residue, String atomName, Atom bondedTo,
134 int key, ForceField forceField, List<Bond> bondList)
135 throws MissingHeavyAtomException, MissingAtomTypeException {
136 Atom atom = (Atom) residue.getAtomNode(atomName);
137 AtomType atomType = findAtomType(key, forceField);
138 if (atomType == null) {
139 Residue res = (Residue) residue;
140 throw new MissingAtomTypeException(res, atom);
141 }
142 if (atom == null) {
143 throw new MissingHeavyAtomException(atomName, atomType, bondedTo);
144 }
145 atom.setAtomType(atomType);
146 if (bondedTo != null) {
147 buildBond(atom, bondedTo, forceField, bondList);
148 }
149 return atom;
150 }
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168 public static Atom buildHeavy(MSGroup residue, String atomName, Atom ia, double bond,
169 Atom ib, double angle1, Atom ic, double angle2, int chiral, int lookUp,
170 ForceField forceField) {
171 AtomType atomType = findAtomType(lookUp, forceField);
172 return buildHeavyAtom(residue, atomName, ia, bond, ib, angle1, ic, angle2, chiral, atomType);
173 }
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192 public static Atom buildHeavy(MSGroup residue, String atomName, Atom ia, double bond, Atom ib,
193 double angle1, Atom ic, double angle2, int chiral, int lookUp, ForceField forceField,
194 List<Bond> bondList) {
195 AtomType atomType = findAtomType(lookUp, forceField);
196 return buildHeavyAtom(residue, atomName, ia, bond, ib, angle1, ic, angle2, chiral, atomType,
197 forceField, bondList);
198 }
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216 public static Atom buildHeavy(MSGroup residue, SideChainType atomName, Atom ia, double bond,
217 Atom ib, double angle1, Atom ic, double angle2, int chiral, ForceField forceField,
218 List<Bond> bondList) {
219 AtomType atomType = findAtomType(atomName.getType(), forceField);
220 return buildHeavyAtom(residue, atomName.name(), ia, bond, ib, angle1, ic, angle2, chiral,
221 atomType, forceField, bondList);
222 }
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241 public static Atom buildH(MSGroup residue, String atomName, Atom ia, double bond, Atom ib,
242 double angle1, Atom ic, double angle2, int chiral, int lookUp, ForceField forceField,
243 List<Bond> bondList) {
244 AtomType atomType = findAtomType(lookUp, forceField);
245 return buildHydrogenAtom(residue, atomName, ia, bond, ib, angle1, ic, angle2, chiral,
246 atomType, forceField, bondList);
247 }
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265 public static Atom buildH(MSGroup residue, SideChainType atomName, Atom ia, double bond,
266 Atom ib, double angle1, Atom ic, double angle2, int chiral, ForceField forceField,
267 List<Bond> bondList) {
268 AtomType atomType = findAtomType(atomName.getType(), forceField);
269 return buildHydrogenAtom(residue, atomName.name(), ia, bond, ib, angle1, ic, angle2,
270 chiral, atomType, forceField, bondList);
271 }
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290 public static Atom buildHydrogenAtom(MSGroup residue, String atomName, Atom ia, double bond,
291 Atom ib, double angle1, Atom ic, double angle2, int chiral, AtomType atomType,
292 ForceField forceField, List<Bond> bondList) {
293 if (atomType == null) {
294 return null;
295 }
296 Atom atom = (Atom) residue.getAtomNode(atomName);
297
298 if (atom == null) {
299 String dAtomName = atomName.replaceFirst("H", "D");
300 atom = (Atom) residue.getAtomNode(dAtomName);
301 }
302
303
304 if (atom == null && residue instanceof Molecule molecule) {
305 if (StringUtils.looksLikeWater(molecule.getName())) {
306 atom = (Atom) molecule.getAtomNode("H");
307 if (atom == null) {
308 atom = (Atom) molecule.getAtomNode("D");
309 }
310 if (atom != null) {
311 atom.setName(atomName);
312 }
313 }
314 }
315
316 if (atom == null) {
317 String resName = ia.getResidueName();
318 int resSeq = ia.getResidueNumber();
319 Character chainID = ia.getChainID();
320 Character altLoc = ia.getAltLoc();
321 String segID = ia.getSegID();
322 double occupancy = ia.getOccupancy();
323 double tempFactor = ia.getTempFactor();
324 atom = new Atom(0, atomName, altLoc, new double[3], resName, resSeq, chainID,
325 occupancy, tempFactor, segID, true);
326 residue.addMSNode(atom);
327 intxyz(atom, ia, bond, ib, angle1, ic, angle2, chiral);
328 }
329 atom.setAtomType(atomType);
330 buildBond(ia, atom, forceField, bondList);
331 return atom;
332 }
333
334
335
336
337
338
339
340
341 public static AtomType findAtomType(int key, ForceField forceField) {
342 BioType bioType = forceField.getBioType(Integer.toString(key));
343 if (bioType != null) {
344 AtomType atomType = forceField.getAtomType(Integer.toString(bioType.atomType));
345 if (atomType != null) {
346 return atomType;
347 } else {
348 logger.severe(format("The atom type %s was not found for biotype %s.",
349 bioType.atomType, bioType));
350 }
351 }
352 return null;
353 }
354
355
356
357
358
359
360
361
362 public static List<Atom> findAtomsOfElement(Residue residue, int element) {
363 return residue.getAtomList().stream()
364 .filter((Atom a) -> a.getAtomType().atomicNumber == element)
365 .collect(Collectors.toList());
366 }
367
368
369
370
371
372
373
374
375 public static List<Atom> findBondedAtoms(Atom atom, int element) {
376 return findBondedAtoms(atom, null, element);
377 }
378
379
380
381
382
383
384
385
386
387
388 public static List<Atom> findBondedAtoms(Atom atom, Atom toExclude, int element) {
389 return atom.getBonds().stream()
390 .map((Bond b) -> b.get1_2(atom))
391 .filter((Atom a) -> a != toExclude)
392 .filter((Atom a) -> a.getAtomType().atomicNumber == element)
393 .collect(Collectors.toList());
394 }
395
396
397
398
399
400
401
402 public static Atom findNitrogenAtom(Residue residue) {
403 assert residue.getResidueType() == Residue.ResidueType.AA;
404
405
406 List<Atom> nitrogenCandidates = new ArrayList<>(2);
407
408 switch (residue.getAminoAcid3()) {
409 case LYS, LYD -> {
410
411 List<Atom> nitrogenList = findAtomsOfElement(residue, 7);
412 for (Atom nitrogen : nitrogenList) {
413 List<Atom> carbons = findBondedAtoms(nitrogen, 6);
414 if (carbons.size() == 2) {
415 nitrogenCandidates.add(nitrogen);
416 } else if (findBondedAtoms(carbons.get(0), 1).size() < 2) {
417 nitrogenCandidates.add(nitrogen);
418 }
419 }
420 if (nitrogenCandidates.isEmpty()) {
421 throw new IllegalArgumentException(
422 format(" Could not identify N atom of residue %s!", residue));
423 }
424 }
425
426 case ARG, HIS, HIE, HID -> {
427
428
429
430
431 List<Atom> nitrogenList = findAtomsOfElement(residue, 7);
432 Atom commonC = findAtomsOfElement(residue, 6).stream()
433 .filter((Atom carbon) -> findBondedAtoms(carbon, 7).size() >= 2)
434 .findAny().get();
435 nitrogenCandidates = nitrogenList.stream()
436 .filter((Atom nitr) -> !atomAttachedToAtom(nitr, commonC))
437 .collect(Collectors.toList());
438 }
439 case ASN, GLN -> {
440
441
442
443
444 List<Atom> nitrogenList = findAtomsOfElement(residue, 7);
445 for (Atom nitrogen : nitrogenList) {
446 List<Atom> bondedCarbs = findBondedAtoms(nitrogen, 6);
447 for (Atom carbon : bondedCarbs) {
448 if (!hasAttachedAtom(carbon, 8)) {
449 nitrogenCandidates.add(nitrogen);
450 }
451 }
452 }
453 if (nitrogenCandidates.isEmpty()) {
454 throw new IllegalArgumentException(
455 format(" Could not identify N atom of residue %s!", residue));
456 }
457 }
458 case TRP -> {
459
460
461
462
463 List<Atom> nitrogenList = findAtomsOfElement(residue, 7);
464 for (Atom nitrogen : nitrogenList) {
465 List<Atom> bondedCarbs = findBondedAtoms(nitrogen, 6);
466 if (bondedCarbs.size() == 1) {
467 nitrogenCandidates.add(nitrogen);
468 }
469 for (Atom carbon : bondedCarbs) {
470 if (hasAttachedAtom(carbon, 8)) {
471 nitrogenCandidates.add(nitrogen);
472 }
473 }
474 }
475 if (nitrogenCandidates.isEmpty()) {
476 throw new IllegalArgumentException(
477 format(" Could not identify N atom of residue %s!", residue));
478 }
479 }
480 case ACE -> {
481 return null;
482 }
483
484 default -> nitrogenCandidates = findAtomsOfElement(residue, 7);
485 }
486
487 switch (nitrogenCandidates.size()) {
488 case 0 -> {
489 logger.warning(
490 " Did not find any atoms that might be the amide nitrogen for residue " + residue);
491 return null;
492 }
493 case 1 -> {
494 return nitrogenCandidates.get(0);
495 }
496 case 2 -> {
497 logger.fine(format(
498 " Probable NME C-terminal cap attached to residue %s, some atom names may be duplicated!",
499 residue));
500 Atom N = null;
501 for (Atom nitro : nitrogenCandidates) {
502 nitro.setName("N");
503 Optional<Atom> capMethyl = findBondedAtoms(nitro, 6).stream()
504 .filter((Atom carb) -> findBondedAtoms(carb, 1).size() == 3)
505 .findAny();
506 if (capMethyl.isPresent()) {
507 findBondedAtoms(nitro, 1).get(0).setName("H");
508 Atom theCap = capMethyl.get();
509 theCap.setName("CH3");
510 List<Atom> capHydrogenList = findBondedAtoms(theCap, 1);
511 for (int i = 0; i < 3; i++) {
512 capHydrogenList.get(i).setName(format("H%d", i + 1));
513 }
514 } else {
515 N = nitro;
516 }
517 }
518 return N;
519 }
520
521
522 default -> {
523 if(logger.isLoggable(Level.FINE)){
524 logger.warning(format(" Nitrogen could not be mapped for amide nitrogen of residue %s ", residue));
525 }
526 return null;
527 }
528 }
529 }
530
531
532
533
534
535
536
537
538 public static Optional<Atom> findNucleotideO4s(Residue residue) {
539 assert residue.getResidueType() == Residue.ResidueType.NA;
540 return findAtomsOfElement(residue, 8).stream()
541 .filter(o -> findBondedAtoms(o, 6).size() == 2)
542 .findAny();
543 }
544
545
546
547
548
549
550
551
552 public static Atom getAlphaCarbon(Residue residue, Atom N) {
553 List<Atom> resAtoms = residue.getAtomList();
554 List<Atom> caCandidates = findBondedAtoms(N, 6).stream().filter(resAtoms::contains).toList();
555
556 if (residue.getAminoAcid3() == AminoAcid3.PRO) {
557 Atom CA = null;
558 Atom CD = null;
559 Atom aceC = null;
560 for (Atom caCand : caCandidates) {
561 if (hasAttachedAtom(caCand, 8)) {
562 aceC = caCand;
563 } else {
564 List<Atom> attachedH = findBondedAtoms(caCand, 1);
565 if (attachedH.size() == 1) {
566 CA = caCand;
567 } else if (attachedH.size() == 2) {
568 CD = caCand;
569 } else {
570 throw new IllegalArgumentException(format(" Error in parsing proline %s", residue));
571 }
572 }
573 }
574 assert CA != null && CD != null;
575 if (aceC != null) {
576 nameAcetylCap(residue, aceC);
577 }
578 return CA;
579 }
580 if (caCandidates.size() == 1) {
581 return caCandidates.get(0);
582 } else {
583 Atom CA = null;
584 Atom aceC = null;
585 for (Atom caCand : caCandidates) {
586 if (hasAttachedAtom(caCand, 8)) {
587 aceC = caCand;
588 } else {
589 CA = caCand;
590 }
591 }
592 nameAcetylCap(residue, aceC);
593 return CA;
594 }
595 }
596
597
598
599
600
601
602
603
604 public static boolean hasAttachedAtom(Atom atom, int element) {
605 return atom.getBonds().stream()
606 .map((Bond b) -> b.get1_2(atom))
607 .anyMatch((Atom a) -> a.getAtomType().atomicNumber == element);
608 }
609
610
611
612
613
614
615
616
617
618
619
620
621
622 public static void intxyz(
623 Atom atom, Atom ia, double bond, Atom ib, double angle1, Atom ic, double angle2, int chiral) {
624 double[] xa = new double[3];
625 xa = (ia == null) ? null : ia.getXYZ(xa);
626 double[] xb = new double[3];
627 xb = (ib == null) ? null : ib.getXYZ(xb);
628 double[] xc = new double[3];
629 xc = (ic == null) ? null : ic.getXYZ(xc);
630 atom.moveTo(determineIntxyz(xa, bond, xb, angle1, xc, angle2, chiral));
631 }
632
633
634
635
636
637
638 public static void numberAtoms(MolecularAssembly molecularAssembly) {
639 int index = 1;
640 for (Atom a : molecularAssembly.getAtomArray()) {
641 a.setXyzIndex(index++);
642 }
643 index--;
644 if (logger.isLoggable(Level.INFO)) {
645 logger.info(format(" Total number of atoms: %d\n", index));
646 }
647
648 Polymer[] polymers = molecularAssembly.getChains();
649 if (polymers != null) {
650 for (Polymer p : polymers) {
651 List<Residue> residues = p.getResidues();
652 for (Residue r : residues) {
653 r.reOrderAtoms();
654 }
655 }
656 }
657 List<MSNode> molecules = molecularAssembly.getMolecules();
658 for (MSNode n : molecules) {
659 MSGroup m = (MSGroup) n;
660 m.reOrderAtoms();
661 }
662 List<MSNode> water = molecularAssembly.getWater();
663 for (MSNode n : water) {
664 MSGroup m = (MSGroup) n;
665 m.reOrderAtoms();
666 }
667 List<MSNode> ions = molecularAssembly.getIons();
668 for (MSNode n : ions) {
669 MSGroup m = (MSGroup) n;
670 m.reOrderAtoms();
671 }
672 }
673
674
675
676
677
678
679
680
681 public static Atom[] sortAtomsByDistance(Atom reference, List<Atom> toCompare) {
682 Atom[] theAtoms = toCompare.toArray(new Atom[0]);
683 sortAtomsByDistance(reference, theAtoms);
684 return theAtoms;
685 }
686
687
688
689
690
691
692
693 public static void sortAtomsByDistance(Atom reference, Atom[] toCompare) {
694 final double[] refXYZ = reference.getXYZ(new double[3]);
695 Arrays.sort(
696 toCompare,
697 Comparator.comparingDouble(
698 a -> {
699 double[] atXYZ = a.getXYZ(new double[3]);
700 return DoubleMath.dist2(refXYZ, atXYZ);
701 }));
702 }
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719 private static Atom buildHeavyAtom(MSGroup residue, String atomName, Atom ia, double bond,
720 Atom ib, double angle1, Atom ic, double angle2, int chiral, AtomType atomType) {
721 Atom atom = (Atom) residue.getAtomNode(atomName);
722 if (atomType == null) {
723 return null;
724 }
725 if (atom == null) {
726 String resName = ia.getResidueName();
727 int resSeq = ia.getResidueNumber();
728 Character chainID = ia.getChainID();
729 Character altLoc = ia.getAltLoc();
730 String segID = ia.getSegID();
731 double occupancy = ia.getOccupancy();
732 double tempFactor = ia.getTempFactor();
733 atom =
734 new Atom(
735 0,
736 atomName,
737 altLoc,
738 new double[3],
739 resName,
740 resSeq,
741 chainID,
742 occupancy,
743 tempFactor,
744 segID,
745 true);
746 residue.addMSNode(atom);
747 intxyz(atom, ia, bond, ib, angle1, ic, angle2, chiral);
748 }
749 atom.setAtomType(atomType);
750 return atom;
751 }
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770 private static Atom buildHeavyAtom(MSGroup residue, String atomName, Atom ia, double bond,
771 Atom ib, double angle1, Atom ic, double angle2, int chiral, AtomType atomType,
772 ForceField forceField, List<Bond> bondList) {
773 Atom atom =
774 buildHeavyAtom(residue, atomName, ia, bond, ib, angle1, ic, angle2, chiral, atomType);
775 buildBond(ia, atom, forceField, bondList);
776 return atom;
777 }
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800 public static double[] determineIntxyz(double[] ia, double bond, double[] ib, double angle1,
801 double[] ic, double angle2, int chiral) {
802 if (ia != null && !Double.isFinite(bond)) {
803 throw new IllegalArgumentException(
804 String.format(" Passed bond length is non-finite %f", bond));
805 } else if (ib != null && !Double.isFinite(angle1)) {
806 throw new IllegalArgumentException(String.format(" Passed angle is non-finite %f", angle1));
807 } else if (ic != null && !Double.isFinite(angle2)) {
808 throw new IllegalArgumentException(
809 String.format(" Passed dihedral/improper is non-finite %f", angle2));
810 }
811
812 if (chiral == 3) {
813 double[] negChiral = determineIntxyz(ia, bond, ib, angle1, ic, angle2, -1);
814 double[] posChiral = determineIntxyz(ia, bond, ib, angle1, ic, angle2, 1);
815 double[] displacement = new double[3];
816 double dispMag = 0;
817
818 for (int i = 0; i < 3; i++) {
819
820 displacement[i] = 0.5 * (posChiral[i] + negChiral[i]);
821
822 displacement[i] -= ia[i];
823
824 dispMag += (displacement[i] * displacement[i]);
825 }
826
827 dispMag = sqrt(dispMag);
828 double extend = bond / dispMag;
829 assert extend > 0.999;
830
831 double[] outXYZ = new double[3];
832 for (int i = 0; i < 3; i++) {
833 displacement[i] *= extend;
834 outXYZ[i] = displacement[i] + ia[i];
835 }
836 return outXYZ;
837 }
838
839 angle1 = toRadians(angle1);
840 angle2 = toRadians(angle2);
841 double zcos0 = cos(angle1);
842 double zcos1 = cos(angle2);
843 double zsin0 = sin(angle1);
844 double zsin1 = sin(angle2);
845 double[] ret = new double[3];
846 double[] x = new double[3];
847
848
849 if (ia == null) {
850 x[0] = x[1] = x[2] = 0.0;
851 } else if (ib == null) {
852 double[] xa = new double[3];
853 arraycopy(ia, 0, xa, 0, ia.length);
854
855 x[0] = xa[0];
856 x[1] = xa[1];
857 x[2] = xa[2] + bond;
858 } else if (ic == null) {
859 double[] xa = new double[3];
860 double[] xb = new double[3];
861 double[] xab = new double[3];
862 for (int i = 0; i < ia.length; i++) {
863 xa[i] = ia[i];
864 xb[i] = ib[i];
865 }
866
867 sub(xa, xb, xab);
868 double rab = length(xab);
869 normalize(xab, xab);
870 double cosb = xab[2];
871 double sinb = sqrt(xab[0] * xab[0] + xab[1] * xab[1]);
872 double cosg, sing;
873 if (sinb == 0.0d) {
874 cosg = 1.0d;
875 sing = 0.0d;
876 } else {
877 cosg = xab[1] / sinb;
878 sing = xab[0] / sinb;
879 }
880 double xtmp = bond * zsin0;
881 double ztmp = rab - bond * zcos0;
882 x[0] = xb[0] + xtmp * cosg + ztmp * sing * sinb;
883 x[1] = xb[1] - xtmp * sing + ztmp * cosg * sinb;
884 x[2] = xb[2] + ztmp * cosb;
885 } else if (chiral == 0) {
886 double[] xa = new double[3];
887 double[] xb = new double[3];
888 double[] xc = new double[3];
889 double[] xab = new double[3];
890 double[] xbc = new double[3];
891 double[] xt = new double[3];
892 double[] xu = new double[3];
893 for (int i = 0; i < ia.length; i++) {
894 xa[i] = ia[i];
895 xb[i] = ib[i];
896 xc[i] = ic[i];
897 }
898
899 sub(xa, xb, xab);
900 normalize(xab, xab);
901 sub(xb, xc, xbc);
902 normalize(xbc, xbc);
903 xt[0] = xab[2] * xbc[1] - xab[1] * xbc[2];
904 xt[1] = xab[0] * xbc[2] - xab[2] * xbc[0];
905 xt[2] = xab[1] * xbc[0] - xab[0] * xbc[1];
906 double cosine = xab[0] * xbc[0] + xab[1] * xbc[1] + xab[2] * xbc[2];
907 double sine = sqrt(max(1.0d - cosine * cosine, eps));
908 if (abs(cosine) >= 1.0d) {
909 logger.warning("Undefined Dihedral");
910 }
911 scale(xt, 1.0d / sine, xt);
912 xu[0] = xt[1] * xab[2] - xt[2] * xab[1];
913 xu[1] = xt[2] * xab[0] - xt[0] * xab[2];
914 xu[2] = xt[0] * xab[1] - xt[1] * xab[0];
915 x[0] = xa[0] + bond * (xu[0] * zsin0 * zcos1 + xt[0] * zsin0 * zsin1 - xab[0] * zcos0);
916 x[1] = xa[1] + bond * (xu[1] * zsin0 * zcos1 + xt[1] * zsin0 * zsin1 - xab[1] * zcos0);
917 x[2] = xa[2] + bond * (xu[2] * zsin0 * zcos1 + xt[2] * zsin0 * zsin1 - xab[2] * zcos0);
918 } else if (abs(chiral) == 1) {
919 double[] xa = new double[3];
920 double[] xb = new double[3];
921 double[] xc = new double[3];
922 double[] xba = new double[3];
923 double[] xac = new double[3];
924 double[] xt = new double[3];
925 for (int i = 0; i < ia.length; i++) {
926 xa[i] = ia[i];
927 xb[i] = ib[i];
928 xc[i] = ic[i];
929 }
930 sub(xb, xa, xba);
931 normalize(xba, xba);
932 sub(xa, xc, xac);
933 normalize(xac, xac);
934 xt[0] = xba[2] * xac[1] - xba[1] * xac[2];
935 xt[1] = xba[0] * xac[2] - xba[2] * xac[0];
936 xt[2] = xba[1] * xac[0] - xba[0] * xac[1];
937 double cosine = xba[0] * xac[0] + xba[1] * xac[1] + xba[2] * xac[2];
938 double sine2 = max(1.0d - cosine * cosine, eps);
939 if (abs(cosine) >= 1.0d) {
940 logger.warning("Defining Atom Colinear");
941 }
942 double a = (-zcos1 - cosine * zcos0) / sine2;
943 double b = (zcos0 + cosine * zcos1) / sine2;
944 double c = (1.0d + a * zcos1 - b * zcos0) / sine2;
945 if (c > eps) {
946 c = chiral * sqrt(c);
947 } else if (c < -eps) {
948 c =
949 sqrt(
950 (a * xac[0] + b * xba[0]) * (a * xac[0] + b * xba[0])
951 + (a * xac[1] + b * xba[1]) * (a * xac[1] + b * xba[1])
952 + (a * xac[2] + b * xba[2]) * (a * xac[2] + b * xba[2]));
953 a /= c;
954 b /= c;
955 c = 0.0d;
956 } else {
957 c = 0.0d;
958 }
959 x[0] = xa[0] + bond * (a * xac[0] + b * xba[0] + c * xt[0]);
960 x[1] = xa[1] + bond * (a * xac[1] + b * xba[1] + c * xt[1]);
961 x[2] = xa[2] + bond * (a * xac[2] + b * xba[2] + c * xt[2]);
962 }
963 arraycopy(x, 0, ret, 0, ret.length);
964 return ret;
965 }
966
967
968 public static class MissingAtomTypeException extends Exception {
969
970 @Serial
971 private static final long serialVersionUID = 1L;
972
973 public final Residue residue;
974 public final Atom atom;
975
976 public MissingAtomTypeException(Residue residue, Atom atom) {
977 this.residue = residue;
978 this.atom = atom;
979 }
980
981 @Override
982 public String toString() {
983 StringBuilder sb = new StringBuilder();
984 sb.append(format(" Atom %s", atom));
985 if (residue != null) {
986 sb.append(format("\n of residue %c-%s", residue.getChainID(), residue));
987 }
988 sb.append("\n could not be assigned an atom type.\n");
989 return sb.toString();
990 }
991 }
992
993
994 public static class MissingHeavyAtomException extends Exception {
995
996 @Serial
997 private static final long serialVersionUID = 1L;
998
999 public final String atomName;
1000 public final AtomType atomType;
1001 final Atom bondedTo;
1002
1003 public MissingHeavyAtomException(String atomName, AtomType atomType, Atom bondedTo) {
1004 this.atomName = atomName;
1005 this.atomType = atomType;
1006 this.bondedTo = bondedTo;
1007 }
1008
1009 @Override
1010 public String toString() {
1011 StringBuilder sb = new StringBuilder();
1012 if (atomType != null) {
1013 sb.append(format("\n An atom of type\n %s\n", atomType));
1014 } else {
1015 sb.append(format("\n Atom %s", atomName));
1016 }
1017 sb.append(" was not found");
1018 if (bondedTo != null) {
1019 sb.append(format(" bonded to atom %s ", bondedTo));
1020 }
1021 sb.append(".\n");
1022 return sb.toString();
1023 }
1024 }
1025 }