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
|
package SmilesScripts::Chirality;
use strict;
use warnings;
use parent Exporter::;
our @EXPORT = qw(
chiral_volume
order_points_around_center
square_planar
trigonal_bipyramidal
);
use List::Util qw( all sum );
use Math::MatrixReal;
use Scalar::Util qw( blessed );
my $PI = 4 * atan2( 1, 1 );
sub chiral_volume
{
my @coords = @_[0..3];
my $normalise_vectors = $_[4];
if( all { blessed $_ && $_->isa( 'Chemistry::Atom' ) } @coords ) {
@coords = map { $_->coords } @coords;
} elsif( all { blessed $_ && $_->isa( 'Math::VectorReal' ) } @coords ) {
# No need to do anything
} else {
die "chiral_volume() called with parameters of unknown type\n";
}
my( $a, $b, $c ) = map { $_ - $coords[0] } @coords[1..3];
( $a, $b, $c ) = map { $_->norm } ( $a, $b, $c ) if $normalise_vectors;
return $a . ( $b x $c );
}
# Project the given 3D points to a minimal distance plane.
# Return the order of points in which they are "encountered" around their center.
sub order_points_around_center
{
my( @points3D ) = @_;
my( $A, $B, $C, $D ) = plane_equation( @points3D );
my @points2D = project_3D_to_plane( $A, $B, $C, $D, @points3D );
my $reference = shift @points2D; # first point is the reference
my @angles = map { $_ >= 0 ? $_ : $_ + 2*$PI }
map { atan2 $_->[1] - $reference->[1], $_->[0] - $reference->[0] }
@points2D;
return sort { $angles[$a] <=> $angles[$b] } 0..$#angles;
}
sub plane_equation
{
my( @points ) = @_;
my( $A, $B, $C, $D );
if( @points == 3 ) {
# From https://math.stackexchange.com/questions/2686606/equation-of-a-plane-passing-through-3-points
$A = $points[0]->[1] * $points[1]->[2] - $points[1]->[2] * $points[2]->[1];
$B = $points[0]->[2] * $points[2]->[0] - $points[0]->[0] * $points[2]->[2];
$C = $points[1]->[0] * $points[2]->[1] - $points[0]->[1] * $points[1]->[0];
$D = $points[0]->[0] * $points[1]->[1] * $points[2]->[2] -
$points[0]->[2] * $points[1]->[1] * $points[2]->[0];
} else {
# Calculate the centroid of points
my @centroid = ( sum( map { $_->[0] } @points ) / @points,
sum( map { $_->[1] } @points ) / @points,
sum( map { $_->[2] } @points ) / @points );
# Calculate the covariance matrix
my $xx = sum map { ($_->[0] - $centroid[0])**2 } @points;
my $yy = sum map { ($_->[1] - $centroid[1])**2 } @points;
my $zz = sum map { ($_->[2] - $centroid[2])**2 } @points;
my $xy = sum map { ($_->[0] - $centroid[0]) * ($_->[1] - $centroid[1]) } @points;
my $xz = sum map { ($_->[0] - $centroid[0]) * ($_->[2] - $centroid[2]) } @points;
my $yz = sum map { ($_->[1] - $centroid[1]) * ($_->[2] - $centroid[2]) } @points;
my @matrix = ( [ $xx, $xy, $xz ],
[ $xy, $yy, $yz ],
[ $xz, $yz, $zz ] );
my $matrix = Math::MatrixReal->new_from_rows( \@matrix );
# Find eigenvalues and eigenvectors, select smallest eigenvalue
my( $l, $V ) = $matrix->sym_diagonalize;
my $min = 1;
for (2..3) {
$min = $_ if $l->element( $min, 1 ) > $l->element( $_, 1 );
}
# Calculate the plane coefficients
( $A, $B, $C ) = @{$V->row( $min )->[0][0]};
$D = $A * $centroid[0] + $B * $centroid[1] + $C * $centroid[2];
}
return ( $A, $B, $C, $D );
}
sub project_3D_to_plane
{
my $A = shift;
my $B = shift;
my $C = shift;
my $D = shift;
my @points = @_;
return map { [ $_->[0] - $A * ( $A * $_->[0] + $B * $_->[1] + $C * $_->[2] + $D ) / ( $A**2 + $B**2 + $C**2 ),
$_->[1] - $B * ( $A * $_->[0] + $B * $_->[1] + $C * $_->[2] + $D ) / ( $A**2 + $B**2 + $C**2 ) ] }
@points;
}
sub square_planar
{
my @points3D = @_;
my @order = order_points_around_center( @points3D );
while( $order[0] != 0 ) {
push @order, shift @order;
}
return '@SP1' if $order[2] == 2;
return '@SP2' if $order[2] == 1;
return '@SP3' if $order[2] == 3;
}
sub trigonal_bipyramidal
{
my @points3D = @_;
my @order = order_points_around_center( @points3D );
while( $order[0] != 0 ) {
push @order, shift @order;
}
return $order[1] == 1 ? '@TB1' : '@TB2';
}
1;
|
