1
/**********************************************************************
2
pointgroup.cpp - Brute force symmetry analyzer.
4
(C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca
5
Some portions Copyright (C) 2007 by Geoffrey R. Hutchison
6
(Ported to C++, integrated with Open Babel)
8
This file is part of the Open Babel project.
9
For more information, see <http://openbabel.sourceforge.net/>
11
This program is free software; you can redistribute it and/or modify
12
it under the terms of the GNU General Public License as published by
13
the Free Software Foundation version 2 of the License.
15
This program is distributed in the hope that it will be useful,
16
but WITHOUT ANY WARRANTY; without even the implied warranty of
17
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
GNU General Public License for more details.
19
***********************************************************************/
21
#include <openbabel/babelconfig.h>
23
#include <openbabel/mol.h>
24
#include <openbabel/atom.h>
25
#include <openbabel/pointgroup.h>
32
#define M_PI 3.1415926535897932384626433832795028841971694
41
* All specific structures should have corresponding elements in the
42
* same position generic structure does.
44
* Planes are characterized by the surface normal direction
45
* (taken in the direction *from* the coordinate origin)
46
* and distance from the coordinate origin to the plane
47
* in the direction of the surface normal.
49
* Inversion is characterized by location of the inversion center.
51
* Rotation is characterized by a vector (distance+direction) from the origin
52
* to the rotation axis, axis direction and rotation order. Rotations
53
* are in the clockwise direction looking opposite to the direction
54
* of the axis. Note that this definition of the rotation axis
55
* is *not* unique, since an arbitrary multiple of the axis direction
56
* can be added to the position vector without changing actual operation.
58
* Mirror rotation is defined by the same parameters as normal rotation,
59
* but the origin is now unambiguous since it defines the position of the
60
* plane associated with the axis.
65
const char * group_name ; /* Canonical group name */
66
const char * symmetry_code ; /* Group symmetry code */
67
int (*check)( void ) ; /* Additional verification routine, not used */
71
* Point groups I know about
73
POINT_GROUP PointGroups[] = {
79
{ "C4", "(C4) (C2) "},
81
{ "C6", "(C6) (C3) (C2) "},
83
{ "C8", "(C8) (C4) (C2) "},
85
{ "D3", "(C3) 3*(C2) "},
86
{ "D4", "(C4) 5*(C2) "},
87
{ "D5", "(C5) 5*(C2) "},
88
{ "D6", "(C6) (C3) 7*(C2) "},
89
{ "D7", "(C7) 7*(C2) "},
90
{ "D8", "(C8) (C4) 9*(C2) "},
91
{ "C2v", "(C2) 2*(sigma) "},
92
{ "C3v", "(C3) 3*(sigma) "},
93
{ "C4v", "(C4) (C2) 4*(sigma) "},
94
{ "C5v", "(C5) 5*(sigma) "},
95
{ "C6v", "(C6) (C3) (C2) 6*(sigma) "},
96
{ "C7v", "(C7) 7*(sigma) "},
97
{ "C8v", "(C8) (C4) (C2) 8*(sigma) "},
98
{ "C2h", "(i) (C2) (sigma) "},
99
{ "C3h", "(C3) (S3) (sigma) "},
100
{ "C4h", "(i) (C4) (C2) (S4) (sigma) "},
101
{ "C5h", "(C5) (S5) (sigma) "},
102
{ "C6h", "(i) (C6) (C3) (C2) (S6) (S3) (sigma) "},
103
{ "C7h", "(C7) (S7) (sigma) "},
104
{ "C8h", "(i) (C8) (C4) (C2) (S8) (S4) (sigma) "},
105
{ "D2h", "(i) 3*(C2) 3*(sigma) "},
106
{ "D3h", "(C3) 3*(C2) (S3) 4*(sigma) "},
107
{ "D4h", "(i) (C4) 5*(C2) (S4) 5*(sigma) "},
108
{ "D5h", "(C5) 5*(C2) (S5) 6*(sigma) "},
109
{ "D6h", "(i) (C6) (C3) 7*(C2) (S6) (S3) 7*(sigma) "},
110
{ "D7h", "(C7) 7*(C2) (S7) 8*(sigma) "},
111
{ "D8h", "(i) (C8) (C4) 9*(C2) (S8) (S4) 9*(sigma) "},
112
{ "D2d", "3*(C2) (S4) 2*(sigma) "},
113
{ "D3d", "(i) (C3) 3*(C2) (S6) 3*(sigma) "},
114
{ "D4d", "(C4) 5*(C2) (S8) 4*(sigma) "},
115
{ "D5d", "(i) (C5) 5*(C2) (S10) 5*(sigma) "},
116
{ "D6d", "(C6) (C3) 7*(C2) (S12) (S4) 6*(sigma) "},
117
{ "D7d", "(i) (C7) 7*(C2) (S14) 7*(sigma) "},
118
{ "D8d", "(C8) (C4) 9*(C2) (S16) 8*(sigma) "},
119
{ "S4", "(C2) (S4) "},
120
{ "S6", "(i) (C3) (S6) "},
121
{ "S8", "(C4) (C2) (S8) "},
122
{ "T", "4*(C3) 3*(C2) "},
123
{ "Th", "(i) 4*(C3) 3*(C2) 4*(S6) 3*(sigma) "},
124
{ "Td", "4*(C3) 3*(C2) 3*(S4) 6*(sigma) "},
125
{ "O", "3*(C4) 4*(C3) 9*(C2) "},
126
{ "Oh", "(i) 3*(C4) 4*(C3) 9*(C2) 4*(S6) 3*(S4) 9*(sigma) "},
127
{ "Cinfv", "(Cinf) (sigma) "},
128
{ "Dinfh", "(i) (Cinf) (C2) 2*(sigma) "},
129
{ "I", "6*(C5) 10*(C3) 15*(C2) "},
130
{ "Ih", "(i) 6*(C5) 10*(C3) 15*(C2) 6*(S10) 10*(S6) 15*(sigma) "},
131
{ "Kh", "(i) (Cinf) (sigma) "},
133
#define PointGroupsCount (sizeof(PointGroups)/sizeof(POINT_GROUP))
135
class PointGroupPrivate
138
struct _SYMMETRY_ELEMENT_ {
139
void (*transform_atom)( struct _SYMMETRY_ELEMENT_ *el, OBAtom *from, OBAtom *to ) ;
140
int * transform ; /* Correspondence table for the transformation */
141
int order ; /* Applying transformation this many times is identity */
142
int nparam ; /* 4 for inversion and planes, 7 for axes */
143
double maxdev ; /* Largest error associated with the element */
145
double normal[ DIMENSION ] ;
146
double direction[ DIMENSION ] ;
149
typedef _SYMMETRY_ELEMENT_ SYMMETRY_ELEMENT;
153
ToleranceSame = 1e-3;
154
TolerancePrimary = 1e-1;
155
ToleranceFinal = 1e-3;
159
OptChangeThreshold = 1e-10;
160
DistanceFromCenter = NULL;
167
MolecularPlane = NULL;
168
InversionCentersCount = 0;
169
InversionCenters = NULL;
172
ImproperAxesCount = 0;
174
NormalAxesCounts = NULL;
175
ImproperAxesCounts = NULL;
178
PointGroupRejectionReason = NULL ;
190
double ToleranceSame ;
191
double TolerancePrimary ;
192
double ToleranceFinal ;
195
double GradientStep ;
196
double OptChangeThreshold ;
197
double CenterOfSomething[ DIMENSION ];
198
double * DistanceFromCenter ;
204
SYMMETRY_ELEMENT ** Planes ;
205
SYMMETRY_ELEMENT * MolecularPlane ;
206
int InversionCentersCount ;
207
SYMMETRY_ELEMENT ** InversionCenters ;
208
int NormalAxesCount ;
209
SYMMETRY_ELEMENT ** NormalAxes ;
210
int ImproperAxesCount ;
211
SYMMETRY_ELEMENT ** ImproperAxes ;
212
int * NormalAxesCounts ;
213
int * ImproperAxesCounts ;
214
int BadOptimization ;
215
const char * SymmetryCode ;
216
char * PointGroupRejectionReason;
229
bool equivalentAtoms(OBAtom &a1, OBAtom &a2)
231
if (a1.GetAtomicNum() != a2.GetAtomicNum())
233
if (a1.GetIsotope() != a2.GetIsotope())
235
if (a1.GetFormalCharge() != a2.GetFormalCharge())
237
if (a1.GetSpinMultiplicity() != a2.GetSpinMultiplicity())
244
establish_pairs( SYMMETRY_ELEMENT *elem )
247
char * atom_used = (char *)calloc( _mol->NumAtoms(), 1 ) ;
248
double distance, best_distance ;
252
if( atom_used == NULL ){
253
fprintf( stderr, "Out of memory for tagging array in establish_pairs()\n" ) ;
254
exit( EXIT_FAILURE ) ;
256
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
257
if( elem->transform[i] >= _mol->NumAtoms() ){ /* No symmetric atom yet */
258
if( verbose > 2 ) printf( " looking for a pair for %d\n", i ) ;
259
elem->transform_atom( elem, _mol->GetAtom(i+1), &symmetric ) ; // ATOM INDEX ISSUE
260
if( verbose > 2 ) printf( " new coordinates are: (%g,%g,%g)\n",
261
symmetric.x(), symmetric.y(), symmetric.z() ) ;
263
best_distance = 2*TolerancePrimary ;/* Performance value we'll reject */
264
for( j = 0 ; j < _mol->NumAtoms() ; j++ ){
266
atom = _mol->GetAtom(j+1);
268
if( atom_used[j] || !equivalentAtoms(*atom, symmetric) )
271
distance = symmetric.GetDistance(atom) ;
272
if( verbose > 2 ) printf( " distance to %d is %g\n", j, distance ) ;
273
if( distance < best_distance ){
275
best_distance = distance ;
278
if( best_distance > TolerancePrimary ){ /* Too bad, there is no symmetric atom */
280
printf( " no pair for atom %d - best was %d with err = %g\n", i, best_j, best_distance ) ;
284
elem->transform[i] = best_j ;
285
atom_used[best_j] = 1 ;
286
if( verbose > 1 ) printf( " atom %d transforms to the atom %d, err = %g\n", i, best_j, best_distance ) ;
294
check_transform_order( SYMMETRY_ELEMENT *elem )
298
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
299
if( elem->transform[i] == i ) /* Identity transform is Ok for any order */
301
if( elem->transform_atom == rotate_reflect_atom ){
302
j = elem->transform[i] ;
303
if( elem->transform[j] == i )
304
continue ; /* Second-order transform is Ok for improper axis */
306
for( j = elem->order - 1, k = elem->transform[i] ; j > 0 ; j--, k = elem->transform[k] ){
308
if( verbose > 0 ) printf( " transform looped %d steps too early from atom %d\n", j, i ) ;
312
if( k != i && elem->transform_atom == rotate_reflect_atom ){
313
/* For improper axes, the complete loop may also take twice the order */
314
for( j = elem->order ; j > 0 ; j--, k = elem->transform[k] ){
316
if( verbose > 0 ) printf( " (improper) transform looped %d steps too early from atom %d\n", j, i ) ;
322
if( verbose > 0 ) printf( " transform failed to loop after %d steps from atom %d\n", elem->order, i ) ;
330
same_transform( SYMMETRY_ELEMENT *a, SYMMETRY_ELEMENT *b )
335
if( ( a->order != b->order ) || ( a->nparam != b->nparam ) || ( a->transform_atom != b->transform_atom ) )
337
for( i = 0, code = 1 ; i < _mol->NumAtoms() ; i++ ){
338
if( a->transform[i] != b->transform[i] ){
343
if( code == 0 && a->order > 2 ){ /* b can also be a reverse transformation for a */
344
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
345
j = a->transform[i] ;
346
if( b->transform[j] != i )
355
alloc_symmetry_element( void )
357
SYMMETRY_ELEMENT * elem = (SYMMETRY_ELEMENT *)calloc( 1, sizeof( SYMMETRY_ELEMENT ) ) ;
361
fprintf( stderr, "Out of memory allocating symmetry element\n" ) ;
362
exit( EXIT_FAILURE ) ;
364
elem->transform = (int*)calloc( _mol->NumAtoms(), sizeof( int ) ) ;
365
if( elem->transform == NULL ){
366
fprintf( stderr, "Out of memory allocating transform table for symmetry element\n" ) ;
367
exit( EXIT_FAILURE ) ;
369
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
370
elem->transform[i] = _mol->NumAtoms() + 1 ; /* An impossible value */
376
destroy_symmetry_element( SYMMETRY_ELEMENT *elem )
379
if( elem->transform != NULL )
380
free( elem->transform ) ;
386
check_transform_quality( SYMMETRY_ELEMENT *elem )
392
for( i = 0, max_r = 0 ; i < _mol->NumAtoms() ; i++ ){
393
j = elem->transform[i] ;
394
elem->transform_atom( elem, _mol->GetAtom(i+1), &symmetric ) ;
396
r = symmetric.GetDistance(_mol->GetAtom(j+1));
397
if( r > ToleranceFinal ){
398
if( verbose > 0 ) printf( " distance to symmetric atom (%g) is too big for %d\n", r, i ) ;
401
if( r > max_r ) max_r = r ;
403
elem->maxdev = max_r ;
408
eval_optimization_target_function( SYMMETRY_ELEMENT *elem, int *finish )
412
double target, r, maxr ;
414
if( elem->nparam >= 4 ){
415
for( k = 0, r = 0 ; k < DIMENSION ; k++ ){
416
r += elem->normal[k]*elem->normal[k] ;
419
if( r < ToleranceSame ){
420
fprintf( stderr, "Normal collapced!\n" ) ;
421
exit( EXIT_FAILURE ) ;
423
for( k = 0 ; k < DIMENSION ; k++ ){
424
elem->normal[k] /= r ;
426
if( elem->distance < 0 ){
427
elem->distance = -elem->distance ;
428
for( k = 0 ; k < DIMENSION ; k++ ){
429
elem->normal[k] = -elem->normal[k] ;
433
if( elem->nparam >= 7 ){
434
for( k = 0, r = 0 ; k < DIMENSION ; k++ ){
435
r += elem->direction[k]*elem->direction[k] ;
438
if( r < ToleranceSame ){
439
fprintf( stderr, "Direction collapced!\n" ) ;
440
exit( EXIT_FAILURE ) ;
442
for( k = 0 ; k < DIMENSION ; k++ ){
443
elem->direction[k] /= r ;
446
for( i = 0, target = maxr = 0 ; i < _mol->NumAtoms() ; i++ ){
447
elem->transform_atom( elem, _mol->GetAtom(i+1), &symmetric ) ;
448
j = elem->transform[i] ;
449
r = symmetric.GetDistance(_mol->GetAtom(j+1));
450
if( r > maxr ) maxr = r ;
453
if( finish != NULL ){
455
if( maxr < ToleranceFinal )
462
get_params( SYMMETRY_ELEMENT *elem, double values[] )
464
memcpy( values, &elem->distance, elem->nparam * sizeof( double ) ) ;
468
set_params( SYMMETRY_ELEMENT *elem, double values[] )
470
memcpy( &elem->distance, values, elem->nparam * sizeof( double ) ) ;
474
optimize_transformation_params( SYMMETRY_ELEMENT *elem )
476
double values[ MAXPARAM ] ;
477
double grad [ MAXPARAM ] ;
478
double force [ MAXPARAM ] ;
479
double step [ MAXPARAM ] ;
480
double f, fold, fnew, fnew2, fdn, fup, snorm ;
482
int vars = elem->nparam ;
487
if( vars > MAXPARAM ){
488
fprintf( stderr, "Catastrophe in optimize_transformation_params()!\n" ) ;
489
exit( EXIT_FAILURE ) ;
494
f = eval_optimization_target_function( elem, &finish ) ;
495
/* Evaluate function, gradient and diagonal force constants */
496
if( verbose > 1 ) printf( " function value = %g\n", f ) ;
498
if( verbose > 1 ) printf( " function value is small enough\n" ) ;
502
if( fabs( f-fold ) > OptChangeThreshold )
505
if( hits >= OptChangeHits ){
506
if( verbose > 1 ) printf( " no progress is made, stop optimization\n" ) ;
510
get_params( elem, values ) ;
511
for( i = 0 ; i < vars ; i++ ){
512
values[i] -= GradientStep ;
513
set_params( elem, values ) ;
514
fdn = eval_optimization_target_function( elem, NULL ) ;
515
values[i] += 2*GradientStep ;
516
set_params( elem, values ) ;
517
fup = eval_optimization_target_function( elem, NULL ) ;
518
values[i] -= GradientStep ;
519
grad[i] = ( fup - fdn ) / ( 2 * GradientStep ) ;
520
force[i] = ( fup + fdn - 2*f ) / ( GradientStep * GradientStep ) ;
521
if( verbose > 1 ) printf( " i = %d, grad = %12.6e, force = %12.6e\n", i, grad[i], force[i] ) ;
523
/* Do a quasy-Newton step */
524
for( i = 0, snorm = 0 ; i < vars ; i++ ){
525
if( force[i] < 0 ) force[i] = -force[i] ;
526
if( force[i] < 1e-3 ) force[i] = 1e-3 ;
527
if( force[i] > 1e3 ) force[i] = 1e3 ;
528
step[i] = - grad[i]/force[i] ;
529
snorm += step[i] * step[i] ;
531
snorm = sqrt( snorm ) ;
532
if( snorm > MaxOptStep ){ /* Renormalize step */
533
for( i = 0 ; i < vars ; i++ )
534
step[i] *= MaxOptStep/snorm ;
538
for( i = 0 ; i < vars ; i++ ){
539
values[i] += step[i] ;
541
set_params( elem, values ) ;
542
fnew = eval_optimization_target_function( elem, NULL ) ;
545
for( i = 0 ; i < vars ; i++ ){
546
values[i] -= step[i] ;
549
set_params( elem, values ) ;
551
} while( snorm > MinOptStep ) ;
552
if( (snorm > MinOptStep) && (snorm < MaxOptStep / 2) ){ /* try to do quadratic interpolation */
553
for( i = 0 ; i < vars ; i++ )
554
values[i] += step[i] ;
555
set_params( elem, values ) ;
556
fnew2 = eval_optimization_target_function( elem, NULL ) ;
557
if( verbose > 1 ) printf( " interpolation base points: %g, %g, %g\n", f, fnew, fnew2 ) ;
558
for( i = 0 ; i < vars ; i++ )
559
values[i] -= 2*step[i] ;
560
a = ( 4*f - fnew2 - 3*fnew ) / 2 ;
561
b = ( f + fnew2 - 2*fnew ) / 2 ;
562
if( verbose > 1 ) printf( " linear interpolation coefficients %g, %g\n", a, b ) ;
565
if( x > 0.2 && x < 1.8 ){
566
if( verbose > 1 ) printf( " interpolated: %g\n", x ) ;
567
for( i = 0 ; i < vars ; i++ )
568
values[i] += x*step[i] ;
574
for( i = 0 ; i < vars ; i++ )
575
values[i] += 2*step[i] ;
578
for( i = 0 ; i < vars ; i++ )
579
values[i] += step[i] ;
582
set_params( elem, values ) ;
584
} while( snorm > MinOptStep && ++cycle < MaxOptCycles ) ;
585
f = eval_optimization_target_function( elem, NULL ) ;
586
if( cycle >= MaxOptCycles ) BadOptimization = 1 ;
588
if( cycle >= MaxOptCycles )
589
printf( " maximum number of optimization cycles made\n" ) ;
590
printf( " optimization completed after %d cycles with f = %g\n", cycle, f ) ;
595
refine_symmetry_element( SYMMETRY_ELEMENT *elem, int build_table )
600
if( build_table && (establish_pairs( elem ) < 0) ){
602
if( verbose > 0 ) printf( " no transformation correspondence table can be constructed\n" ) ;
605
for( i = 0 ; i < PlanesCount ; i++ ){
606
if( same_transform( Planes[i], elem ) ){
608
if( verbose > 0 ) printf( " transformation is identical to plane %d\n", i ) ;
612
for( i = 0 ; i < InversionCentersCount ; i++ ){
613
if( same_transform( InversionCenters[i], elem ) ){
615
if( verbose > 0 ) printf( " transformation is identical to inversion center %d\n", i ) ;
619
for( i = 0 ; i < NormalAxesCount ; i++ ){
620
if( same_transform( NormalAxes[i], elem ) ){
622
if( verbose > 0 ) printf( " transformation is identical to normal axis %d\n", i ) ;
626
for( i = 0 ; i < ImproperAxesCount ; i++ ){
627
if( same_transform( ImproperAxes[i], elem ) ){
629
if( verbose > 0 ) printf( " transformation is identical to improper axis %d\n", i ) ;
633
if( check_transform_order( elem ) < 0 ){
635
if( verbose > 0 ) printf( " incorrect transformation order\n" ) ;
638
optimize_transformation_params( elem ) ;
639
if( check_transform_quality( elem ) < 0 ){
641
if( verbose > 0 ) printf( " refined transformation does not pass the numeric threshold\n" ) ;
649
* Plane-specific functions
652
static void mirror_atom( SYMMETRY_ELEMENT *plane, OBAtom *from, OBAtom *to )
654
double r = plane->distance;
656
r -= from->x() * plane->normal[0];
657
r -= from->y() * plane->normal[1];
658
r -= from->z() * plane->normal[2];
660
// copy the "type" of from into to
661
to->SetAtomicNum(from->GetAtomicNum());
662
to->SetIsotope(from->GetIsotope());
663
to->SetFormalCharge(from->GetFormalCharge());
664
to->SetSpinMultiplicity(from->GetSpinMultiplicity());
666
to->SetVector(from->x() + 2*r*plane->normal[0],
667
from->y() + 2*r*plane->normal[1],
668
from->z() + 2*r*plane->normal[2]);
672
init_mirror_plane( int i, int j )
674
SYMMETRY_ELEMENT * plane = alloc_symmetry_element() ;
675
double dx[ DIMENSION ], midpoint[ DIMENSION ], rab, r ;
678
if( verbose > 0 ) printf( "Trying mirror plane for atoms %d,%d\n", i, j ) ;
680
plane->transform_atom = mirror_atom;
684
dx[0] = _mol->GetAtom(i+1)->x() - _mol->GetAtom(j+1)->x();
685
dx[1] = _mol->GetAtom(i+1)->y() - _mol->GetAtom(j+1)->y();
686
dx[2] = _mol->GetAtom(i+1)->z() - _mol->GetAtom(j+1)->z();
688
midpoint[0] = ( _mol->GetAtom(i+1)->x() + _mol->GetAtom(j+1)->x() ) / 2.0 ;
689
midpoint[1] = ( _mol->GetAtom(i+1)->y() + _mol->GetAtom(j+1)->y() ) / 2.0 ;
690
midpoint[2] = ( _mol->GetAtom(i+1)->z() + _mol->GetAtom(j+1)->z() ) / 2.0 ;
692
rab = _mol->GetAtom(i+1)->GetDistance(_mol->GetAtom(j+1));
694
if( rab < ToleranceSame ){
695
fprintf( stderr, "Atoms %d and %d coincide (r = %g)\n", i, j, rab ) ;
696
exit( EXIT_FAILURE ) ;
698
for( k = 0, r = 0 ; k < DIMENSION ; k++ ){
699
plane->normal[k] = dx[k]/rab ;
700
r += midpoint[k]*plane->normal[k] ;
702
if( r < 0 ){ /* Reverce normal direction, distance is always positive! */
704
for( k = 0 ; k < DIMENSION ; k++ ){
705
plane->normal[k] = -plane->normal[k] ;
708
plane->distance = r ;
709
if( verbose > 0 ) printf( " initial plane is at %g from the origin\n", r ) ;
710
if( refine_symmetry_element( plane, 1 ) < 0 ){
711
if( verbose > 0 ) printf( " refinement failed for the plane\n" ) ;
712
destroy_symmetry_element( plane ) ;
719
init_ultimate_plane( void )
721
SYMMETRY_ELEMENT * plane = alloc_symmetry_element() ;
722
double d0[ DIMENSION ], d1[ DIMENSION ], d2[ DIMENSION ] ;
723
double p[ DIMENSION ] ;
724
double r, s0, s1, s2 ;
728
if( verbose > 0 ) printf( "Trying whole-molecule mirror plane\n" ) ;
730
plane->transform_atom = mirror_atom ;
733
for( k = 0 ; k < DIMENSION ; k++ )
734
d0[k] = d1[k] = d2[k] = 0 ;
735
d0[0] = 1 ; d1[1] = 1 ; d2[2] = 1 ;
736
for( i = 1 ; i < _mol->NumAtoms() ; i++ ){
737
for( j = 0 ; j < i ; j++ ){
738
p[0] = _mol->GetAtom(i+1)->x() - _mol->GetAtom(j+1)->x();
739
p[1] = _mol->GetAtom(i+1)->y() - _mol->GetAtom(j+1)->y();
740
p[2] = _mol->GetAtom(i+1)->z() - _mol->GetAtom(j+1)->z();
742
r = sqrt(SQUARE(p[0]) + SQUARE(p[1]) + SQUARE(p[2])); // distance between atoms i and j
744
for( k = 0, s0=s1=s2=0 ; k < DIMENSION ; k++ ){
750
for( k = 0 ; k < DIMENSION ; k++ ){
757
for( k = 0, s0=s1=s2=0 ; k < DIMENSION ; k++ ){
763
if( s0 >= s1 && s0 >= s2 ) d = d0 ;
764
if( s1 >= s0 && s1 >= s2 ) d = d1 ;
765
if( s2 >= s0 && s2 >= s1 ) d = d2 ;
767
fprintf( stderr, "Catastrophe in init_ultimate_plane(): %g, %g and %g have no ordering!\n", s0, s1, s2 ) ;
768
exit( EXIT_FAILURE ) ;
770
for( k = 0, r = 0 ; k < DIMENSION ; k++ )
774
for( k = 0 ; k < DIMENSION ; k++ )
775
plane->normal[k] = d[k]/r ;
778
for( k = 1 ; k < DIMENSION ; k++ )
779
plane->normal[k] = 0 ;
780
plane->normal[0] = 1 ;
782
for( k = 0, r = 0 ; k < DIMENSION ; k++ )
783
r += CenterOfSomething[k]*plane->normal[k] ;
784
plane->distance = r ;
785
for( k = 0 ; k < _mol->NumAtoms() ; k++ )
786
plane->transform[k] = k ;
787
if( refine_symmetry_element( plane, 0 ) < 0 ){
788
if( verbose > 0 ) printf( " refinement failed for the plane\n" ) ;
789
destroy_symmetry_element( plane ) ;
796
* Inversion-center specific functions
799
invert_atom( SYMMETRY_ELEMENT *center, OBAtom *from, OBAtom *to )
802
// copy the "type" of from into to
803
to->SetAtomicNum(from->GetAtomicNum());
804
to->SetIsotope(from->GetIsotope());
805
to->SetFormalCharge(from->GetFormalCharge());
806
to->SetSpinMultiplicity(from->GetSpinMultiplicity());
808
to->SetVector(2*center->distance*center->normal[0] - from->x(),
809
2*center->distance*center->normal[1] - from->y(),
810
2*center->distance*center->normal[2] - from->z());
814
init_inversion_center( void )
816
SYMMETRY_ELEMENT * center = alloc_symmetry_element() ;
820
if( verbose > 0 ) printf( "Trying inversion center at the center of something\n" ) ;
822
center->transform_atom = invert_atom ;
825
for( k = 0, r = 0 ; k < DIMENSION ; k++ )
826
r += CenterOfSomething[k]*CenterOfSomething[k] ;
829
for( k = 0 ; k < DIMENSION ; k++ )
830
center->normal[k] = CenterOfSomething[k]/r ;
833
center->normal[0] = 1 ;
834
for( k = 1 ; k < DIMENSION ; k++ )
835
center->normal[k] = 0 ;
837
center->distance = r ;
838
if( verbose > 0 ) printf( " initial inversion center is at %g from the origin\n", r ) ;
839
if( refine_symmetry_element( center, 1 ) < 0 ){
840
if( verbose > 0 ) printf( " refinement failed for the inversion center\n" ) ;
841
destroy_symmetry_element( center ) ;
848
* Normal rotation axis-specific routines.
851
rotate_atom( SYMMETRY_ELEMENT *axis, OBAtom *from, OBAtom *to )
853
double x[3], y[3], a[3], b[3], c[3] ;
854
double angle = axis->order ? 2*M_PI/axis->order : 1.0 ;
855
double a_sin = sin( angle ) ;
856
double a_cos = cos( angle ) ;
860
if( DIMENSION != 3 ){
861
fprintf( stderr, "Catastrophe in rotate_atom!\n" ) ;
862
exit( EXIT_FAILURE ) ;
865
x[0] = from->x() - axis->distance * axis->normal[0];
866
x[1] = from->y() - axis->distance * axis->normal[1];
867
x[2] = from->z() - axis->distance * axis->normal[2];
869
for( i = 0, dot = 0 ; i < 3 ; i++ )
870
dot += x[i] * axis->direction[i] ;
871
for( i = 0 ; i < 3 ; i++ )
872
a[i] = axis->direction[i] * dot ;
873
for( i = 0 ; i < 3 ; i++ )
875
c[0] = b[1]*axis->direction[2] - b[2]*axis->direction[1] ;
876
c[1] = b[2]*axis->direction[0] - b[0]*axis->direction[2] ;
877
c[2] = b[0]*axis->direction[1] - b[1]*axis->direction[0] ;
878
for( i = 0 ; i < 3 ; i++ )
879
y[i] = a[i] + b[i]*a_cos + c[i]*a_sin ;
881
to->SetVector(y[0] + axis->distance * axis->normal[0],
882
y[1] + axis->distance * axis->normal[1],
883
y[2] + axis->distance * axis->normal[2]);
885
// copy the "type" of from into to
886
to->SetAtomicNum(from->GetAtomicNum());
887
to->SetIsotope(from->GetIsotope());
888
to->SetFormalCharge(from->GetFormalCharge());
889
to->SetSpinMultiplicity(from->GetSpinMultiplicity());
893
init_ultimate_axis(void)
895
SYMMETRY_ELEMENT * axis = alloc_symmetry_element() ;
896
double dir[ DIMENSION ], rel[ DIMENSION ] ;
900
if( verbose > 0 ) printf( "Trying infinity axis\n" ) ;
902
axis->transform_atom = rotate_atom ;
905
for( k = 0 ; k < DIMENSION ; k++ )
907
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
909
rel[0] = _mol->GetAtom(i+1)->x() - CenterOfSomething[0];
910
rel[1] = _mol->GetAtom(i+1)->z() - CenterOfSomething[1];
911
rel[2] = _mol->GetAtom(i+1)->y() - CenterOfSomething[2];
913
s = rel[0]*dir[0] + rel[1]*dir[1] + rel[2]*dir[2];
916
for( k = 0 ; k < DIMENSION ; k++ )
918
else for( k = 0 ; k < DIMENSION ; k++ )
921
for( k = 0, s = 0 ; k < DIMENSION ; k++ )
922
s += SQUARE( dir[k] ) ;
925
for( k = 0 ; k < DIMENSION ; k++ )
928
for( k = 0 ; k < DIMENSION ; k++ )
929
axis->direction[k] = dir[k] ;
930
for( k = 0, s = 0 ; k < DIMENSION ; k++ )
931
s += SQUARE( CenterOfSomething[k] ) ;
934
for( k = 0 ; k < DIMENSION ; k++ )
935
axis->normal[k] = CenterOfSomething[k]/s ;
937
for( k = 1 ; k < DIMENSION ; k++ )
938
axis->normal[k] = 0 ;
939
axis->normal[0] = 1 ;
942
for( k = 0 ; k < _mol->NumAtoms() ; k++ )
943
axis->transform[k] = k ;
944
if( refine_symmetry_element( axis, 0 ) < 0 ){
945
if( verbose > 0 ) printf( " refinement failed for the infinity axis\n" ) ;
946
destroy_symmetry_element( axis ) ;
954
init_c2_axis( int i, int j, const double support[ DIMENSION ] )
956
SYMMETRY_ELEMENT * axis ;
959
double r, center[ DIMENSION ] ;
962
printf( "Trying c2 axis for the pair (%d,%d) with the support (%g,%g,%g)\n",
963
i, j, support[0], support[1], support[2] ) ;
965
/* First, do a quick sanity check */
967
vector3 supportVec(support[0], support[1], support[2]);
968
ris = vector3(_mol->GetAtom(i+1)->GetVector() - supportVec).length();
969
rjs = vector3(_mol->GetAtom(j+1)->GetVector() - supportVec).length();
971
if( fabs( ris - rjs ) > TolerancePrimary ){
973
if( verbose > 0 ) printf( " Support can't actually define a rotation axis\n" ) ;
976
axis = alloc_symmetry_element() ;
977
axis->transform_atom = rotate_atom ;
980
for( k = 0, r = 0 ; k < DIMENSION ; k++ )
981
r += CenterOfSomething[k]*CenterOfSomething[k] ;
984
for( k = 0 ; k < DIMENSION ; k++ )
985
axis->normal[k] = CenterOfSomething[k]/r ;
988
axis->normal[0] = 1 ;
989
for( k = 1 ; k < DIMENSION ; k++ )
990
axis->normal[k] = 0 ;
994
center[0] = ( _mol->GetAtom(i+1)->x() + _mol->GetAtom(j+1)->x() ) / 2 - support[0] ;
995
center[1] = ( _mol->GetAtom(i+1)->y() + _mol->GetAtom(j+1)->y() ) / 2 - support[1] ;
996
center[2] = ( _mol->GetAtom(i+1)->z() + _mol->GetAtom(j+1)->z() ) / 2 - support[2] ;
997
r = sqrt(SQUARE(center[0]) + SQUARE(center[1]) + SQUARE(center[2]));
999
if( r <= TolerancePrimary ){ /* c2 is underdefined, let's do something special */
1000
if( MolecularPlane != NULL ){
1001
if( verbose > 0 ) printf( " c2 is underdefined, but there is a molecular plane\n" ) ;
1002
for( k = 0 ; k < DIMENSION ; k++ )
1003
axis->direction[k] = MolecularPlane->normal[k] ;
1006
if( verbose > 0 ) printf( " c2 is underdefined, trying random direction\n" ) ;
1008
center[0] = _mol->GetAtom(i+1)->x() - _mol->GetAtom(j+1)->x();
1009
center[1] = _mol->GetAtom(i+1)->y() - _mol->GetAtom(j+1)->y();
1010
center[2] = _mol->GetAtom(i+1)->z() - _mol->GetAtom(j+1)->z();
1012
if( fabs( center[2] ) + fabs( center[1] ) > ToleranceSame ){
1013
axis->direction[0] = 0 ;
1014
axis->direction[1] = center[2] ;
1015
axis->direction[2] = -center[1] ;
1018
axis->direction[0] = -center[2] ;
1019
axis->direction[1] = 0 ;
1020
axis->direction[2] = center[0] ;
1022
for( k = 0, r = 0 ; k < DIMENSION ; k++ )
1023
r += axis->direction[k] * axis->direction[k] ;
1025
for( k = 0 ; k < DIMENSION ; k++ )
1026
axis->direction[k] /= r ;
1029
else { /* direction is Ok, renormalize it */
1030
for( k = 0 ; k < DIMENSION ; k++ )
1031
axis->direction[k] = center[k]/r ;
1033
if( refine_symmetry_element( axis, 1 ) < 0 ){
1034
if( verbose > 0 ) printf( " refinement failed for the c2 axis\n" ) ;
1035
destroy_symmetry_element( axis ) ;
1042
init_axis_parameters( double a[3], double b[3], double c[3] )
1044
SYMMETRY_ELEMENT * axis ;
1045
int i, order, sign ;
1046
double ra, rb, rc, rab, rbc, rac, r ;
1049
ra = rb = rc = rab = rbc = rac = 0 ;
1050
for( i = 0 ; i < DIMENSION ; i++ ){
1055
ra = sqrt(ra) ; rb = sqrt(rb) ; rc = sqrt(rc) ;
1056
if( fabs( ra - rb ) > TolerancePrimary || fabs( ra - rc ) > TolerancePrimary || fabs( rb - rc ) > TolerancePrimary ){
1058
if( verbose > 0 ) printf( " points are not on a sphere\n" ) ;
1061
for( i = 0 ; i < DIMENSION ; i++ ){
1062
rab += (a[i]-b[i])*(a[i]-b[i]) ;
1063
rac += (a[i]-c[i])*(a[i]-c[i]) ;
1064
rbc += (c[i]-b[i])*(c[i]-b[i]) ;
1069
if( fabs( rab - rbc ) > TolerancePrimary ){
1071
if( verbose > 0 ) printf( " points can't be rotation-equivalent\n" ) ;
1074
if( rab <= ToleranceSame || rbc <= ToleranceSame || rac <= ToleranceSame ){
1076
if( verbose > 0 ) printf( " rotation is underdefined by these points\n" ) ;
1080
angle = M_PI - 2*asin( rac/(2*rab) ) ;
1081
if( verbose > 1 ) printf( " rotation angle is %f\n", angle ) ;
1082
if( fabs(angle) <= M_PI/(MaxAxisOrder+1) ){
1084
if( verbose > 0 ) printf( " atoms are too close to a straight line\n" ) ;
1087
order = floor( (2*M_PI)/angle + 0.5 ) ;
1088
if( order <= 2 || order > MaxAxisOrder ){
1090
if( verbose > 0 ) printf( " rotation axis order (%d) is not from 3 to %d\n", order, MaxAxisOrder ) ;
1093
axis = alloc_symmetry_element() ;
1094
axis->order = order ;
1096
for( i = 0, r = 0 ; i < DIMENSION ; i++ )
1097
r += CenterOfSomething[i]*CenterOfSomething[i] ;
1100
for( i = 0 ; i < DIMENSION ; i++ )
1101
axis->normal[i] = CenterOfSomething[i]/r ;
1104
axis->normal[0] = 1 ;
1105
for( i = 1 ; i < DIMENSION ; i++ )
1106
axis->normal[i] = 0 ;
1108
axis->distance = r ;
1109
axis->direction[0] = (b[1]-a[1])*(c[2]-b[2]) - (b[2]-a[2])*(c[1]-b[1]) ;
1110
axis->direction[1] = (b[2]-a[2])*(c[0]-b[0]) - (b[0]-a[0])*(c[2]-b[2]) ;
1111
axis->direction[2] = (b[0]-a[0])*(c[1]-b[1]) - (b[1]-a[1])*(c[0]-b[0]) ;
1113
* Arbitrarily select axis direction so that first non-zero component
1114
* or the direction is positive.
1117
if( axis->direction[0] <= 0 )
1118
if( axis->direction[0] < 0 )
1120
else if( axis->direction[1] <= 0 )
1121
if( axis->direction[1] < 0 )
1123
else if( axis->direction[2] < 0 )
1126
for( i = 0 ; i < DIMENSION ; i++ )
1127
axis->direction[i] = -axis->direction[i] ;
1128
for( i = 0, r = 0 ; i < DIMENSION ; i++ )
1129
r += axis->direction[i]*axis->direction[i] ;
1131
for( i = 0 ; i < DIMENSION ; i++ )
1132
axis->direction[i] /= r ;
1134
printf( " axis origin is at (%g,%g,%g)\n",
1135
axis->normal[0]*axis->distance, axis->normal[1]*axis->distance, axis->normal[2]*axis->distance ) ;
1136
printf( " axis is in the direction (%g,%g,%g)\n", axis->direction[0], axis->direction[1], axis->direction[2] ) ;
1142
init_higher_axis( int ia, int ib, int ic )
1144
SYMMETRY_ELEMENT * axis ;
1145
double a[ DIMENSION ], b[ DIMENSION ], c[ DIMENSION ] ;
1147
if( verbose > 0 ) printf( "Trying cn axis for the triplet (%d,%d,%d)\n", ia, ib, ic ) ;
1149
/* Do a quick check of geometry validity */
1151
a[0] = _mol->GetAtom(ia+1)->x() - CenterOfSomething[0];
1152
a[1] = _mol->GetAtom(ia+1)->y() - CenterOfSomething[1];
1153
a[2] = _mol->GetAtom(ia+1)->z() - CenterOfSomething[2];
1155
b[0] = _mol->GetAtom(ib+1)->x() - CenterOfSomething[0];
1156
b[1] = _mol->GetAtom(ib+1)->y() - CenterOfSomething[1];
1157
b[2] = _mol->GetAtom(ib+1)->z() - CenterOfSomething[2];
1159
c[0] = _mol->GetAtom(ic+1)->x() - CenterOfSomething[0];
1160
c[1] = _mol->GetAtom(ic+1)->y() - CenterOfSomething[1];
1161
c[2] = _mol->GetAtom(ic+1)->z() - CenterOfSomething[2];
1163
if( ( axis = init_axis_parameters( a, b, c ) ) == NULL ){
1164
if( verbose > 0 ) printf( " no coherrent axis is defined by the points\n" ) ;
1167
axis->transform_atom = rotate_atom ;
1168
if( refine_symmetry_element( axis, 1 ) < 0 ){
1169
if( verbose > 0 ) printf( " refinement failed for the c%d axis\n", axis->order ) ;
1170
destroy_symmetry_element( axis ) ;
1177
* Improper axes-specific routines.
1178
* These are obtained by slight modifications of normal rotation
1182
rotate_reflect_atom( SYMMETRY_ELEMENT *axis, OBAtom *from, OBAtom *to )
1184
double x[3], y[3], a[3], b[3], c[3] ;
1185
double angle = 2*M_PI/axis->order ;
1186
double a_sin = sin( angle ) ;
1187
double a_cos = cos( angle ) ;
1191
if( DIMENSION != 3 ){
1192
fprintf( stderr, "Catastrophe in rotate_reflect_atom!\n" ) ;
1193
exit( EXIT_FAILURE ) ;
1196
x[0] = from->x() - axis->distance * axis->normal[0];
1197
x[1] = from->y() - axis->distance * axis->normal[1];
1198
x[2] = from->z() - axis->distance * axis->normal[2];
1200
for( i = 0, dot = 0 ; i < 3 ; i++ )
1201
dot += x[i] * axis->direction[i] ;
1202
for( i = 0 ; i < 3 ; i++ )
1203
a[i] = axis->direction[i] * dot ;
1204
for( i = 0 ; i < 3 ; i++ )
1205
b[i] = x[i] - a[i] ;
1206
c[0] = b[1]*axis->direction[2] - b[2]*axis->direction[1] ;
1207
c[1] = b[2]*axis->direction[0] - b[0]*axis->direction[2] ;
1208
c[2] = b[0]*axis->direction[1] - b[1]*axis->direction[0] ;
1209
for( i = 0 ; i < 3 ; i++ )
1210
y[i] = -a[i] + b[i]*a_cos + c[i]*a_sin ;
1212
to->SetVector(y[0] + axis->distance * axis->normal[0],
1213
y[1] + axis->distance * axis->normal[1],
1214
y[2] + axis->distance * axis->normal[2]);
1216
// copy the "type" of from into to
1217
to->SetAtomicNum(from->GetAtomicNum());
1218
to->SetIsotope(from->GetIsotope());
1219
to->SetFormalCharge(from->GetFormalCharge());
1220
to->SetSpinMultiplicity(from->GetSpinMultiplicity());
1224
init_improper_axis( int ia, int ib, int ic )
1226
SYMMETRY_ELEMENT * axis ;
1227
double a[ DIMENSION ], b[ DIMENSION ], c[ DIMENSION ] ;
1228
double centerpoint[ DIMENSION ] ;
1232
if( verbose > 0 ) printf( "Trying sn axis for the triplet (%d,%d,%d)\n", ia, ib, ic ) ;
1234
/* First, reduce the problem to Cn case */
1235
a[0] = _mol->GetAtom(ia+1)->x() - CenterOfSomething[0];
1236
a[1] = _mol->GetAtom(ia+1)->y() - CenterOfSomething[1];
1237
a[2] = _mol->GetAtom(ia+1)->z() - CenterOfSomething[2];
1239
b[0] = _mol->GetAtom(ib+1)->x() - CenterOfSomething[0];
1240
b[1] = _mol->GetAtom(ib+1)->y() - CenterOfSomething[1];
1241
b[2] = _mol->GetAtom(ib+1)->z() - CenterOfSomething[2];
1243
c[0] = _mol->GetAtom(ic+1)->x() - CenterOfSomething[0];
1244
c[1] = _mol->GetAtom(ic+1)->y() - CenterOfSomething[1];
1245
c[2] = _mol->GetAtom(ic+1)->z() - CenterOfSomething[2];
1247
for( i = 0, r = 0 ; i < DIMENSION ; i++ ){
1248
centerpoint[i] = a[i] + c[i] + 2*b[i] ;
1249
r += centerpoint[i]*centerpoint[i] ;
1252
if( r <= ToleranceSame ){
1254
if( verbose > 0 ) printf( " atoms can not define improper axis of the order more than 2\n" ) ;
1257
for( i = 0 ; i < DIMENSION ; i++ )
1258
centerpoint[i] /= r ;
1259
for( i = 0, r = 0 ; i < DIMENSION ; i++ )
1260
r += centerpoint[i] * b[i] ;
1261
for( i = 0 ; i < DIMENSION ; i++ )
1262
b[i] = 2*r*centerpoint[i] - b[i] ;
1263
/* Do a quick check of geometry validity */
1264
if( ( axis = init_axis_parameters( a, b, c ) ) == NULL ){
1265
if( verbose > 0 ) printf( " no coherrent improper axis is defined by the points\n" ) ;
1268
axis->transform_atom = rotate_reflect_atom ;
1269
if( refine_symmetry_element( axis, 1 ) < 0 ){
1270
if( verbose > 0 ) printf( " refinement failed for the s%d axis\n", axis->order ) ;
1271
destroy_symmetry_element( axis ) ;
1282
find_center_of_something( void )
1285
double coord_sum[ DIMENSION ] ;
1289
for( j = 0 ; j < DIMENSION ; j++ )
1291
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
1292
atom = _mol->GetAtom(i+1);
1293
coord_sum[0] = atom->x() + atom->y() + atom->z();
1295
for( j = 0 ; j < DIMENSION ; j++ )
1296
CenterOfSomething[j] = coord_sum[j]/_mol->NumAtoms() ;
1298
printf( "Center of something is at %15.10f, %15.10f, %15.10f\n",
1299
CenterOfSomething[0], CenterOfSomething[1], CenterOfSomething[2] ) ;
1300
DistanceFromCenter = (double *) calloc( _mol->NumAtoms(), sizeof( double ) ) ;
1301
if( DistanceFromCenter == NULL ){
1302
fprintf( stderr, "Unable to allocate array for the distances\n" ) ;
1303
exit( EXIT_FAILURE ) ;
1305
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
1306
atom = _mol->GetAtom(i+1);
1307
r = SQUARE(atom->x() - CenterOfSomething[0])
1308
+ SQUARE(atom->y() - CenterOfSomething[1])
1309
+ SQUARE(atom->z() - CenterOfSomething[2]);
1311
DistanceFromCenter[i] = r ;
1319
SYMMETRY_ELEMENT * plane ;
1321
plane = init_ultimate_plane() ;
1322
if( plane != NULL ){
1323
MolecularPlane = plane ;
1325
Planes = (SYMMETRY_ELEMENT **) realloc( Planes, sizeof( SYMMETRY_ELEMENT* ) * PlanesCount ) ;
1326
if( Planes == NULL ){
1327
perror( "Out of memory in find_planes" ) ;
1328
exit( EXIT_FAILURE ) ;
1330
Planes[ PlanesCount - 1 ] = plane ;
1332
for( i = 1 ; i < _mol->NumAtoms() ; i++ ){
1333
for( j = 0 ; j < i ; j++ ){
1334
if( !equivalentAtoms(*_mol->GetAtom(i+1), *_mol->GetAtom(j+1)) )
1336
if( ( plane = init_mirror_plane( i, j ) ) != NULL ){
1338
Planes = (SYMMETRY_ELEMENT **) realloc( Planes, sizeof( SYMMETRY_ELEMENT* ) * PlanesCount ) ;
1339
if( Planes == NULL ){
1340
perror( "Out of memory in find_planes" ) ;
1341
exit( EXIT_FAILURE ) ;
1343
Planes[ PlanesCount - 1 ] = plane ;
1350
find_inversion_centers(void)
1352
SYMMETRY_ELEMENT * center ;
1354
if( ( center = init_inversion_center() ) != NULL ){
1355
InversionCenters = (SYMMETRY_ELEMENT **) calloc( 1, sizeof( SYMMETRY_ELEMENT* ) ) ;
1356
InversionCenters[0] = center ;
1357
InversionCentersCount = 1 ;
1362
find_infinity_axis(void)
1364
SYMMETRY_ELEMENT * axis ;
1366
if( ( axis = init_ultimate_axis() ) != NULL ){
1368
NormalAxes = (SYMMETRY_ELEMENT **) realloc( NormalAxes, sizeof( SYMMETRY_ELEMENT* ) * NormalAxesCount ) ;
1369
if( NormalAxes == NULL ){
1370
perror( "Out of memory in find_infinity_axes()" ) ;
1371
exit( EXIT_FAILURE ) ;
1373
NormalAxes[ NormalAxesCount - 1 ] = axis ;
1381
double center[ DIMENSION ] ;
1382
double * distances = (double*)calloc( _mol->NumAtoms(), sizeof( double ) ) ;
1384
SYMMETRY_ELEMENT * axis ;
1385
OBAtom *a1, *a2, *a3, *a4;
1387
if( distances == NULL ){
1388
fprintf( stderr, "Out of memory in find_c2_axes()\n" ) ;
1389
exit( EXIT_FAILURE ) ;
1391
for( i = 1 ; i < _mol->NumAtoms() ; i++ ){
1392
for( j = 0 ; j < i ; j++ ){
1393
if( !equivalentAtoms(*_mol->GetAtom(i+1), *_mol->GetAtom(j+1)) )
1395
if( fabs( DistanceFromCenter[i] - DistanceFromCenter[j] ) > TolerancePrimary )
1396
continue ; /* A very cheap, but quite effective check */
1398
* First, let's try to get it cheap and use CenterOfSomething
1400
a1 = _mol->GetAtom(i+1);
1401
a2 = _mol->GetAtom(j+1);
1402
center[0] = (a1->x() + a2->x()) / 2.0;
1403
center[1] = (a1->y() + a2->z()) / 2.0;
1404
center[2] = (a1->z() + a2->y()) / 2.0;
1406
r = (vector3(center[0], center[1], center[2])
1407
- vector3(CenterOfSomething[0], CenterOfSomething[1], CenterOfSomething[2])).length();
1409
if( r > 5*TolerancePrimary ){ /* It's Ok to use CenterOfSomething */
1410
if( ( axis = init_c2_axis( i, j, CenterOfSomething ) ) != NULL ){
1412
NormalAxes = (SYMMETRY_ELEMENT **) realloc( NormalAxes, sizeof( SYMMETRY_ELEMENT* ) * NormalAxesCount ) ;
1413
if( NormalAxes == NULL ){
1414
perror( "Out of memory in find_c2_axes" ) ;
1415
exit( EXIT_FAILURE ) ;
1417
NormalAxes[ NormalAxesCount - 1 ] = axis ;
1422
* Now, C2 axis can either pass through an atom, or through the
1423
* middle of the other pair.
1425
for( k = 0 ; k < _mol->NumAtoms() ; k++ ){
1426
if( ( axis = init_c2_axis( i, j, _mol->GetAtom(k+1)->GetVector().AsArray() ) ) != NULL ){
1428
NormalAxes = (SYMMETRY_ELEMENT **) realloc( NormalAxes, sizeof( SYMMETRY_ELEMENT* ) * NormalAxesCount ) ;
1429
if( NormalAxes == NULL ){
1430
perror( "Out of memory in find_c2_axes" ) ;
1431
exit( EXIT_FAILURE ) ;
1433
NormalAxes[ NormalAxesCount - 1 ] = axis ;
1437
* Prepare data for an additional pre-screening check
1439
for( k = 0 ; k < _mol->NumAtoms() ; k++ ){
1440
r = SQUARE(_mol->GetAtom(k+1)->x() - center[0])
1441
+ SQUARE(_mol->GetAtom(k+1)->y() - center[1])
1442
+ SQUARE(_mol->GetAtom(k+1)->z() - center[2]);
1443
distances[k] = sqrt(r) ;
1445
for( k = 0 ; k < _mol->NumAtoms() ; k++ ){
1446
a3 = _mol->GetAtom(k+1);
1447
for( l = 0 ; l < _mol->NumAtoms() ; l++ ){
1448
a4 = _mol->GetAtom(l+1);
1449
if( !equivalentAtoms(*a3, *a4) )
1451
if( fabs( DistanceFromCenter[k] - DistanceFromCenter[l] ) > TolerancePrimary ||
1452
fabs( distances[k] - distances[l] ) > TolerancePrimary )
1453
continue ; /* We really need this one to run reasonably fast! */
1455
center[0] = (a3->x() + a4->x()) / 2.0;
1456
center[1] = (a3->y() + a4->y()) / 2.0;
1457
center[2] = (a3->z() + a4->z()) / 2.0;
1459
if( ( axis = init_c2_axis( i, j, center ) ) != NULL ){
1461
NormalAxes = (SYMMETRY_ELEMENT **) realloc( NormalAxes, sizeof( SYMMETRY_ELEMENT* ) * NormalAxesCount ) ;
1462
if( NormalAxes == NULL ){
1463
perror( "Out of memory in find_c2_axes" ) ;
1464
exit( EXIT_FAILURE ) ;
1466
NormalAxes[ NormalAxesCount - 1 ] = axis ;
1476
find_higher_axes(void)
1479
SYMMETRY_ELEMENT * axis ;
1481
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
1482
for( j = i + 1 ; j < _mol->NumAtoms() ; j++ ){
1483
if( !equivalentAtoms(*_mol->GetAtom(i+1), *_mol->GetAtom(j+1)) )
1485
if( fabs( DistanceFromCenter[i] - DistanceFromCenter[j] ) > TolerancePrimary )
1486
continue ; /* A very cheap, but quite effective check */
1487
for( k = 0 ; k < _mol->NumAtoms() ; k++ ){
1488
if( !equivalentAtoms(*_mol->GetAtom(i+1), *_mol->GetAtom(k+1)) )
1490
if( ( fabs( DistanceFromCenter[i] - DistanceFromCenter[k] ) > TolerancePrimary ) ||
1491
( fabs( DistanceFromCenter[j] - DistanceFromCenter[k] ) > TolerancePrimary ) )
1493
if( ( axis = init_higher_axis( i, j, k ) ) != NULL ){
1495
NormalAxes = (SYMMETRY_ELEMENT **) realloc( NormalAxes, sizeof( SYMMETRY_ELEMENT* ) * NormalAxesCount ) ;
1496
if( NormalAxes == NULL ){
1497
perror( "Out of memory in find_higher_axes" ) ;
1498
exit( EXIT_FAILURE ) ;
1500
NormalAxes[ NormalAxesCount - 1 ] = axis ;
1508
find_improper_axes(void)
1511
SYMMETRY_ELEMENT * axis ;
1513
for( i = 0 ; i < _mol->NumAtoms() ; i++ ){
1514
for( j = i + 1 ; j < _mol->NumAtoms() ; j++ ){
1515
for( k = 0 ; k < _mol->NumAtoms() ; k++ ){
1516
if( ( axis = init_improper_axis( i, j, k ) ) != NULL ){
1517
ImproperAxesCount++ ;
1518
ImproperAxes = (SYMMETRY_ELEMENT **) realloc( ImproperAxes, sizeof( SYMMETRY_ELEMENT* ) * ImproperAxesCount ) ;
1519
if( ImproperAxes == NULL ){
1520
perror( "Out of memory in find_higher_axes" ) ;
1521
exit( EXIT_FAILURE ) ;
1523
ImproperAxes[ ImproperAxesCount - 1 ] = axis ;
1531
report_planes( void )
1535
if( PlanesCount == 0 )
1536
printf( "There are no planes of symmetry in the molecule\n" ) ;
1538
if( PlanesCount == 1 )
1539
printf( "There is a plane of symmetry in the molecule\n" ) ;
1540
else printf( "There are %d planes of symmetry in the molecule\n", PlanesCount ) ;
1541
printf( " Residual Direction of the normal Distance\n" ) ;
1542
for( i = 0 ; i < PlanesCount ; i++ ){
1543
printf( "%3d %8.4e ", i, Planes[i]->maxdev ) ;
1544
printf( "(%11.8f,%11.8f,%11.8f) ", Planes[i]->normal[0], Planes[i]->normal[1], Planes[i]->normal[2] ) ;
1545
printf( "%14.8f\n", Planes[i]->distance ) ;
1551
report_inversion_centers( void )
1553
if( InversionCentersCount == 0 )
1554
printf( "There is no inversion center in the molecule\n" ) ;
1556
printf( "There in an inversion center in the molecule\n" ) ;
1557
printf( " Residual Position\n" ) ;
1558
printf( " %8.4e ", InversionCenters[0]->maxdev ) ;
1559
printf( "(%14.8f,%14.8f,%14.8f)\n",
1560
InversionCenters[0]->distance * InversionCenters[0]->normal[0],
1561
InversionCenters[0]->distance * InversionCenters[0]->normal[1],
1562
InversionCenters[0]->distance * InversionCenters[0]->normal[2] ) ;
1571
if( NormalAxesCount == 0 )
1572
printf( "There are no normal axes in the molecule\n" ) ;
1574
if( NormalAxesCount == 1 )
1575
printf( "There is a normal axis in the molecule\n" ) ;
1576
else printf( "There are %d normal axes in the molecule\n", NormalAxesCount ) ;
1577
printf( " Residual Order Direction of the axis Supporting point\n" ) ;
1578
for( i = 0 ; i < NormalAxesCount ; i++ ){
1579
printf( "%3d %8.4e ", i, NormalAxes[i]->maxdev ) ;
1580
if( NormalAxes[i]->order == 0 )
1582
else printf( "%3d ", NormalAxes[i]->order ) ;
1583
printf( "(%11.8f,%11.8f,%11.8f) ",
1584
NormalAxes[i]->direction[0], NormalAxes[i]->direction[1], NormalAxes[i]->direction[2] ) ;
1585
printf( "(%14.8f,%14.8f,%14.8f)\n",
1586
NormalAxes[0]->distance * NormalAxes[0]->normal[0],
1587
NormalAxes[0]->distance * NormalAxes[0]->normal[1],
1588
NormalAxes[0]->distance * NormalAxes[0]->normal[2] ) ;
1594
report_improper_axes( void )
1598
if( ImproperAxesCount == 0 )
1599
printf( "There are no improper axes in the molecule\n" ) ;
1601
if( ImproperAxesCount == 1 )
1602
printf( "There is an improper axis in the molecule\n" ) ;
1603
else printf( "There are %d improper axes in the molecule\n", ImproperAxesCount ) ;
1604
printf( " Residual Order Direction of the axis Supporting point\n" ) ;
1605
for( i = 0 ; i < ImproperAxesCount ; i++ ){
1606
printf( "%3d %8.4e ", i, ImproperAxes[i]->maxdev ) ;
1607
if( ImproperAxes[i]->order == 0 )
1609
else printf( "%3d ", ImproperAxes[i]->order ) ;
1610
printf( "(%11.8f,%11.8f,%11.8f) ",
1611
ImproperAxes[i]->direction[0], ImproperAxes[i]->direction[1], ImproperAxes[i]->direction[2] ) ;
1612
printf( "(%14.8f,%14.8f,%14.8f)\n",
1613
ImproperAxes[0]->distance * ImproperAxes[0]->normal[0],
1614
ImproperAxes[0]->distance * ImproperAxes[0]->normal[1],
1615
ImproperAxes[0]->distance * ImproperAxes[0]->normal[2] ) ;
1621
* General symmetry handling
1624
report_and_reset_counters( void )
1626
printf( " %10ld candidates examined\n"
1627
" %10ld removed early\n"
1628
" %10ld removed during initial mating stage\n"
1629
" %10ld removed as duplicates\n"
1630
" %10ld removed because of the wrong transformation order\n"
1631
" %10ld removed after unsuccessful optimization\n"
1632
" %10ld accepted\n",
1633
StatTotal, StatEarly, StatPairs, StatDups, StatOrder, StatOpt, StatAccept ) ;
1634
StatTotal = StatEarly = StatPairs = StatDups = StatOrder = StatOpt = StatAccept = 0 ;
1638
find_symmetry_elements( void )
1640
find_center_of_something() ;
1642
printf( "Looking for the inversion center\n" ) ;
1644
find_inversion_centers() ;
1646
report_and_reset_counters() ;
1647
printf( "Looking for the planes of symmetry\n" ) ;
1651
report_and_reset_counters() ;
1652
printf( "Looking for infinity axis\n" ) ;
1654
find_infinity_axis() ;
1656
report_and_reset_counters() ;
1657
printf( "Looking for C2 axes\n" ) ;
1661
report_and_reset_counters() ;
1662
printf( "Looking for higher axes\n" ) ;
1664
find_higher_axes() ;
1666
report_and_reset_counters() ;
1667
printf( "Looking for the improper axes\n" ) ;
1669
find_improper_axes() ;
1671
report_and_reset_counters() ;
1676
compare_axes( const void *a, const void *b )
1678
SYMMETRY_ELEMENT * axis_a = *(SYMMETRY_ELEMENT**) a ;
1679
SYMMETRY_ELEMENT * axis_b = *(SYMMETRY_ELEMENT**) b ;
1680
int i, order_a, order_b ;
1682
order_a = axis_a->order ; if( order_a == 0 ) order_a = 10000 ;
1683
order_b = axis_b->order ; if( order_b == 0 ) order_b = 10000 ;
1684
if( ( i = order_b - order_a ) != 0 ) return i ;
1685
if( axis_a->maxdev > axis_b->maxdev ) return -1 ;
1686
if( axis_a->maxdev < axis_b->maxdev ) return 1 ;
1691
sort_symmetry_elements( void )
1693
if( PlanesCount > 1 ){
1694
qsort( Planes, PlanesCount, sizeof( SYMMETRY_ELEMENT * ), compare_axes ) ;
1696
if( NormalAxesCount > 1 ){
1697
qsort( NormalAxes, NormalAxesCount, sizeof( SYMMETRY_ELEMENT * ), compare_axes ) ;
1699
if( ImproperAxesCount > 1 ){
1700
qsort( ImproperAxes, ImproperAxesCount, sizeof( SYMMETRY_ELEMENT * ), compare_axes ) ;
1705
report_symmetry_elements_verbose( void )
1707
report_inversion_centers() ;
1709
report_improper_axes() ;
1714
summarize_symmetry_elements( void )
1718
NormalAxesCounts = (int*) calloc( MaxAxisOrder+1, sizeof( int ) ) ;
1719
ImproperAxesCounts = (int*) calloc( MaxAxisOrder+1, sizeof( int ) ) ;
1720
for( i = 0 ; i < NormalAxesCount ; i++ )
1721
NormalAxesCounts[ NormalAxes[i]->order ]++ ;
1722
for( i = 0 ; i < ImproperAxesCount ; i++ )
1723
ImproperAxesCounts[ ImproperAxes[i]->order ]++ ;
1727
report_symmetry_elements_brief( void )
1730
char * symmetry_code = (char*)calloc( 1, 10*(PlanesCount+NormalAxesCount+ImproperAxesCount+InversionCentersCount+2) ) ;
1733
if( symmetry_code == NULL ){
1734
fprintf( stderr, "Unable to allocate memory for symmetry ID code in report_symmetry_elements_brief()\n" ) ;
1735
exit( EXIT_FAILURE ) ;
1737
if( PlanesCount + NormalAxesCount + ImproperAxesCount + InversionCentersCount == 0 )
1738
printf( "Molecule has no symmetry elements\n" ) ;
1740
printf( "Molecule has the following symmetry elements: " ) ;
1741
if( InversionCentersCount > 0 ) strcat( symmetry_code, "(i) " ) ;
1742
if( NormalAxesCounts[0] == 1 )
1743
strcat( symmetry_code, "(Cinf) " ) ;
1744
if( NormalAxesCounts[0] > 1 ) {
1745
snprintf( buf, 100, "%d*(Cinf) ", NormalAxesCounts[0] ) ;
1746
strcat( symmetry_code, buf ) ;
1748
for( i = MaxAxisOrder ; i >= 2 ; i-- ){
1749
if( NormalAxesCounts[i] == 1 ){ snprintf( buf, 100, "(C%d) ", i ) ; strcat( symmetry_code, buf ) ; }
1750
if( NormalAxesCounts[i] > 1 ){ snprintf( buf, 100, "%d*(C%d) ", NormalAxesCounts[i], i ) ; strcat( symmetry_code, buf ) ; }
1752
for( i = MaxAxisOrder ; i >= 2 ; i-- ){
1753
if( ImproperAxesCounts[i] == 1 ){ snprintf( buf, 100, "(S%d) ", i ) ; strcat( symmetry_code, buf ) ; }
1754
if( ImproperAxesCounts[i] > 1 ){ snprintf( buf, 100, "%d*(S%d) ", ImproperAxesCounts[i], i ) ; strcat( symmetry_code, buf ) ; }
1756
if( PlanesCount == 1 ) strcat( symmetry_code, "(sigma) " ) ;
1757
if( PlanesCount > 1 ){ snprintf( buf, 100, "%d*(sigma) ", PlanesCount ) ; strcat( symmetry_code, buf ) ; }
1758
printf( "%s\n", symmetry_code ) ;
1760
SymmetryCode = symmetry_code ;
1763
const char *identify_point_group( void )
1766
int last_matching = -1;
1767
int matching_count = 0;
1769
for( i = 0 ; i < PointGroupsCount ; i++ ){
1770
if( strcmp( SymmetryCode, PointGroups[i].symmetry_code ) == 0 ){
1775
if( matching_count == 0 ){
1776
printf( "These symmetry elements match no point group I know of. Sorry.\n" ) ;
1778
if( matching_count > 1 ){
1779
printf( "These symmetry elements match more than one group I know of.\n"
1780
"SOMETHING IS VERY WRONG\n" ) ;
1781
printf( "Matching groups are:\n" ) ;
1782
for( i = 0 ; i < PointGroupsCount ; i++ ){
1783
if( ( strcmp( SymmetryCode, PointGroups[i].symmetry_code ) == 0 )) {
1784
printf( " %s\n", PointGroups[i].group_name ) ;
1788
if( matching_count == 1 ){
1789
printf( "It seems to be the %s point group\n", PointGroups[last_matching].group_name ) ;
1791
return PointGroups[last_matching].group_name;
1794
}; // end class PointGroupPrivate
1796
OBPointGroup::OBPointGroup()
1798
d = new PointGroupPrivate;
1801
OBPointGroup::~OBPointGroup()
1806
void OBPointGroup::Setup(OBMol *mol)
1811
const char* OBPointGroup::IdentifyPointGroup()
1813
d->find_symmetry_elements();
1814
d->sort_symmetry_elements();
1815
d->summarize_symmetry_elements();
1816
if( d->BadOptimization ) {
1820
if( d->verbose >= 0 )
1821
d->report_symmetry_elements_verbose();
1822
d->report_symmetry_elements_brief();
1823
return d->identify_point_group();
1826
} // end namespace OpenBabel
1828
//! \file pointgroup.cpp
1829
//! \brief Brute-force point group detection
added
removed
src/pointgroup.cpp
