-
Notifications
You must be signed in to change notification settings - Fork 0
/
Molecule.java
1192 lines (1058 loc) · 47.7 KB
/
Molecule.java
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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import java.util.*;
import java.io.*;
import org.apache.commons.math3.geometry.euclidean.threed.*;
import org.jgrapht.*;
import org.jgrapht.graph.*;
import org.jgrapht.alg.*;
import com.google.common.collect.*;
import Jama.*;
/**
* Represents a molecule. This class is effectively immutable.
*/
public class Molecule implements Immutable, Serializable
{
/** For serialization. */
public static final long serialVersionUID = 1L;
/** Name of the molecule. */
public final String name;
/** The atoms in the molecule. */
public final List<Atom> contents;
/** The connectivity graph. */
public final SimpleWeightedGraph<Atom, DefaultWeightedEdge> connectivity;
/**
* Factory method to create a molecule given a map of old atoms to new atoms. Should be used
* to move atoms.
* @param atomMap a map from old atoms to new atoms (does not have to include all atoms)
*/
public Molecule moveAtoms(Map<Atom,Atom> atomMap)
{
// copy the list of vertices
List<Atom> newContents = new ArrayList<Atom>(contents.size());
for (Atom a : contents)
{
if ( atomMap.containsKey(a) )
newContents.add(atomMap.get(a));
else
newContents.add(a);
}
// populate a new connectivity graph
SimpleWeightedGraph<Atom,DefaultWeightedEdge> newConnectivity = new SimpleWeightedGraph<Atom,DefaultWeightedEdge>(DefaultWeightedEdge.class);
for (Atom newAtom : newContents)
newConnectivity.addVertex(newAtom);
for (DefaultWeightedEdge e : connectivity.edgeSet())
{
// get old edge data
Double bondOrder = connectivity.getEdgeWeight(e);
Atom fromAtom = connectivity.getEdgeSource(e);
Atom toAtom = connectivity.getEdgeTarget(e);
// replace any changes
if ( atomMap.containsKey(fromAtom) )
fromAtom = atomMap.get(fromAtom);
if ( atomMap.containsKey(toAtom) )
toAtom = atomMap.get(toAtom);
// create new edge
if (! newContents.contains(fromAtom) || ! newContents.contains(toAtom))
throw new IllegalArgumentException("edge not in graph FromAtom: " + getAtomString(fromAtom) + " ToAtom: " + getAtomString(toAtom));
DefaultWeightedEdge newEdge = newConnectivity.addEdge(fromAtom,toAtom);
newConnectivity.setEdgeWeight(newEdge, bondOrder);
}
// return result
return new Molecule(name, newContents, newConnectivity);
}
/** Adds the specified shift to all the atoms. */
public Molecule shift(Vector3D shift)
{
Map<Atom,Atom> atomMap = new HashMap<>();
for (Atom a : contents)
atomMap.put(a, a.moveAtom(a.position.add(shift)));
return moveAtoms(atomMap);
}
/**
* A shallow clone of a connectivity graph.
* @param oldConnectivity the old connectivity graph
* @return the clone
*/
public static SimpleWeightedGraph<Atom,DefaultWeightedEdge> cloneConnectivity(SimpleWeightedGraph<Atom,DefaultWeightedEdge> oldConnectivity)
{
// populate a new connectivity graph
SimpleWeightedGraph<Atom,DefaultWeightedEdge> newConnectivity = new SimpleWeightedGraph<Atom,DefaultWeightedEdge>(DefaultWeightedEdge.class);
for (Atom oldAtom : oldConnectivity.vertexSet())
newConnectivity.addVertex(oldAtom);
for (DefaultWeightedEdge e : oldConnectivity.edgeSet())
{
Double bondOrder = oldConnectivity.getEdgeWeight(e);
Atom fromAtom = oldConnectivity.getEdgeSource(e);
Atom toAtom = oldConnectivity.getEdgeTarget(e);
DefaultWeightedEdge newEdge = newConnectivity.addEdge(fromAtom,toAtom);
newConnectivity.setEdgeWeight(newEdge, bondOrder);
}
return newConnectivity;
}
/**
* Factory method to create a new molecule by shifting this one to have its center at the origin.
* @return this minus the barycenter of this
*/
public Molecule normalize()
{
Vector3D barycenter = Vector3D.ZERO;
for (Atom a : contents) barycenter = barycenter.add(a.position);
barycenter = barycenter.scalarMultiply(1.0/contents.size());
return shift(barycenter.negate());
}
/**
* Returns the centroid of this molecule.
* @return the centroid vector
*/
public Vector3D getCentroid()
{
Vector3D barycenter = Vector3D.ZERO;
for (Atom a : contents) barycenter = barycenter.add(a.position);
barycenter = barycenter.scalarMultiply(1.0/contents.size());
return barycenter;
}
/**
* Constructs a new molecule. The connectivity graph will be used directly (i.e., no defensive copy is made).
* @param name the name of the molecule
* @param contents the atoms in the molecule
* @param connectivity the connectivity graph
*/
public Molecule(String name, List<Atom> contents, SimpleWeightedGraph<Atom,DefaultWeightedEdge> connectivity)
{
this.name = name;
this.contents = ImmutableList.copyOf(contents);
this.connectivity = connectivity;
}
/**
* Returns the number of the atom. Answers are given in 1,2,...n where n is
* the number of atoms in the molecule.
* @param atom an Atom in this Molecule
* @return the requested atom number (0 if it is not in the molecule)
*/
public int getAtomNumber(Atom atom)
{
return contents.indexOf(atom) + 1;
}
/**
* Returns the element and atom number of an atom. e.g., C5.
* Throws an IllegalArgumentException if the Atom is not in the molecule.
* @param atom the atom whose string is desired
* @return the description of the atom
*/
public String getAtomString(Atom atom)
{
if ( atom == null)
throw new NullPointerException("atom cannot be null");
int atomNumber = getAtomNumber(atom);
if ( atomNumber == 0 )
throw new IllegalArgumentException("cannot print out atom because it is not in the molecule!\n" + atom.toString());
return atom.element.symbol + atomNumber;
}
/** For debugging. */
public String getClosestAtomString(Atom atom)
{
double closestDistance = 100.0;
Atom closestAtom = null;
for (Atom a : contents)
{
if ( a.type1 != atom.type1 )
continue;
double distance = Vector3D.distance(atom.position, a.position);
if ( distance < closestDistance )
{
closestDistance = distance;
closestAtom = a;
}
}
if ( closestDistance > 0.1 || closestAtom == null )
return atom.toString();
else
return getAtomString(closestAtom);
}
/**
* Given a single atom, returns a set of all connected atoms, including the
* specified atom. Uses a breadth-first search algorithm.
* @param startingAtom the atom to start exploring the connectivity from
* @return the atoms in the subgraph of the startingAtom
*/
public Set<Atom> exploreGraph(Atom startingAtom)
{
if ( ! contents.contains(startingAtom) )
throw new IllegalArgumentException("cannot search the connectivity graph because the specified atom is not in this molecule");
Set<Atom> returnSet = new HashSet<>();
LinkedList<Atom> searchQueue = new LinkedList<>();
searchQueue.add(startingAtom);
// breadth-first search
while (searchQueue.size() > 0)
{
Atom currentNode = searchQueue.remove();
for (Atom a : getAdjacentAtoms(currentNode))
{
if ( ! returnSet.contains(a) )
{
returnSet.add(a);
searchQueue.add(a);
}
}
}
return returnSet;
}
/**
* Given a bond between includeAtom and excludeAtom, returns a set of atoms
* containing all the atoms on the includeAtom side of the bond, including
* includeAtom. Will throw an exception if includeAtom and excludeAtom are not
* directly bonded. Will throw an exception if includeAtom and excludeAtom
* form a ring.
* @param excludeAtom this atom will not be included in the result
* @param includeAtom this atom will be included in the result
* @return the atoms on the includeAtom side of the graph
*/
public Set<Atom> getHalfGraph(Atom excludeAtom, Atom includeAtom)
{
// if these atoms are not directly bonded, then return an empty set
if ( !directlyConnected(excludeAtom, includeAtom) )
throw new IllegalArgumentException(String.format("atoms %s and %s not directly connected", getAtomString(excludeAtom), getAtomString(includeAtom)));
// preform a breadth-first search of one branch of the graph only
Set<Atom> returnSet = new HashSet<Atom>();
LinkedList<Atom> searchQueue = new LinkedList<Atom>();
Set<Atom> searched = new HashSet<>();
for (Atom a : getAdjacentAtoms(includeAtom))
{
searchQueue.add(a);
returnSet.add(a);
}
searchQueue.remove(excludeAtom);
returnSet.remove(excludeAtom);
returnSet.add(includeAtom);
searched.add(includeAtom);
Atom lastNode = includeAtom;
while (searchQueue.size() > 0)
{
Atom currentNode = searchQueue.remove();
Set<Atom> adjacent = getAdjacentAtoms(currentNode);
adjacent.remove(lastNode);
for (Atom a : adjacent)
{
if ( a == excludeAtom)
{
// if the excluded atom is found, this is a ring!
GaussianInputFile gjf = new GaussianInputFile(this);
gjf.write("error.gjf");
throw new IllegalArgumentException("includeAtom " + getAtomString(includeAtom) +
" and excludeAtom " + getAtomString(excludeAtom) + " cannot form a ring!");
}
// if this is an atom we haven't already searched, mark it
if ( ! searched.contains(a) )
{
searchQueue.add(a);
searched.add(a);
returnSet.add(a);
}
}
}
return(returnSet);
}
/**
* {@link #getHalfGraph(Atom,Atom)}
* @param index1 the index of the excludeAtom
* @param index2 the index of the includeAtom
* @return the atoms on the includeAtom side of the graph
*/
public Set<Atom> getHalfGraph(int index1, int index2)
{
Atom atom1 = contents.get(index1);
Atom atom2 = contents.get(index2);
if ( atom1 == null || atom2 == null )
throw new NullPointerException("atoms not found in this molecule for half graph call");
return getHalfGraph(atom1, atom2);
}
/**
* Returns atom indices instead of atoms.
* {@link #getHalfGraph(Atom,Atom)}
* @param excludeAtom this atom will not be included in the result
* @param includeAtom this atom will be included in the result
* @return the atom indices on the includeAtom side of the graph
*/
public Set<Integer> getHalfGraphIndices(Atom excludeAtom, Atom includeAtom)
{
Set<Atom> set = getHalfGraph(excludeAtom, includeAtom);
Set<Integer> returnSet = new HashSet<>();
for (Atom a : set)
returnSet.add(contents.indexOf(a));
return returnSet;
}
/**
* Returns atom indices instead of atoms.
* {@link #getHalfGraph(Atom,Atom)}
* @param index1 the excludeAtom index
* @param index2 the includeAtom index
* @return the atom indices on the includeAtom side of the graph
*/
public Set<Integer> getHalfGraphIndices(int index1, int index2)
{
return getHalfGraphIndices(contents.get(index1), contents.get(index2));
}
/**
* Determines whether atom1 and atom2 share an edge (i.e., are bonded).
* No exception is thrown if these atoms aren't in the graph.
* @param atom1 test whether this atom is connected to the other atom
* @param atom2 test whether this atom is connected to the other atom
* @return true if the atoms are bonded
*/
public boolean directlyConnected(Atom atom1, Atom atom2)
{
DefaultWeightedEdge e = connectivity.getEdge(atom1,atom2);
if ( e == null )
return false;
return true;
}
/**
* {@link #directlyConnected(Atom,Atom)}
* @param i index of atom1
* @param j index of atom2
* @return true if atoms are bonded
*/
public boolean directlyConnected(int i, int j)
{
return directlyConnected(contents.get(i), contents.get(j));
}
/**
* Returns the distance between two atoms.
* @param atom1 one of the two atoms
* @param atom2 one of the two atoms
* @return the distance between the atoms in angstroms
*/
public static double getDistance(Atom atom1, Atom atom2)
{
return Vector3D.distance(atom1.position, atom2.position);
}
/**
* Returns the angle between three atoms.
* @param atom1 one of the three atoms
* @param atom2 one of the three atoms
* @param atom3 one of the three atoms
* @return the angle between the atoms in degrees
*/
public static double getAngle(Atom atom1, Atom atom2, Atom atom3)
{
return getAngle(atom1.position, atom2.position, atom3.position);
}
/**
* Returns the v1-v2-v3 angle.
* @param v1 the first point
* @param v2 the second point
* @param v3 the third point
* @return the angle in degrees
*/
public static double getAngle(Vector3D v1, Vector3D v2, Vector3D v3)
{
Vector3D v1prime = v1.subtract(v2);
Vector3D v3prime = v3.subtract(v2);
return Math.toDegrees(Vector3D.angle(v1prime, v3prime));
}
/**
* Returns the angle between three atoms.
* @param index1 the index of atom1
* @param index2 the index of atom2
* @param index3 the index of atom3
* @return the angle in degrees
*/
public double getAngle(int index1, int index2, int index3)
{
return getAngle(contents.get(index1), contents.get(index2), contents.get(index3));
}
/**
* Moves the group associated with atom2 to the specified distance.
* Motion occurs along the atom1-atom2 bond vector. Note that this returns a new
* molecule. No checks are made.
* @param atom1 this atom will be held fixed
* @param atom2 this atom and anything connected to it will be moved
* @param requestedDistance the requested distance in Angstroms
* @return a new Molecule containing the same connectivity but new positions
*/
public Molecule setDistance(Atom atom1, Atom atom2, double requestedDistance)
{
// determine which atoms have to be moved
Set<Atom> toBeMoved = getHalfGraph(atom1, atom2);
// determine how much to move the atoms
Vector3D oldPosition1 = atom1.position;
Vector3D oldPosition2 = atom2.position;
double currentDistance = getDistance(atom1, atom2);
Vector3D translateVector = oldPosition2.subtract(oldPosition1);
double scaling = (requestedDistance - currentDistance)/currentDistance;
Vector3D requiredTranslation = translateVector.scalarMultiply(scaling);
Map<Atom,Atom> atomMap = new HashMap<>();
for (Atom oldAtom : toBeMoved)
{
Vector3D oldPosition = oldAtom.position;
Vector3D newPosition = oldPosition.add(requiredTranslation);
Atom newAtom = oldAtom.moveAtom(newPosition);
atomMap.put(oldAtom, newAtom);
}
return moveAtoms(atomMap);
}
/**
* Moves the group associated with atom2 to the specified distance.
* Motion occurs along the atom1-atom2 bond vector. Note that this returns a new
* molecule. No checks are made. Indices are 0, 1, ..., n-1 (the as the contents field).
* @param index1 this atom index will be held fixed
* @param index2 this atom indexand anything connected to it will be moved
* @param requestedDistance the requested distance in Angstroms
* @return a new Molecule containing the same connectivity but new positions
*/
public Molecule setDistance(int index1, int index2, double requestedDistance)
{
return setDistance(contents.get(index1), contents.get(index2), requestedDistance);
}
/**
* Rotates the atom1-atom2-atom3 angle, moving only atom3 and anything in its
* attached subgraph. No checks.
* @param atom1 will not be moved
* @param atom2 will not be moved
* @param atom3 will be moved
* @param theta rotation in degrees
* @return a new Molecule containing the same connectivity but new positions
*/
public Molecule rotateAngle(Atom atom1, Atom atom2, Atom atom3, double theta)
{
// figure out which atoms to move
Set<Atom> toBeMoved = getHalfGraph(atom2, atom3);
// create atom map
LinkedHashMap<Atom,Atom> atomMap = new LinkedHashMap<Atom,Atom>();
// move everything to put atom2 at the origin
Vector3D v1 = null;
Vector3D v3 = null;
for (Atom a : contents)
{
Vector3D oldPosition = a.position;
Vector3D newPosition = oldPosition.subtract(atom2.position);
if ( toBeMoved.contains(a) )
atomMap.put(a,a.moveAtom(newPosition));
if ( a == atom1 )
v1 = newPosition;
else if ( a == atom3 )
v3 = newPosition;
}
// form the rotation axis and matrix
Vector3D rotationAxis = Vector3D.crossProduct(v1, v3);
Rotation rotation = new Rotation(rotationAxis, Math.toRadians(theta));
// apply rotation and undo translation
LinkedHashMap<Atom,Atom> atomMap2 = new LinkedHashMap<Atom,Atom>();
for (Atom a : atomMap.keySet())
{
Vector3D oldPosition = atomMap.get(a).position;
Vector3D newPosition = rotation.applyTo(oldPosition);
newPosition = newPosition.add(atom2.position);
atomMap2.put(a,a.moveAtom(newPosition));
}
// create new Molecule
return moveAtoms(atomMap2);
}
/**
* Method alias. Indices are 0, 1, ..., n-1. No checks.
*/
public Molecule rotateAngle(int i, int j, int k, double theta)
{
return rotateAngle(contents.get(i), contents.get(j), contents.get(k), theta);
}
/**
* Set the atom1-atom2-atom3 angle to theta degrees, moving atom3 and its subgraph only.
* New molecule returned. No checks.
* @param atom1 not moved
* @param atom2 not moved
* @param atom3 moved
* @param theta desired angle in degrees
*/
public Molecule setAngle(Atom atom1, Atom atom2, Atom atom3, double theta)
{
double currentAngle = getAngle(atom1, atom2, atom3);
double requiredRotation = theta - currentAngle;
return rotateAngle(atom1, atom2, atom3, requiredRotation);
}
/**
* Method alias. Indices are 1,2,...,n. No checks.
*/
public Molecule setAngle(int i, int j, int k, double theta)
{
return setAngle(contents.get(i), contents.get(j), contents.get(k), theta);
}
/**
* Returns a new Molecule with a rotated dihedral.
* Note that the old ProtoTorsion will no longer point to the new Molecule.
* @param protoTorsion the torsion to rotate
* @param theta the desired dihedral angle in degrees
* @return the new molecule
*/
public Molecule setDihedral(ProtoTorsion protoTorsion, double theta)
{
// get fields
Atom atom1 = protoTorsion.atom1;
Atom atom2 = protoTorsion.atom2;
Atom atom3 = protoTorsion.atom3;
Atom atom4 = protoTorsion.atom4;
Set<Atom> atomsToRotate = getHalfGraph(atom2,atom3);
// determine how much rotation is needed
double currentDihedralAngle = protoTorsion.getDihedralAngle();
double requiredRotation = currentDihedralAngle - theta;
if ( Math.abs(requiredRotation) < 0.001 )
return this;
Map<Atom,Atom> atomMap = new HashMap<>();
// move atom 3 to the origin
// define the rotation axis as the vector from atom3 (now at origin) to atom2
Vector3D rotationAxis = null;
for (Atom a : contents)
{
Vector3D oldPosition = a.position;
Vector3D newPosition = oldPosition.subtract(atom3.position);
if ( atomsToRotate.contains(a) )
atomMap.put(a,a.moveAtom(newPosition));
if ( a == atom2 )
rotationAxis = newPosition;
}
// rotate the atoms and make a new atom map
Map<Atom,Atom> atomMap2 = new HashMap<>();
Rotation rotation = new Rotation(rotationAxis, Math.toRadians(requiredRotation));
for (Atom a : atomMap.keySet())
{
// update rotation
Vector3D oldPosition = atomMap.get(a).position;
Vector3D newPosition = rotation.applyTo(oldPosition);
// undo translation
newPosition = newPosition.add(atom3.position);
// update map
atomMap2.put(a, a.moveAtom(newPosition));
}
// return new Molecule
return moveAtoms(atomMap2);
}
/**
* Returns a new Molecule with a rotated dihedral. Alias method.
* Note that the old IndexTorsion will still be valid for the new Molecule.
* @param indexTorsion the torsion to rotate
* @param theta the angle to set the torsion to in degrees
*/
public Molecule setDihedral(IndexTorsion indexTorsion, double theta)
{
// figure out which atoms to move
Set<Atom> atomsToRotate = new HashSet<>();
for (Integer i : indexTorsion.atomIndicesToRotate)
atomsToRotate.add(contents.get(i));
// get prototorsion
ProtoTorsion protoTorsion = indexTorsion.getProtoTorsion(this);
Atom atom1 = protoTorsion.atom1;
Atom atom2 = protoTorsion.atom2;
Atom atom3 = protoTorsion.atom3;
Atom atom4 = protoTorsion.atom4;
// determine how much rotation is needed
double currentDihedralAngle = protoTorsion.getDihedralAngle();
double requiredRotation = currentDihedralAngle - theta;
if ( Math.abs(requiredRotation) < 0.001 )
return this;
Map<Atom,Atom> atomMap = new HashMap<>();
// move atom 3 to the origin
// define the rotation axis as the vector from atom3 (now at origin) to atom2
Vector3D rotationAxis = null;
for (Atom a : contents)
{
Vector3D oldPosition = a.position;
Vector3D newPosition = oldPosition.subtract(atom3.position);
if ( atomsToRotate.contains(a) )
atomMap.put(a,a.moveAtom(newPosition));
if ( a == atom2 )
rotationAxis = newPosition;
}
// rotate the atoms and make a new atom map
Map<Atom,Atom> atomMap2 = new HashMap<>();
Rotation rotation = new Rotation(rotationAxis, Math.toRadians(requiredRotation));
for (Atom a : atomMap.keySet())
{
// update rotation
Vector3D oldPosition = atomMap.get(a).position;
Vector3D newPosition = rotation.applyTo(oldPosition);
// undo translation
newPosition = newPosition.add(atom3.position);
// update map
atomMap2.put(a, a.moveAtom(newPosition));
}
// return new Molecule
return moveAtoms(atomMap2);
}
/**
* Creates a new Molecule where atom2 and its subgraph have been moved to
* make atom1 sp2-hybridized (bond angles set at 120 degrees).
* Note that the new center will be sp2, but could have distorted torsion angles.
* @param atom1 the atom to be adjusted to sp2
* @param atom2 the group to be moved
* @param forceAngle true if we want to force the atom1alpha-atom1-atom1beta angle to 120 (safe for non-prolines)
* @return a new Molecule with adjusted hybridization
*/
public Molecule set_sp2(Atom atom1, Atom atom2, boolean forceAngle)
{
// note current bond length
double currentLength = getDistance(atom1, atom2);
// get the neighbors of atom1
// call them atom1alpha and atom1beta
List<Atom> atom1neighbors = new LinkedList<>(getAdjacentAtoms(atom1));
if ( atom1neighbors.size() != 3 )
throw new IllegalArgumentException("expected 3 neighbors for atom 1, found " + atom1neighbors.size());
else if ( ! atom1neighbors.contains(atom2) )
throw new IllegalArgumentException("atoms are not adjacent");
atom1neighbors.remove(atom2);
Atom atom1alpha = atom1neighbors.get(0);
Atom atom1beta = atom1neighbors.get(1);
int atom1alphaIndex = contents.indexOf(atom1alpha);
int atom1betaIndex = contents.indexOf(atom1beta);
int atom1index = contents.indexOf(atom1);
int atom2index = contents.indexOf(atom2);
// force the existing bond angle to 120 degrees
Molecule newMolecule = this;
if (forceAngle)
newMolecule = setAngle(atom1alpha, atom1, atom1beta, 120.0);
//System.out.println(getAtomString(atom1alpha));
//System.out.println(getAtomString(atom1));
//System.out.println(getAtomString(atom1beta));
// translate atom1 to the origin
List<Vector3D> newPositions = new LinkedList<>();
for (Atom a : newMolecule.contents)
{
Vector3D oldPosition = a.position;
Vector3D newPosition = oldPosition.subtract(atom1.position);
newPositions.add(newPosition);
}
// get unit vectors for atom1alpha and atom1beta
Vector3D atom1alphaPosition = newPositions.get(atom1alphaIndex).normalize();
Vector3D atom1betaPosition = newPositions.get(atom1betaIndex).normalize();
// get the cross product of atom1alpha and atom1beta
Vector3D atom1crossPosition = Vector3D.crossProduct(atom1alphaPosition, atom1betaPosition).normalize();
// find the linear transformation matrix that rotates atom1alphaPosition to a=(-sqrt(3)/2, 0.5, 0.0)
// and atom1betaPosition to b=(-0.5, sqrt(3)/2, 0.0). If these are two points of an equilateral triangle,
// the third is at c=(1,0,0).
// solve the simultaneous matrix equations:
// T atom1alphaPosition = a
// T atom1betaPosition = b
// T atom1cross = c
//
// call atom1alphaPosition A, atom1betaPosition B, and atom1cross C.
// this is equivalent to solving:
//
// T [ ABC ] = [ abc ]
//
// where this means concatenated column vectors. Therefore,
// T = [ abc ] [ ABC ]^-1
double[][] preMatrix_abc = { { -0.5, Math.sqrt(3.0)/2.0, 0.0 }, {-0.5, -1.0 * Math.sqrt(3.0)/2.0, 0.0}, {0.0, 0.0, 1.0} };
Matrix matrix_abc = new Matrix(preMatrix_abc);
matrix_abc = matrix_abc.transpose();
double[][] preMatrix_ABC = { { atom1alphaPosition.getX(), atom1alphaPosition.getY(), atom1alphaPosition.getZ() },
{ atom1betaPosition.getX(), atom1betaPosition.getY(), atom1betaPosition.getZ() },
{ atom1crossPosition.getX(), atom1crossPosition.getY(), atom1crossPosition.getZ() } };
Matrix matrix_ABC = new Matrix(preMatrix_ABC);
matrix_ABC = matrix_ABC.transpose();
Matrix matrix_ABC_inverse = matrix_ABC.inverse();
Matrix T = matrix_abc.times(matrix_ABC_inverse);
Matrix Tinverse = T.inverse();
// apply the inverse of T to (1,0,0) to get the third vertex of the triangle
double[][] preMatrix_c = { { 1.0, 0.0, 0.0 } };
Matrix matrix_c = new Matrix(preMatrix_c);
matrix_c = matrix_c.transpose();
Matrix thirdVertex = Tinverse.times(matrix_c);
Vector3D thirdVertexPosition = new Vector3D( thirdVertex.get(0, 0), thirdVertex.get(1, 0), thirdVertex.get(2, 0) );
//double angle1 = getAngle(newPositions.get(atom1alphaNumber-1), newPositions.get(atom1number-1), newPositions.get(atom1betaNumber-1));
//double angle2 = getAngle(newPositions.get(atom1betaNumber-1), newPositions.get(atom1number-1), newPositions.get(atom2number-1));
//double angle3 = getAngle(newPositions.get(atom2number-1), newPositions.get(atom1number-1), newPositions.get(atom1alphaNumber-1));
//System.out.println(angle1);
//System.out.println(angle2);
//System.out.println(angle3);
// calculate the necessary rotation
Vector3D atom2position = newPositions.get(atom2index);
Vector3D rotationAxis = Vector3D.crossProduct( thirdVertexPosition, atom2position );
double requiredTheta = Vector3D.angle( thirdVertexPosition, atom2position );
Rotation fixRotation = new Rotation( rotationAxis, -1.0 * requiredTheta );
// determine which atoms should be moved
Set<Integer> atomIndicesToMove = getHalfGraphIndices(atom1, atom2);
List<Vector3D> newPositions2 = new LinkedList<>();
for (int i=0; i < newPositions.size(); i++)
{
if ( atomIndicesToMove.contains(i) )
{
// rotate this atom
Vector3D oldPosition = newPositions.get(i);
Vector3D newPosition = fixRotation.applyTo(oldPosition);
newPositions2.add(newPosition);
}
else
{
// do not rotate this atom
newPositions2.add( newPositions.get(i) );
}
}
// undo translation
List<Vector3D> newPositions3 = new LinkedList<>();
for (Vector3D v : newPositions2)
newPositions3.add(v.add(atom1.position));
/*double angle1 = getAngle(newPositions3.get(atom1alphaIndex), newPositions3.get(atom1index), newPositions3.get(atom1betaIndex));
double angle2 = getAngle(newPositions3.get(atom1betaIndex), newPositions3.get(atom1index), newPositions3.get(atom2index));
double angle3 = getAngle(newPositions3.get(atom2index), newPositions3.get(atom1index), newPositions3.get(atom1alphaIndex));
System.out.println(angle1);
System.out.println(angle2);
System.out.println(angle3);*/
//System.out.println(atom1alphaNumber + " " + atom1number + " " + atom1betaNumber);
//System.out.println(atom1betaNumber + " " + atom1number + " " + atom2number);
//System.out.println(atom2number + " " + atom1number + " " + atom1alphaNumber);
// create new atom map
Map<Atom,Atom> newAtomMap = new HashMap<>();
for (int i=0; i < contents.size(); i++)
{
Atom oldAtom = contents.get(i);
Atom newAtom = oldAtom.moveAtom( newPositions3.get(i) );
if ( !oldAtom.equals(newAtom) )
newAtomMap.put(oldAtom, newAtom);
}
Molecule rotatedMolecule = moveAtoms(newAtomMap);
//for (Atom key : newAtomMap.keySet())
// System.out.println(String.format("%s %s : %s %s", getAtomString(key), key, rotatedMolecule.getAtomString(newAtomMap.get(key)), newAtomMap.get(key)));
// set bond length
Molecule returnMolecule = rotatedMolecule.setDistance(atom1index, atom2index, currentLength);
//System.out.println(returnMolecule.getAngle(atom1alphaNumber,atom1number,atom1betaNumber));
//System.out.println(returnMolecule.getAngle(atom1betaNumber,atom1number,atom2number));
//System.out.println(returnMolecule.getAngle(atom2number,atom1number,atom1alphaNumber));
return returnMolecule;
}
public Molecule set_sp2(Atom atom1, Atom atom2)
{
return set_sp2(atom1, atom2, true);
}
/**
* Returns the bonded neighbors of includeAtom. Does not include includeAtom itself.
* @param includeAtom the atom whose neighbors are to be searched
* @return the Atoms adjacent to includeAtom
*/
public Set<Atom> getAdjacentAtoms(Atom includeAtom)
{
Set<Atom> returnSet = new HashSet<Atom>();
if ( ! connectivity.containsVertex(includeAtom) )
throw new IllegalArgumentException("includeAtom must be within this connectivity graph!");
for (DefaultWeightedEdge e : connectivity.edgesOf(includeAtom))
{
returnSet.add(connectivity.getEdgeSource(e));
returnSet.add(connectivity.getEdgeTarget(e));
}
returnSet.remove(includeAtom);
return(returnSet);
}
/**
* Checks if atom1 and atom2 are more than two bonds apart.
* @param atom1 the first atom
* @param atom2 the second atom
* @return true if atom1 and atom2 are separated by three or more bonds
*/
public boolean areSeparated(Atom atom1, Atom atom2)
{
Set<Atom> atom1neighbors = getAdjacentAtoms(atom1);
// check if direct neighbors
if (atom1neighbors.contains(atom2))
return false;
Set<Atom> atom2neighbors = getAdjacentAtoms(atom2);
// check if geminal
for (Atom a : atom1neighbors)
if (atom2neighbors.contains(a))
return false;
return true;
}
/**
* Checks if atoms are too close in a molecule, given another molecule
* whose atoms we know are not too close. Intended for assessing the result
* of dihedral changes. Does not assume the molecules have the same composition.
* @param oldMolecule the molecule this molecule was modified from
* @return true if there is at least one atom that is too close to another atom
*/
public boolean tooClose(Molecule oldMolecule)
{
// the atoms that have not changed in the new peptide
ArrayList<Atom> oldAtoms = new ArrayList<>();
// the atoms that are new in newPeptide
ArrayList<Atom> newAtoms = new ArrayList<>();
// populate lists
for (Atom a : contents)
{
if ( oldMolecule.contents.contains(a) )
oldAtoms.add(a);
else
newAtoms.add(a);
}
// compare distances between old atoms and new atoms
for (Atom oldAtom : oldAtoms)
{
Vector3D oldPosition = oldAtom.position;
for (Atom newAtom : newAtoms)
{
Vector3D newPosition = newAtom.position;
double distance = Vector3D.distance(oldPosition,newPosition);
if ( distance < Settings.MINIMUM_INTERATOMIC_DISTANCE &&
connectivity.getEdge(oldAtom, newAtom) == null )
return true;
}
}
return false;
}
/**
* Checks if the atoms are too close in a molecule. The minimum distance
* is controlled by Settings.MINIMUM_DISTANCE. Computations are performed
* on the triangular distance matrix.
* @return true if there is at least one atom that is too close to another atom
*/
public boolean tooClose()
{
// compute the upper triangle of distance contacts
for (int i=0; i < contents.size(); i++)
{
Atom atom1 = contents.get(i);
Vector3D position1 = atom1.position;
for (int j=i+1; j < contents.size(); j++)
{
Atom atom2 = contents.get(j);
Vector3D position2 = atom2.position;
// ignores distances between directly connected atoms
if ( Vector3D.distance(position1, position2) < Settings.MINIMUM_INTERATOMIC_DISTANCE &&
connectivity.getEdge(atom1,atom2) == null )
return true;
}
}
return false;
}
/**
* Calculates the "distance" between two molecules by comparing their positions.
* @param molecule1 the first molecule
* @param molecule2 the second molecule
* @return the RMS average distance between the atoms in molecule1 and molecule2 in angstroms
*/
public static double calculateRMSD(Molecule molecule1, Molecule molecule2)
{
if ( molecule1.contents.size() != molecule2.contents.size() )
throw new IllegalArgumentException("molecules are not the same size");
double RMSD = 0.0;
for (int i=0; i < molecule1.contents.size(); i++)
{
Vector3D fromVector = molecule1.contents.get(i).position;
Vector3D toVector = molecule2.contents.get(i).position;
RMSD += Math.pow(Vector3D.distance(fromVector,toVector), 2);
}
RMSD = RMSD * ( 1.0 / molecule1.contents.size() );
RMSD = Math.sqrt(RMSD);
return RMSD;
}
/**
* Calculates the "distance" between two molecules by comparing their positions.
* Molecules are superimposed first.
* @param molecule1 the first molecule
* @param molecule2 the second molecule
* @param atomIndices the indices (0,1,...,n-1) of the atoms to be compared
* @return the RMS average distance between the superimposed molecules, counting only the atoms in the lsit
*/
public static double calculateRMSD(Molecule molecule1, Molecule molecule2, List<Integer> atomIndices)
{
Molecule m2 = superimpose(molecule1,molecule2,atomIndices);
double RMSD = 0.0;
for (Integer atomIndex : atomIndices)
{
Atom fromAtom = molecule1.contents.get(atomIndex);
Atom toAtom = m2.contents.get(atomIndex);
Vector3D fromVector = fromAtom.position;
Vector3D toVector = toAtom.position;
RMSD += Math.pow(Vector3D.distance(fromVector,toVector), 2);
}
RMSD = RMSD * ( 1.0 / molecule1.contents.size() );
RMSD = Math.sqrt(RMSD);
return RMSD;
}
/**
* Superimposes two molecules based on a list of atom numbers.
* Superposition is performed using the
* <a href="http://en.wikipedia.org/wiki/Kabsch_algorithm">Kabsch algorithm</a>.
* @param molecule1 one of the molecules
* @param molecule2 another molecule (superposition transformation will be applied to this vector)
* @param atomIndices the atom indices to use for the superposition
* @return a new Molecule, which is molecule2 rotated to be maximally superimposed with molecule1
*/
public static Molecule superimpose(Molecule molecule1, Molecule molecule2, List<Integer> atomIndices)
{
// check validity of molecule arguments
if ( molecule1.contents.size() != molecule2.contents.size() )
throw new IllegalArgumentException("molecule size mismatch");
// check validity of atom numbers list
// does not check the validity of the numbers themselves
if ( atomIndices == null || atomIndices.size() < 3 || atomIndices.size() > molecule1.contents.size() )
throw new IllegalArgumentException("invalid atom numbers list");
// collect centroids
Vector3D centroid1 = molecule1.getCentroid();
Vector3D centroid2 = molecule2.getCentroid();
// move molecules to origin
Molecule m1 = molecule1.normalize();
Molecule m2 = molecule2.normalize();