Octopus
symmetries_finite.c
Go to the documentation of this file.
1/*
2 * Brute force symmetry analyzer.
3 * This is actually C++ program, masquerading as a C one!
4 *
5 * (C) 2011 X. Andrade Converted to a Fortran library.
6 * (C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
21 *
22 * Revision 1.16 2003/04/04 13:05:03 patchkov
23 * Revision 1.15 2000/01/25 16:47:17 patchkov
24 * Revision 1.14 2000/01/25 16:39:08 patchkov
25 * Revision 1.13 1996/05/24 12:32:08 ps
26 * Revision 1.12 1996/05/23 16:10:47 ps
27 * First reasonably stable version.
28 *
29 */
30
31#include <config.h>
32
33#include <math.h>
34#include <stdio.h>
35#include <stdlib.h>
36#include <string.h>
37
38#define DIMENSION 3
39#define MAXPARAM 7
40
41typedef struct {
42 int type;
43 double x[DIMENSION];
44} ATOM;
45
46/*
47 * All specific structures should have corresponding elements in the
48 * same position generic structure does.
49 *
50 * Planes are characterized by the surface normal direction
51 * (taken in the direction *from* the coordinate origin)
52 * and distance from the coordinate origin to the plane
53 * in the direction of the surface normal.
54 *
55 * Inversion is characterized by location of the inversion center.
56 *
57 * Rotation is characterized by a vector (distance+direction) from the origin
58 * to the rotation axis, axis direction and rotation order. Rotations
59 * are in the clockwise direction looking opposite to the direction
60 * of the axis. Note that this definition of the rotation axis
61 * is *not* unique, since an arbitrary multiple of the axis direction
62 * can be added to the position vector without changing actual operation.
63 *
64 * Mirror rotation is defined by the same parameters as normal rotation,
65 * but the origin is now unambiguous since it defines the position of the
66 * plane associated with the axis.
67 *
68 */
69
70typedef struct _SYMMETRY_ELEMENT_ {
71 void (*transform_atom)(struct _SYMMETRY_ELEMENT_ *el, ATOM *from, ATOM *to);
72 int *transform; /* Correspondence table for the transformation */
73 int order; /* Applying transformation this many times is identity */
74 int nparam; /* 4 for inversion and planes, 7 for axes */
75 double maxdev; /* Larges error associated with the element */
76 double distance;
77 double normal[DIMENSION];
78 double direction[DIMENSION];
80
81typedef struct {
82 char *group_name; /* Canonical group name */
83 char *symmetry_code; /* Group symmetry code */
84 int (*check)(void); /* Additional verification routine, not used */
87static double ToleranceSame = 1e-3;
88static double TolerancePrimary = 5e-2;
89static double ToleranceFinal = 1e-4;
90static double MaxOptStep = 5e-1;
91static double MinOptStep = 1e-7;
92static double GradientStep = 1e-7;
93static double OptChangeThreshold = 1e-10;
94static double CenterOfSomething[DIMENSION];
95static double *DistanceFromCenter = NULL;
96static int verbose = 0;
97static int MaxOptCycles = 200;
98static int OptChangeHits = 5;
99static int MaxAxisOrder = 20;
100static int AtomsCount = 0;
101static ATOM *Atoms = NULL;
102static int PlanesCount = 0;
103static SYMMETRY_ELEMENT **Planes = NULL;
104static SYMMETRY_ELEMENT *MolecularPlane = NULL;
105static int InversionCentersCount = 0;
107static int NormalAxesCount = 0;
108static SYMMETRY_ELEMENT **NormalAxes = NULL;
109static int ImproperAxesCount = 0;
110static SYMMETRY_ELEMENT **ImproperAxes = NULL;
111static int *NormalAxesCounts = NULL;
112static int *ImproperAxesCounts = NULL;
113static int BadOptimization = 0;
114static char *SymmetryCode = NULL; /*"" ;*/
116 * Statistics
117 */
118static long StatTotal = 0;
119static long StatEarly = 0;
120static long StatPairs = 0;
121static long StatDups = 0;
122static long StatOrder = 0;
123static long StatOpt = 0;
124static long StatAccept = 0;
125
127 * Point groups I know about
128 */
129int true(void) { return 1; }
131 {"C1", "", true},
132 {"Cs", "(sigma) ", true},
133 {"Ci", "(i) ", true},
134 {"C2", "(C2) ", true},
135 {"C3", "(C3) ", true},
136 {"C4", "(C4) (C2) ", true},
137 {"C5", "(C5) ", true},
138 {"C6", "(C6) (C3) (C2) ", true},
139 {"C7", "(C7) ", true},
140 {"C8", "(C8) (C4) (C2) ", true},
141 {"D2", "3*(C2) ", true},
142 {"D3", "(C3) 3*(C2) ", true},
143 {"D4", "(C4) 5*(C2) ", true},
144 {"D5", "(C5) 5*(C2) ", true},
145 {"D6", "(C6) (C3) 7*(C2) ", true},
146 {"D7", "(C7) 7*(C2) ", true},
147 {"D8", "(C8) (C4) 9*(C2) ", true},
148 {"C2v", "(C2) 2*(sigma) ", true},
149 {"C3v", "(C3) 3*(sigma) ", true},
150 {"C4v", "(C4) (C2) 4*(sigma) ", true},
151 {"C5v", "(C5) 5*(sigma) ", true},
152 {"C6v", "(C6) (C3) (C2) 6*(sigma) ", true},
153 {"C7v", "(C7) 7*(sigma) ", true},
154 {"C8v", "(C8) (C4) (C2) 8*(sigma) ", true},
155 {"C2h", "(i) (C2) (sigma) ", true},
156 {"C3h", "(C3) (S3) (sigma) ", true},
157 {"C4h", "(i) (C4) (C2) (S4) (sigma) ", true},
158 {"C5h", "(C5) (S5) (sigma) ", true},
159 {"C6h", "(i) (C6) (C3) (C2) (S6) (S3) (sigma) ", true},
160 {"C7h", "(C7) (S7) (sigma) ", true},
161 {"C8h", "(i) (C8) (C4) (C2) (S8) (S4) (sigma) ", true},
162 {"D2h", "(i) 3*(C2) 3*(sigma) ", true},
163 {"D3h", "(C3) 3*(C2) (S3) 4*(sigma) ", true},
164 {"D4h", "(i) (C4) 5*(C2) (S4) 5*(sigma) ", true},
165 {"D5h", "(C5) 5*(C2) (S5) 6*(sigma) ", true},
166 {"D6h", "(i) (C6) (C3) 7*(C2) (S6) (S3) 7*(sigma) ", true},
167 {"D7h", "(C7) 7*(C2) (S7) 8*(sigma) ", true},
168 {"D8h", "(i) (C8) (C4) 9*(C2) (S8) (S4) 9*(sigma) ", true},
169 {"D2d", "3*(C2) (S4) 2*(sigma) ", true},
170 {"D3d", "(i) (C3) 3*(C2) (S6) 3*(sigma) ", true},
171 {"D4d", "(C4) 5*(C2) (S8) 4*(sigma) ", true},
172 {"D5d", "(i) (C5) 5*(C2) (S10) 5*(sigma) ", true},
173 {"D6d", "(C6) (C3) 7*(C2) (S12) (S4) 6*(sigma) ", true},
174 {"D7d", "(i) (C7) 7*(C2) (S14) 7*(sigma) ", true},
175 {"D8d", "(C8) (C4) 9*(C2) (S16) 8*(sigma) ", true},
176 {"S4", "(C2) (S4) ", true},
177 {"S6", "(i) (C3) (S6) ", true},
178 {"S8", "(C4) (C2) (S8) ", true},
179 {"T", "4*(C3) 3*(C2) ", true},
180 {"Th", "(i) 4*(C3) 3*(C2) 4*(S6) 3*(sigma) ", true},
181 {"Td", "4*(C3) 3*(C2) 3*(S4) 6*(sigma) ", true},
182 {"O", "3*(C4) 4*(C3) 9*(C2) ", true},
183 {"Oh", "(i) 3*(C4) 4*(C3) 9*(C2) 4*(S6) 3*(S4) 9*(sigma) ", true},
184 {"Cinfv", "(Cinf) (sigma) ", true},
185 {"Dinfh", "(i) (Cinf) (C2) 2*(sigma) ", true},
186 {"I", "6*(C5) 10*(C3) 15*(C2) ", true},
187 {"Ih", "(i) 6*(C5) 10*(C3) 15*(C2) 6*(S10) 10*(S6) 15*(sigma) ", true},
188 {"Kh", "(i) (Cinf) (sigma) ", true},
189};
190#define PointGroupsCount (sizeof(PointGroups) / sizeof(POINT_GROUP))
191char *PointGroupRejectionReason = NULL;
192
193/*
194 * Generic functions
195 */
196
197double pow2(double x) { return x * x; }
198
200 int i, j, k, best_j;
201 char *atom_used = calloc(AtomsCount, 1);
202 double distance, best_distance;
203 ATOM symmetric;
204
205 if (atom_used == NULL) {
206 fprintf(stderr, "Out of memory for tagging array in establish_pairs()\n");
207 exit(EXIT_FAILURE);
209 for (i = 0; i < AtomsCount; i++) {
210 if (elem->transform[i] >= AtomsCount) { /* No symmetric atom yet */
211 if (verbose > 2)
212 printf(" looking for a pair for %d\n", i);
213 elem->transform_atom(elem, Atoms + i, &symmetric);
214 if (verbose > 2)
215 printf(" new coordinates are: (%g,%g,%g)\n", symmetric.x[0],
216 symmetric.x[1], symmetric.x[2]);
217 best_j = i;
218 best_distance = 2 * TolerancePrimary; /* Performance value we'll reject */
219 for (j = 0; j < AtomsCount; j++) {
220 if (Atoms[j].type != symmetric.type || atom_used[j])
221 continue;
222 for (k = 0, distance = 0; k < DIMENSION; k++) {
223 distance += pow2(symmetric.x[k] - Atoms[j].x[k]);
224 }
226 if (verbose > 2)
227 printf(" distance to %d is %g\n", j, distance);
228 if (distance < best_distance) {
229 best_j = j;
230 best_distance = distance;
231 }
233 if (best_distance >
234 TolerancePrimary) { /* Too bad, there is no symmetric atom */
235 if (verbose > 0)
236 printf(" no pair for atom %d - best was %d with err = %g\n", i,
237 best_j, best_distance);
238 free(atom_used);
239 return -1;
241 elem->transform[i] = best_j;
242 atom_used[best_j] = 1;
243 if (verbose > 1)
244 printf(" atom %d transforms to the atom %d, err = %g\n", i,
245 best_j, best_distance);
246 }
247 }
248 free(atom_used);
249 return 0;
251
253 int i, j, k;
255
256 for (i = 0; i < AtomsCount; i++) {
257 if (elem->transform[i] == i) /* Identity transform is Ok for any order */
258 continue;
260 j = elem->transform[i];
261 if (elem->transform[j] == i)
262 continue; /* Second-order transform is Ok for improper axis */
263 }
264 for (j = elem->order - 1, k = elem->transform[i]; j > 0;
265 j--, k = elem->transform[k]) {
266 if (k == i) {
267 if (verbose > 0)
268 printf(" transform looped %d steps too early from atom %d\n",
269 j, i);
270 return -1;
271 }
273 if (k != i && elem->transform_atom == rotate_reflect_atom) {
274 /* For improper axes, the complete loop may also take twice the order */
275 for (j = elem->order; j > 0; j--, k = elem->transform[k]) {
276 if (k == i) {
277 if (verbose > 0)
278 printf(" (improper) transform looped %d steps too early "
279 "from atom %d\n",
280 j, i);
281 return -1;
283 }
284 }
285 if (k != i) {
286 if (verbose > 0)
287 printf(" transform failed to loop after %d steps from atom %d\n",
288 elem->order, i);
289 return -1;
290 }
292 return 0;
293}
296 int i, j;
297 int code;
298
299 if ((a->order != b->order) || (a->nparam != b->nparam) ||
301 return 0;
302 for (i = 0, code = 1; i < AtomsCount; i++) {
303 if (a->transform[i] != b->transform[i]) {
304 code = 0;
305 break;
306 }
307 }
308 if (code == 0 &&
309 a->order > 2) { /* b can also be a reverse transformation for a */
310 for (i = 0; i < AtomsCount; i++) {
311 j = a->transform[i];
312 if (b->transform[j] != i)
313 return 0;
314 }
315 return 1;
316 }
317 return code;
318}
319
321 SYMMETRY_ELEMENT *elem = calloc(1, sizeof(SYMMETRY_ELEMENT));
322 int i;
324 if (elem == NULL) {
325 fprintf(stderr, "Out of memory allocating symmetry element\n");
326 exit(EXIT_FAILURE);
328 elem->transform = calloc(AtomsCount, sizeof(int));
329 if (elem->transform == NULL) {
330 fprintf(stderr,
331 "Out of memory allocating transform table for symmetry element\n");
332 exit(EXIT_FAILURE);
333 }
334 for (i = 0; i < AtomsCount; i++) {
335 elem->transform[i] = AtomsCount + 1; /* An impossible value */
336 }
337 return elem;
338}
341 if (elem != NULL) {
342 if (elem->transform != NULL)
343 free(elem->transform);
344 free(elem);
346}
347
349 int i, j, k;
350 ATOM symmetric;
351 double r, max_r;
352
353 for (i = 0, max_r = 0; i < AtomsCount; i++) {
354 j = elem->transform[i];
355 elem->transform_atom(elem, Atoms + i, &symmetric);
356 for (k = 0, r = 0; k < DIMENSION; k++) {
357 r += pow2(symmetric.x[k] - Atoms[j].x[k]);
358 }
359 r = sqrt(r);
360 if (r > ToleranceFinal) {
361 if (verbose > 0)
362 printf(" distance to symmetric atom (%g) is too big for %d\n", r,
363 i);
364 return -1;
365 }
366 if (r > max_r)
367 max_r = r;
368 }
369 elem->maxdev = max_r;
370 return 0;
372
373double eval_optimization_target_function(SYMMETRY_ELEMENT *elem, int *finish) {
374 int i, j, k;
375 ATOM symmetric;
376 double target, r, maxr;
377
378 if (elem->nparam >= 4) {
379 for (k = 0, r = 0; k < DIMENSION; k++) {
380 r += elem->normal[k] * elem->normal[k];
381 }
382 r = sqrt(r);
383 if (r < ToleranceSame) {
384 fprintf(stderr, "Normal collapsed!\n");
385 exit(EXIT_FAILURE);
386 }
387 for (k = 0; k < DIMENSION; k++) {
388 elem->normal[k] /= r;
389 }
390 if (elem->distance < 0) {
391 elem->distance = -elem->distance;
392 for (k = 0; k < DIMENSION; k++) {
393 elem->normal[k] = -elem->normal[k];
394 }
395 }
396 }
397 if (elem->nparam >= 7) {
398 for (k = 0, r = 0; k < DIMENSION; k++) {
399 r += elem->direction[k] * elem->direction[k];
400 }
401 r = sqrt(r);
402 if (r < ToleranceSame) {
403 fprintf(stderr, "Direction collapsed!\n");
404 exit(EXIT_FAILURE);
405 }
406 for (k = 0; k < DIMENSION; k++) {
407 elem->direction[k] /= r;
408 }
409 }
410 for (i = 0, target = maxr = 0; i < AtomsCount; i++) {
411 elem->transform_atom(elem, Atoms + i, &symmetric);
412 j = elem->transform[i];
413 for (k = 0, r = 0; k < DIMENSION; k++) {
414 r += pow2(Atoms[j].x[k] - symmetric.x[k]);
415 }
416 if (r > maxr)
417 maxr = r;
418 target += r;
419 }
420 if (finish != NULL) {
421 *finish = 0;
422 if (sqrt(maxr) < ToleranceFinal)
423 *finish = 1;
424 }
425 return target;
426}
427
428void get_params(SYMMETRY_ELEMENT *elem, double values[]) {
429 memcpy(values, &elem->distance, elem->nparam * sizeof(double));
430}
431
432void set_params(SYMMETRY_ELEMENT *elem, double values[]) {
433 memcpy(&elem->distance, values, elem->nparam * sizeof(double));
434}
435
437 double values[MAXPARAM];
438 double grad[MAXPARAM];
439 double force[MAXPARAM];
440 double step[MAXPARAM];
441 double f, fold, fnew, fnew2, fdn, fup, snorm;
442 double a, b, x;
443 int vars = elem->nparam;
444 int cycle = 0;
445 int i, finish;
446 int hits = 0;
447
448 if (vars > MAXPARAM) {
449 fprintf(stderr, "Catastrophe in optimize_transformation_params()!\n");
450 exit(EXIT_FAILURE);
451 }
452 f = 0;
453 do {
454 fold = f;
455 f = eval_optimization_target_function(elem, &finish);
456 /* Evaluate function, gradient and diagonal force constants */
457 if (verbose > 1)
458 printf(" function value = %g\n", f);
459 if (finish) {
460 if (verbose > 1)
461 printf(" function value is small enough\n");
462 break;
463 }
464 if (cycle > 0) {
465 if (fabs(f - fold) > OptChangeThreshold)
466 hits = 0;
467 else
468 hits++;
469 if (hits >= OptChangeHits) {
470 if (verbose > 1)
471 printf(" no progress is made, stop optimization\n");
472 break;
473 }
474 }
475 get_params(elem, values);
476 for (i = 0; i < vars; i++) {
478 set_params(elem, values);
479 fdn = eval_optimization_target_function(elem, NULL);
480 values[i] += 2 * GradientStep;
481 set_params(elem, values);
482 fup = eval_optimization_target_function(elem, NULL);
484 grad[i] = (fup - fdn) / (2 * GradientStep);
485 force[i] = (fup + fdn - 2 * f) / (GradientStep * GradientStep);
486 if (verbose > 1)
487 printf(" i = %d, grad = %12.6e, force = %12.6e\n", i, grad[i],
488 force[i]);
489 }
490 /* Do a quasy-Newton step */
491 for (i = 0, snorm = 0; i < vars; i++) {
492 if (force[i] < 0)
493 force[i] = -force[i];
494 if (force[i] < 1e-3)
495 force[i] = 1e-3;
496 if (force[i] > 1e3)
497 force[i] = 1e3;
498 step[i] = -grad[i] / force[i];
499 snorm += step[i] * step[i];
500 }
501 snorm = sqrt(snorm);
502 if (snorm > MaxOptStep) { /* Renormalize step */
503 for (i = 0; i < vars; i++)
504 step[i] *= MaxOptStep / snorm;
505 snorm = MaxOptStep;
506 }
507 do {
508 for (i = 0; i < vars; i++) {
509 values[i] += step[i];
510 }
511 set_params(elem, values);
512 fnew = eval_optimization_target_function(elem, NULL);
513 if (fnew < f)
514 break;
515 for (i = 0; i < vars; i++) {
516 values[i] -= step[i];
517 step[i] /= 2;
518 }
519 set_params(elem, values);
520 snorm /= 2;
521 } while (snorm > MinOptStep);
522 if ((snorm > MinOptStep) &&
523 (snorm < MaxOptStep / 2)) { /* try to do quadratic interpolation */
524 for (i = 0; i < vars; i++)
525 values[i] += step[i];
526 set_params(elem, values);
527 fnew2 = eval_optimization_target_function(elem, NULL);
528 if (verbose > 1)
529 printf(" interpolation base points: %g, %g, %g\n", f, fnew,
530 fnew2);
531 for (i = 0; i < vars; i++)
532 values[i] -= 2 * step[i];
533 a = (4 * f - fnew2 - 3 * fnew) / 2;
534 b = (f + fnew2 - 2 * fnew) / 2;
535 if (verbose > 1)
536 printf(" linear interpolation coefficients %g, %g\n", a, b);
537 if (b > 0) {
538 x = -a / (2 * b);
539 if (x > 0.2 && x < 1.8) {
540 if (verbose > 1)
541 printf(" interpolated: %g\n", x);
542 for (i = 0; i < vars; i++)
543 values[i] += x * step[i];
544 } else
545 b = 0;
546 }
547 if (b <= 0) {
548 if (fnew2 < fnew) {
549 for (i = 0; i < vars; i++)
550 values[i] += 2 * step[i];
551 } else {
552 for (i = 0; i < vars; i++)
553 values[i] += step[i];
554 }
555 }
556 set_params(elem, values);
557 }
558 } while (snorm > MinOptStep && ++cycle < MaxOptCycles);
560 if (cycle >= MaxOptCycles)
561 BadOptimization = 1;
562 if (verbose > 0) {
563 if (cycle >= MaxOptCycles)
564 printf(" maximum number of optimization cycles made\n");
565 printf(" optimization completed after %d cycles with f = %g\n",
566 cycle, f);
567 }
568}
569
570int refine_symmetry_element(SYMMETRY_ELEMENT *elem, int build_table) {
571 int i;
572
573 if (build_table && (establish_pairs(elem) < 0)) {
574 StatPairs++;
575 if (verbose > 0)
576 printf(" no transformation correspondence table can be "
577 "constructed\n");
578 return -1;
579 }
580 for (i = 0; i < PlanesCount; i++) {
581 if (same_transform(Planes[i], elem)) {
582 StatDups++;
583 if (verbose > 0)
584 printf(" transformation is identical to plane %d\n", i);
585 return -1;
586 }
587 }
588 for (i = 0; i < InversionCentersCount; i++) {
589 if (same_transform(InversionCenters[i], elem)) {
590 StatDups++;
591 if (verbose > 0)
592 printf(" transformation is identical to inversion center %d\n",
593 i);
594 return -1;
595 }
596 }
597 for (i = 0; i < NormalAxesCount; i++) {
598 if (same_transform(NormalAxes[i], elem)) {
599 StatDups++;
600 if (verbose > 0)
601 printf(" transformation is identical to normal axis %d\n", i);
602 return -1;
603 }
604 }
605 for (i = 0; i < ImproperAxesCount; i++) {
606 if (same_transform(ImproperAxes[i], elem)) {
607 StatDups++;
608 if (verbose > 0)
609 printf(" transformation is identical to improper axis %d\n", i);
610 return -1;
611 }
612 }
613 if (check_transform_order(elem) < 0) {
614 StatOrder++;
615 if (verbose > 0)
616 printf(" incorrect transformation order\n");
617 return -1;
618 }
620 if (check_transform_quality(elem) < 0) {
621 StatOpt++;
622 if (verbose > 0)
623 printf(" refined transformation does not pass the numeric "
624 "threshold\n");
625 return -1;
626 }
627 StatAccept++;
628 return 0;
629}
630
631/*
632 * Plane-specific functions
633 */
634
635void mirror_atom(SYMMETRY_ELEMENT *plane, ATOM *from, ATOM *to) {
636 int i;
637 double r;
638
639 for (i = 0, r = plane->distance; i < DIMENSION; i++) {
640 r -= from->x[i] * plane->normal[i];
641 }
642 to->type = from->type;
643 for (i = 0; i < DIMENSION; i++) {
644 to->x[i] = from->x[i] + 2 * r * plane->normal[i];
645 }
646}
647
650 double dx[DIMENSION], midpoint[DIMENSION], rab, r;
651 int k;
652
653 if (verbose > 0)
654 printf("Trying mirror plane for atoms %d,%d\n", i, j);
655 StatTotal++;
657 plane->order = 2;
658 plane->nparam = 4;
659 for (k = 0, rab = 0; k < DIMENSION; k++) {
660 dx[k] = Atoms[i].x[k] - Atoms[j].x[k];
661 midpoint[k] = (Atoms[i].x[k] + Atoms[j].x[k]) / 2.0;
662 rab += dx[k] * dx[k];
663 }
664 rab = sqrt(rab);
665 if (rab < ToleranceSame) {
666 fprintf(stderr, "Atoms %d and %d coincide (r = %g)\n", i, j, rab);
667 exit(EXIT_FAILURE);
668 }
669 for (k = 0, r = 0; k < DIMENSION; k++) {
670 plane->normal[k] = dx[k] / rab;
671 r += midpoint[k] * plane->normal[k];
672 }
673 if (r < 0) { /* Reverce normal direction, distance is always positive! */
674 r = -r;
675 for (k = 0; k < DIMENSION; k++) {
676 plane->normal[k] = -plane->normal[k];
677 }
678 }
679 plane->distance = r;
680 if (verbose > 0)
681 printf(" initial plane is at %g from the origin\n", r);
682 if (refine_symmetry_element(plane, 1) < 0) {
683 if (verbose > 0)
684 printf(" refinement failed for the plane\n");
686 return NULL;
687 }
688 return plane;
689}
690
693 double d0[DIMENSION], d1[DIMENSION], d2[DIMENSION];
694 double p[DIMENSION];
695 double r, s0, s1, s2;
696 double *d;
697 int i, j, k;
698
699 if (verbose > 0)
700 printf("Trying whole-molecule mirror plane\n");
701 StatTotal++;
703 plane->order = 1;
704 plane->nparam = 4;
705 for (k = 0; k < DIMENSION; k++)
706 d0[k] = d1[k] = d2[k] = 0;
707 d0[0] = 1;
708 d1[1] = 1;
709 d2[2] = 1;
710 for (i = 1; i < AtomsCount; i++) {
711 for (j = 0; j < i; j++) {
712 for (k = 0, r = 0; k < DIMENSION; k++) {
713 p[k] = Atoms[i].x[k] - Atoms[j].x[k];
714 r += p[k] * p[k];
715 }
716 r = sqrt(r);
717 for (k = 0, s0 = s1 = s2 = 0; k < DIMENSION; k++) {
718 p[k] /= r;
719 s0 += p[k] * d0[k];
720 s1 += p[k] * d1[k];
721 s2 += p[k] * d2[k];
722 }
723 for (k = 0; k < DIMENSION; k++) {
724 d0[k] -= s0 * p[k];
725 d1[k] -= s1 * p[k];
726 d2[k] -= s2 * p[k];
727 }
728 }
729 }
730 for (k = 0, s0 = s1 = s2 = 0; k < DIMENSION; k++) {
731 s0 += d0[k];
732 s1 += d1[k];
733 s2 += d2[k];
734 }
735 d = NULL;
736 if (s0 >= s1 && s0 >= s2)
737 d = d0;
738 if (s1 >= s0 && s1 >= s2)
739 d = d1;
740 if (s2 >= s0 && s2 >= s1)
741 d = d2;
742 if (d == NULL) {
743 fprintf(stderr,
744 "Catastrophe in init_ultimate_plane(): %g, %g and %g have no "
745 "ordering!\n",
746 s0, s1, s2);
747 exit(EXIT_FAILURE);
748 }
749 for (k = 0, r = 0; k < DIMENSION; k++)
750 r += d[k] * d[k];
751 r = sqrt(r);
752 if (r > 0) {
753 for (k = 0; k < DIMENSION; k++)
754 plane->normal[k] = d[k] / r;
755 } else {
756 for (k = 1; k < DIMENSION; k++)
757 plane->normal[k] = 0;
758 plane->normal[0] = 1;
759 }
760 for (k = 0, r = 0; k < DIMENSION; k++)
761 r += CenterOfSomething[k] * plane->normal[k];
762 plane->distance = r;
763 for (k = 0; k < AtomsCount; k++)
764 plane->transform[k] = k;
765 if (refine_symmetry_element(plane, 0) < 0) {
766 if (verbose > 0)
767 printf(" refinement failed for the plane\n");
769 return NULL;
770 }
771 return plane;
772}
773/*
774 * Inversion-center specific functions
775 */
776void invert_atom(SYMMETRY_ELEMENT *center, ATOM *from, ATOM *to) {
777 int i;
778
779 to->type = from->type;
780 for (i = 0; i < DIMENSION; i++) {
781 to->x[i] = 2 * center->distance * center->normal[i] - from->x[i];
782 }
783}
784
787 int k;
788 double r;
789
790 if (verbose > 0)
791 printf("Trying inversion center at the center of something\n");
792 StatTotal++;
793 center->transform_atom = invert_atom;
794 center->order = 2;
795 center->nparam = 4;
796 for (k = 0, r = 0; k < DIMENSION; k++)
798 r = sqrt(r);
799 if (r > 0) {
800 for (k = 0; k < DIMENSION; k++)
801 center->normal[k] = CenterOfSomething[k] / r;
802 } else {
803 center->normal[0] = 1;
804 for (k = 1; k < DIMENSION; k++)
805 center->normal[k] = 0;
806 }
807 center->distance = r;
808 if (verbose > 0)
809 printf(" initial inversion center is at %g from the origin\n", r);
810 if (refine_symmetry_element(center, 1) < 0) {
811 if (verbose > 0)
812 printf(" refinement failed for the inversion center\n");
814 return NULL;
815 }
816 return center;
817}
818
819/*
820 * Normal rotation axis-specific routines.
821 */
822void rotate_atom(SYMMETRY_ELEMENT *axis, ATOM *from, ATOM *to) {
823 double x[3], y[3], a[3], b[3], c[3];
824 double angle = axis->order ? 2 * M_PI / axis->order : 1.0;
825 double a_sin = sin(angle);
826 double a_cos = cos(angle);
827 double dot;
828 int i;
829
830 if (DIMENSION != 3) {
831 fprintf(stderr, "Catastrophe in rotate_atom!\n");
832 exit(EXIT_FAILURE);
833 }
834 for (i = 0; i < 3; i++)
835 x[i] = from->x[i] - axis->distance * axis->normal[i];
836 for (i = 0, dot = 0; i < 3; i++)
837 dot += x[i] * axis->direction[i];
838 for (i = 0; i < 3; i++)
839 a[i] = axis->direction[i] * dot;
840 for (i = 0; i < 3; i++)
841 b[i] = x[i] - a[i];
842 c[0] = b[1] * axis->direction[2] - b[2] * axis->direction[1];
843 c[1] = b[2] * axis->direction[0] - b[0] * axis->direction[2];
844 c[2] = b[0] * axis->direction[1] - b[1] * axis->direction[0];
845 for (i = 0; i < 3; i++)
846 y[i] = a[i] + b[i] * a_cos + c[i] * a_sin;
847 for (i = 0; i < 3; i++)
848 to->x[i] = y[i] + axis->distance * axis->normal[i];
849 to->type = from->type;
850}
851
854 double dir[DIMENSION], rel[DIMENSION];
855 double s;
856 int i, k;
857
858 if (verbose > 0)
859 printf("Trying infinity axis\n");
860 StatTotal++;
862 axis->order = 0;
863 axis->nparam = 7;
864 for (k = 0; k < DIMENSION; k++)
865 dir[k] = 0;
866 for (i = 0; i < AtomsCount; i++) {
867 for (k = 0, s = 0; k < DIMENSION; k++) {
868 rel[k] = Atoms[i].x[k] - CenterOfSomething[k];
869 s += rel[k] * dir[k];
870 }
871 if (s >= 0)
872 for (k = 0; k < DIMENSION; k++)
873 dir[k] += rel[k];
874 else
875 for (k = 0; k < DIMENSION; k++)
876 dir[k] -= rel[k];
877 }
878 for (k = 0, s = 0; k < DIMENSION; k++)
879 s += pow2(dir[k]);
880 s = sqrt(s);
881 if (s > 0)
882 for (k = 0; k < DIMENSION; k++)
883 dir[k] /= s;
884 else
885 dir[0] = 1;
886 for (k = 0; k < DIMENSION; k++)
887 axis->direction[k] = dir[k];
888 for (k = 0, s = 0; k < DIMENSION; k++)
889 s += pow2(CenterOfSomething[k]);
890 s = sqrt(s);
891 if (s > 0)
892 for (k = 0; k < DIMENSION; k++)
893 axis->normal[k] = CenterOfSomething[k] / s;
894 else {
895 for (k = 1; k < DIMENSION; k++)
896 axis->normal[k] = 0;
897 axis->normal[0] = 1;
898 }
899 axis->distance = s;
900 for (k = 0; k < AtomsCount; k++)
901 axis->transform[k] = k;
902 if (refine_symmetry_element(axis, 0) < 0) {
903 if (verbose > 0)
904 printf(" refinement failed for the infinity axis\n");
906 return NULL;
907 }
908 return axis;
909}
910
911SYMMETRY_ELEMENT *init_c2_axis(int i, int j, double support[DIMENSION]) {
912 SYMMETRY_ELEMENT *axis;
913 int k;
914 double ris, rjs;
915 double r, center[DIMENSION];
916
917 if (verbose > 0)
918 printf("Trying c2 axis for the pair (%d,%d) with the support (%g,%g,%g)\n",
919 i, j, support[0], support[1], support[2]);
920 StatTotal++;
921 /* First, do a quick sanity check */
922 for (k = 0, ris = rjs = 0; k < DIMENSION; k++) {
923 ris += pow2(Atoms[i].x[k] - support[k]);
924 rjs += pow2(Atoms[j].x[k] - support[k]);
925 }
926 ris = sqrt(ris);
927 rjs = sqrt(rjs);
928 if (fabs(ris - rjs) > TolerancePrimary) {
929 StatEarly++;
930 if (verbose > 0)
931 printf(" Support can't actually define a rotation axis\n");
932 return NULL;
933 }
934 axis = alloc_symmetry_element();
936 axis->order = 2;
937 axis->nparam = 7;
938 for (k = 0, r = 0; k < DIMENSION; k++)
940 r = sqrt(r);
941 if (r > 0) {
942 for (k = 0; k < DIMENSION; k++)
943 axis->normal[k] = CenterOfSomething[k] / r;
944 } else {
945 axis->normal[0] = 1;
946 for (k = 1; k < DIMENSION; k++)
947 axis->normal[k] = 0;
948 }
949 axis->distance = r;
950 for (k = 0, r = 0; k < DIMENSION; k++) {
951 center[k] = (Atoms[i].x[k] + Atoms[j].x[k]) / 2 - support[k];
952 r += center[k] * center[k];
953 }
954 r = sqrt(r);
955 if (r <=
956 TolerancePrimary) { /* c2 is underdefined, let's do something special */
957 if (MolecularPlane != NULL) {
958 if (verbose > 0)
959 printf(" c2 is underdefined, but there is a molecular plane\n");
960 for (k = 0; k < DIMENSION; k++)
961 axis->direction[k] = MolecularPlane->normal[k];
962 } else {
963 if (verbose > 0)
964 printf(" c2 is underdefined, trying random direction\n");
965 for (k = 0; k < DIMENSION; k++)
966 center[k] = Atoms[i].x[k] - Atoms[j].x[k];
967 if (fabs(center[2]) + fabs(center[1]) > ToleranceSame) {
968 axis->direction[0] = 0;
969 axis->direction[1] = center[2];
970 axis->direction[2] = -center[1];
971 } else {
972 axis->direction[0] = -center[2];
973 axis->direction[1] = 0;
974 axis->direction[2] = center[0];
975 }
976 for (k = 0, r = 0; k < DIMENSION; k++)
977 r += axis->direction[k] * axis->direction[k];
978 r = sqrt(r);
979 for (k = 0; k < DIMENSION; k++)
980 axis->direction[k] /= r;
981 }
982 } else { /* direction is Ok, renormalize it */
983 for (k = 0; k < DIMENSION; k++)
984 axis->direction[k] = center[k] / r;
985 }
986 if (refine_symmetry_element(axis, 1) < 0) {
987 if (verbose > 0)
988 printf(" refinement failed for the c2 axis\n");
990 return NULL;
991 }
992 return axis;
993}
994
995SYMMETRY_ELEMENT *init_axis_parameters(double a[3], double b[3], double c[3]) {
996 SYMMETRY_ELEMENT *axis;
997 int i, order, sign;
998 double ra, rb, rc, rab, rbc, rac, r;
999 double angle;
1000
1001 ra = rb = rc = rab = rbc = rac = 0;
1002 for (i = 0; i < DIMENSION; i++) {
1003 ra += a[i] * a[i];
1004 rb += b[i] * b[i];
1005 rc += c[i] * c[i];
1006 }
1007 ra = sqrt(ra);
1008 rb = sqrt(rb);
1009 rc = sqrt(rc);
1010 if (fabs(ra - rb) > TolerancePrimary || fabs(ra - rc) > TolerancePrimary ||
1011 fabs(rb - rc) > TolerancePrimary) {
1012 StatEarly++;
1013 if (verbose > 0)
1014 printf(" points are not on a sphere\n");
1015 return NULL;
1016 }
1017 for (i = 0; i < DIMENSION; i++) {
1018 rab += (a[i] - b[i]) * (a[i] - b[i]);
1019 rac += (a[i] - c[i]) * (a[i] - c[i]);
1020 rbc += (c[i] - b[i]) * (c[i] - b[i]);
1021 }
1022 rab = sqrt(rab);
1023 rac = sqrt(rac);
1024 rbc = sqrt(rbc);
1025 if (fabs(rab - rbc) > TolerancePrimary) {
1026 StatEarly++;
1027 if (verbose > 0)
1028 printf(" points can't be rotation-equivalent\n");
1029 return NULL;
1030 }
1031 if (rab <= ToleranceSame || rbc <= ToleranceSame || rac <= ToleranceSame) {
1032 StatEarly++;
1033 if (verbose > 0)
1034 printf(" rotation is underdefined by these points\n");
1035 return NULL;
1036 }
1037 rab = (rab + rbc) / 2;
1038 angle = M_PI - 2 * asin(rac / (2 * rab));
1039 if (verbose > 1)
1040 printf(" rotation angle is %f\n", angle);
1041 if (fabs(angle) <= M_PI / (MaxAxisOrder + 1)) {
1042 StatEarly++;
1043 if (verbose > 0)
1044 printf(" atoms are too close to a straight line\n");
1045 return NULL;
1046 }
1047 order = floor((2 * M_PI) / angle + 0.5);
1048 if (order <= 2 || order > MaxAxisOrder) {
1049 StatEarly++;
1050 if (verbose > 0)
1051 printf(" rotation axis order (%d) is not from 3 to %d\n", order,
1052 MaxAxisOrder);
1053 return NULL;
1054 }
1055 axis = alloc_symmetry_element();
1056 axis->order = order;
1057 axis->nparam = 7;
1058 for (i = 0, r = 0; i < DIMENSION; i++)
1060 r = sqrt(r);
1061 if (r > 0) {
1062 for (i = 0; i < DIMENSION; i++)
1063 axis->normal[i] = CenterOfSomething[i] / r;
1064 } else {
1065 axis->normal[0] = 1;
1066 for (i = 1; i < DIMENSION; i++)
1067 axis->normal[i] = 0;
1068 }
1069 axis->distance = r;
1070 axis->direction[0] =
1071 (b[1] - a[1]) * (c[2] - b[2]) - (b[2] - a[2]) * (c[1] - b[1]);
1072 axis->direction[1] =
1073 (b[2] - a[2]) * (c[0] - b[0]) - (b[0] - a[0]) * (c[2] - b[2]);
1074 axis->direction[2] =
1075 (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]);
1076 /*
1077 * Arbitrarily select axis direction so that first non-zero component
1078 * of the direction is positive.
1079 */
1080 sign = 0;
1081 if (axis->direction[0] <= 0) {
1082 if (axis->direction[0] < 0) {
1083 sign = 1;
1084 } else if (axis->direction[1] <= 0) {
1085 if (axis->direction[1] < 0) {
1086 sign = 1;
1087 } else if (axis->direction[2] < 0) {
1088 sign = 1;
1089 }
1090 }
1091 }
1092 if (sign)
1093 for (i = 0; i < DIMENSION; i++)
1094 axis->direction[i] = -axis->direction[i];
1095 for (i = 0, r = 0; i < DIMENSION; i++)
1096 r += axis->direction[i] * axis->direction[i];
1097 r = sqrt(r);
1098 for (i = 0; i < DIMENSION; i++)
1099 axis->direction[i] /= r;
1100 if (verbose > 1) {
1101 printf(" axis origin is at (%g,%g,%g)\n",
1102 axis->normal[0] * axis->distance, axis->normal[1] * axis->distance,
1103 axis->normal[2] * axis->distance);
1104 printf(" axis is in the direction (%g,%g,%g)\n", axis->direction[0],
1105 axis->direction[1], axis->direction[2]);
1106 }
1107 return axis;
1108}
1110SYMMETRY_ELEMENT *init_higher_axis(int ia, int ib, int ic) {
1111 SYMMETRY_ELEMENT *axis;
1112 double a[DIMENSION], b[DIMENSION], c[DIMENSION];
1113 int i;
1114
1115 if (verbose > 0)
1116 printf("Trying cn axis for the triplet (%d,%d,%d)\n", ia, ib, ic);
1117 StatTotal++;
1118 /* Do a quick check of geometry validity */
1119 for (i = 0; i < DIMENSION; i++) {
1120 a[i] = Atoms[ia].x[i] - CenterOfSomething[i];
1121 b[i] = Atoms[ib].x[i] - CenterOfSomething[i];
1122 c[i] = Atoms[ic].x[i] - CenterOfSomething[i];
1123 }
1124 if ((axis = init_axis_parameters(a, b, c)) == NULL) {
1125 if (verbose > 0)
1126 printf(" no coherent axis is defined by the points\n");
1127 return NULL;
1128 }
1130 if (refine_symmetry_element(axis, 1) < 0) {
1131 if (verbose > 0)
1132 printf(" refinement failed for the c%d axis\n", axis->order);
1134 return NULL;
1135 }
1136 return axis;
1137}
1138
1139/*
1140 * Improper axes-specific routines.
1141 * These are obtained by slight modifications of normal rotation
1142 * routines.
1143 */
1145 double x[3], y[3], a[3], b[3], c[3];
1146 double angle = 2 * M_PI / axis->order;
1147 double a_sin = sin(angle);
1148 double a_cos = cos(angle);
1149 double dot;
1150 int i;
1151
1152 if (DIMENSION != 3) {
1153 fprintf(stderr, "Catastrophe in rotate_reflect_atom!\n");
1154 exit(EXIT_FAILURE);
1155 }
1156 for (i = 0; i < 3; i++)
1157 x[i] = from->x[i] - axis->distance * axis->normal[i];
1158 for (i = 0, dot = 0; i < 3; i++)
1159 dot += x[i] * axis->direction[i];
1160 for (i = 0; i < 3; i++)
1161 a[i] = axis->direction[i] * dot;
1162 for (i = 0; i < 3; i++)
1163 b[i] = x[i] - a[i];
1164 c[0] = b[1] * axis->direction[2] - b[2] * axis->direction[1];
1165 c[1] = b[2] * axis->direction[0] - b[0] * axis->direction[2];
1166 c[2] = b[0] * axis->direction[1] - b[1] * axis->direction[0];
1167 for (i = 0; i < 3; i++)
1168 y[i] = -a[i] + b[i] * a_cos + c[i] * a_sin;
1169 for (i = 0; i < 3; i++)
1170 to->x[i] = y[i] + axis->distance * axis->normal[i];
1171 to->type = from->type;
1172}
1173
1174SYMMETRY_ELEMENT *init_improper_axis(int ia, int ib, int ic) {
1175 SYMMETRY_ELEMENT *axis;
1176 double a[DIMENSION], b[DIMENSION], c[DIMENSION];
1177 double centerpoint[DIMENSION];
1178 double r;
1179 int i;
1181 if (verbose > 0)
1182 printf("Trying sn axis for the triplet (%d,%d,%d)\n", ia, ib, ic);
1183 StatTotal++;
1184 /* First, reduce the problem to Cn case */
1185 for (i = 0; i < DIMENSION; i++) {
1186 a[i] = Atoms[ia].x[i] - CenterOfSomething[i];
1187 b[i] = Atoms[ib].x[i] - CenterOfSomething[i];
1188 c[i] = Atoms[ic].x[i] - CenterOfSomething[i];
1189 }
1190 for (i = 0, r = 0; i < DIMENSION; i++) {
1191 centerpoint[i] = a[i] + c[i] + 2 * b[i];
1192 r += centerpoint[i] * centerpoint[i];
1193 }
1194 r = sqrt(r);
1195 if (r <= ToleranceSame) {
1196 StatEarly++;
1197 if (verbose > 0)
1198 printf(
1199 " atoms can not define improper axis of the order more than 2\n");
1200 return NULL;
1201 }
1202 for (i = 0; i < DIMENSION; i++)
1203 centerpoint[i] /= r;
1204 for (i = 0, r = 0; i < DIMENSION; i++)
1205 r += centerpoint[i] * b[i];
1206 for (i = 0; i < DIMENSION; i++)
1207 b[i] = 2 * r * centerpoint[i] - b[i];
1208 /* Do a quick check of geometry validity */
1209 if ((axis = init_axis_parameters(a, b, c)) == NULL) {
1210 if (verbose > 0)
1211 printf(" no coherent improper axis is defined by the points\n");
1212 return NULL;
1213 }
1215 if (refine_symmetry_element(axis, 1) < 0) {
1216 if (verbose > 0)
1217 printf(" refinement failed for the s%d axis\n", axis->order);
1219 return NULL;
1220 }
1221 return axis;
1222}
1223
1224/*
1225 * Control routines
1226 */
1227
1228void find_center_of_something(void) {
1229 int i, j;
1230 double coord_sum[DIMENSION];
1231 double r;
1232
1233 for (j = 0; j < DIMENSION; j++)
1234 coord_sum[j] = 0;
1235 for (i = 0; i < AtomsCount; i++) {
1236 for (j = 0; j < DIMENSION; j++)
1237 coord_sum[j] += Atoms[i].x[j];
1238 }
1239 for (j = 0; j < DIMENSION; j++)
1240 CenterOfSomething[j] = coord_sum[j] / AtomsCount;
1241 if (verbose > 0)
1242 printf("Center of something is at %15.10f, %15.10f, %15.10f\n",
1244 DistanceFromCenter = (double *)calloc(AtomsCount, sizeof(double));
1245 if (DistanceFromCenter == NULL) {
1246 fprintf(stderr, "Unable to allocate array for the distances\n");
1247 exit(EXIT_FAILURE);
1248 }
1249 for (i = 0; i < AtomsCount; i++) {
1250 for (j = 0, r = 0; j < DIMENSION; j++)
1251 r += pow2(Atoms[i].x[j] - CenterOfSomething[j]);
1252 DistanceFromCenter[i] = r;
1253 }
1254}
1255
1256void find_planes(void) {
1257 int i, j;
1258 SYMMETRY_ELEMENT *plane = NULL;
1259
1260 plane = init_ultimate_plane();
1261 if (plane != NULL) {
1262 MolecularPlane = plane;
1263 PlanesCount++;
1265 PlanesCount);
1266 if (Planes == NULL) {
1267 perror("Out of memory in find_planes");
1268 exit(EXIT_FAILURE);
1269 }
1270 Planes[PlanesCount - 1] = plane;
1271 }
1272 for (i = 1; i < AtomsCount; i++) {
1273 for (j = 0; j < i; j++) {
1274 if (Atoms[i].type != Atoms[j].type)
1275 continue;
1276 if ((plane = init_mirror_plane(i, j)) != NULL) {
1277 PlanesCount++;
1278 Planes = (SYMMETRY_ELEMENT **)realloc(
1279 Planes, sizeof(SYMMETRY_ELEMENT *) * PlanesCount);
1280 if (Planes == NULL) {
1281 perror("Out of memory in find_planes");
1282 exit(EXIT_FAILURE);
1283 }
1284 Planes[PlanesCount - 1] = plane;
1285 }
1286 }
1287 }
1288}
1289
1290void destroy_planes(void) {
1291 int i;
1292 for (i = 0; i < PlanesCount; i++)
1294 PlanesCount = 0;
1295 free(Planes);
1296 MolecularPlane = NULL;
1297 Planes = NULL;
1298}
1299
1300void find_inversion_centers(void) {
1302
1303 if ((center = init_inversion_center()) != NULL) {
1305 (SYMMETRY_ELEMENT **)calloc(1, sizeof(SYMMETRY_ELEMENT *));
1308 }
1309}
1310
1311void destroy_inversion_centers(void) {
1312 int i;
1313 for (i = 0; i < InversionCentersCount; i++)
1316 free(InversionCenters);
1317 InversionCenters = NULL;
1318}
1319
1320void find_infinity_axis(void) {
1321 SYMMETRY_ELEMENT *axis;
1322
1323 if ((axis = init_ultimate_axis()) != NULL) {
1325 NormalAxes = (SYMMETRY_ELEMENT **)realloc(
1327 if (NormalAxes == NULL) {
1328 perror("Out of memory in find_infinity_axes()");
1329 exit(EXIT_FAILURE);
1330 }
1331 NormalAxes[NormalAxesCount - 1] = axis;
1332 }
1333}
1334
1335void find_c2_axes(void) {
1336 int i, j, k, l, m;
1337 double center[DIMENSION];
1338 double *distances = calloc(AtomsCount, sizeof(double));
1339 double r;
1340 SYMMETRY_ELEMENT *axis;
1341
1342 if (distances == NULL) {
1343 fprintf(stderr, "Out of memory in find_c2_axes()\n");
1344 exit(EXIT_FAILURE);
1345 }
1346 for (i = 1; i < AtomsCount; i++) {
1347 for (j = 0; j < i; j++) {
1348 if (Atoms[i].type != Atoms[j].type)
1349 continue;
1352 continue; /* A very cheap, but quite effective check */
1353 /*
1354 * First, let's try to get it cheap and use CenterOfSomething
1355 */
1356 for (k = 0, r = 0; k < DIMENSION; k++) {
1357 center[k] = (Atoms[i].x[k] + Atoms[j].x[k]) / 2;
1358 r += pow2(center[k] - CenterOfSomething[k]);
1359 }
1360 r = sqrt(r);
1361 if (r > 5 * TolerancePrimary) { /* It's Ok to use CenterOfSomething */
1362 if ((axis = init_c2_axis(i, j, CenterOfSomething)) != NULL) {
1364 NormalAxes = (SYMMETRY_ELEMENT **)realloc(
1366 if (NormalAxes == NULL) {
1367 perror("Out of memory in find_c2_axes");
1368 exit(EXIT_FAILURE);
1369 }
1370 NormalAxes[NormalAxesCount - 1] = axis;
1371 }
1372 continue;
1373 }
1374 /*
1375 * Now, C2 axis can either pass through an atom, or through the
1376 * middle of the other pair.
1377 */
1378 for (k = 0; k < AtomsCount; k++) {
1379 if ((axis = init_c2_axis(i, j, Atoms[k].x)) != NULL) {
1381 NormalAxes = (SYMMETRY_ELEMENT **)realloc(
1383 if (NormalAxes == NULL) {
1384 perror("Out of memory in find_c2_axes");
1385 exit(EXIT_FAILURE);
1386 }
1387 NormalAxes[NormalAxesCount - 1] = axis;
1388 }
1389 }
1390 /*
1391 * Prepare data for an additional pre-screening check
1392 */
1393 for (k = 0; k < AtomsCount; k++) {
1394 for (l = 0, r = 0; l < DIMENSION; l++)
1395 r += pow2(Atoms[k].x[l] - center[l]);
1396 distances[k] = sqrt(r);
1397 }
1398 for (k = 0; k < AtomsCount; k++) {
1399 for (l = 0; l < AtomsCount; l++) {
1400 if (Atoms[k].type != Atoms[l].type)
1401 continue;
1404 fabs(distances[k] - distances[l]) > TolerancePrimary)
1405 continue; /* We really need this one to run reasonably fast! */
1406 for (m = 0; m < DIMENSION; m++)
1407 center[m] = (Atoms[k].x[m] + Atoms[l].x[m]) / 2;
1408 if ((axis = init_c2_axis(i, j, center)) != NULL) {
1410 NormalAxes = (SYMMETRY_ELEMENT **)realloc(
1412 if (NormalAxes == NULL) {
1413 perror("Out of memory in find_c2_axes");
1414 exit(EXIT_FAILURE);
1415 }
1416 NormalAxes[NormalAxesCount - 1] = axis;
1417 }
1418 }
1419 }
1420 }
1421 }
1422 free(distances);
1423}
1424
1425void destroy_normal_axes() {
1426 int i = 0;
1427 for (i = 0; i < NormalAxesCount; i++)
1429 NormalAxesCount = 0;
1430 free(NormalAxes);
1431 NormalAxes = NULL;
1432}
1433
1434void destroy_improper_axes() {
1435 int i = 0;
1436 for (i = 0; i < ImproperAxesCount; i++)
1439 free(ImproperAxes);
1440 ImproperAxes = NULL;
1441}
1442
1443void find_higher_axes(void) {
1444 int i, j, k;
1445 SYMMETRY_ELEMENT *axis;
1446
1447 for (i = 0; i < AtomsCount; i++) {
1448 for (j = i + 1; j < AtomsCount; j++) {
1449 if (Atoms[i].type != Atoms[j].type)
1450 continue;
1453 continue; /* A very cheap, but quite effective check */
1454 for (k = 0; k < AtomsCount; k++) {
1455 if (Atoms[i].type != Atoms[k].type)
1456 continue;
1461 continue;
1462 if ((axis = init_higher_axis(i, j, k)) != NULL) {
1464 NormalAxes = (SYMMETRY_ELEMENT **)realloc(
1466 if (NormalAxes == NULL) {
1467 perror("Out of memory in find_higher_axes");
1468 exit(EXIT_FAILURE);
1469 }
1470 NormalAxes[NormalAxesCount - 1] = axis;
1471 }
1472 }
1473 }
1474 }
1475}
1476
1477void find_improper_axes(void) {
1478 int i, j, k;
1479 SYMMETRY_ELEMENT *axis;
1480
1481 for (i = 0; i < AtomsCount; i++) {
1482 for (j = i + 1; j < AtomsCount; j++) {
1483 for (k = 0; k < AtomsCount; k++) {
1484 if ((axis = init_improper_axis(i, j, k)) != NULL) {
1486 ImproperAxes = (SYMMETRY_ELEMENT **)realloc(
1488 if (ImproperAxes == NULL) {
1489 perror("Out of memory in find_higher_axes");
1490 exit(EXIT_FAILURE);
1491 }
1492 ImproperAxes[ImproperAxesCount - 1] = axis;
1493 }
1494 }
1495 }
1496 }
1497}
1498
1499void report_planes(void) {
1500 int i;
1501
1502 if (PlanesCount == 0)
1503 printf("There are no planes of symmetry in the molecule\n");
1504 else {
1505 if (PlanesCount == 1)
1506 printf("There is a plane of symmetry in the molecule\n");
1507 else
1508 printf("There are %d planes of symmetry in the molecule\n", PlanesCount);
1509 printf(
1510 " Residual Direction of the normal Distance\n");
1511 for (i = 0; i < PlanesCount; i++) {
1512 printf("%3d %8.4e ", i, Planes[i]->maxdev);
1513 printf("(%11.8f,%11.8f,%11.8f) ", Planes[i]->normal[0],
1514 Planes[i]->normal[1], Planes[i]->normal[2]);
1515 printf("%14.8f\n", Planes[i]->distance);
1516 }
1517 }
1518}
1519
1520void report_inversion_centers(void) {
1521 if (InversionCentersCount == 0)
1522 printf("There is no inversion center in the molecule\n");
1523 else {
1524 printf("There in an inversion center in the molecule\n");
1525 printf(" Residual Position\n");
1526 printf(" %8.4e ", InversionCenters[0]->maxdev);
1527 printf("(%14.8f,%14.8f,%14.8f)\n",
1528 InversionCenters[0]->distance * InversionCenters[0]->normal[0],
1529 InversionCenters[0]->distance * InversionCenters[0]->normal[1],
1530 InversionCenters[0]->distance * InversionCenters[0]->normal[2]);
1531 }
1532}
1533
1534void report_axes(void) {
1535 int i;
1536
1537 if (NormalAxesCount == 0)
1538 printf("There are no normal axes in the molecule\n");
1539 else {
1540 if (NormalAxesCount == 1)
1541 printf("There is a normal axis in the molecule\n");
1542 else
1543 printf("There are %d normal axes in the molecule\n", NormalAxesCount);
1544 printf(" Residual Order Direction of the axis "
1545 " Supporting point\n");
1546 for (i = 0; i < NormalAxesCount; i++) {
1547 printf("%3d %8.4e ", i, NormalAxes[i]->maxdev);
1548 if (NormalAxes[i]->order == 0)
1549 printf("Inf ");
1550 else
1551 printf("%3d ", NormalAxes[i]->order);
1552 printf("(%11.8f,%11.8f,%11.8f) ", NormalAxes[i]->direction[0],
1553 NormalAxes[i]->direction[1], NormalAxes[i]->direction[2]);
1554 printf("(%14.8f,%14.8f,%14.8f)\n",
1555 NormalAxes[0]->distance * NormalAxes[0]->normal[0],
1556 NormalAxes[0]->distance * NormalAxes[0]->normal[1],
1557 NormalAxes[0]->distance * NormalAxes[0]->normal[2]);
1558 }
1559 }
1560}
1561
1562void report_improper_axes(void) {
1563 int i;
1564
1565 if (ImproperAxesCount == 0)
1566 printf("There are no improper axes in the molecule\n");
1567 else {
1568 if (ImproperAxesCount == 1)
1569 printf("There is an improper axis in the molecule\n");
1570 else
1571 printf("There are %d improper axes in the molecule\n", ImproperAxesCount);
1572 printf(" Residual Order Direction of the axis "
1573 " Supporting point\n");
1574 for (i = 0; i < ImproperAxesCount; i++) {
1575 printf("%3d %8.4e ", i, ImproperAxes[i]->maxdev);
1576 if (ImproperAxes[i]->order == 0)
1577 printf("Inf ");
1578 else
1579 printf("%3d ", ImproperAxes[i]->order);
1580 printf("(%11.8f,%11.8f,%11.8f) ", ImproperAxes[i]->direction[0],
1581 ImproperAxes[i]->direction[1], ImproperAxes[i]->direction[2]);
1582 printf("(%14.8f,%14.8f,%14.8f)\n",
1583 ImproperAxes[0]->distance * ImproperAxes[0]->normal[0],
1584 ImproperAxes[0]->distance * ImproperAxes[0]->normal[1],
1585 ImproperAxes[0]->distance * ImproperAxes[0]->normal[2]);
1586 }
1587 }
1588}
1589
1590/*
1591 * General symmetry handling
1592 */
1593void report_and_reset_counters(void) {
1594 printf(" %10ld candidates examined\n"
1595 " %10ld removed early\n"
1596 " %10ld removed during initial mating stage\n"
1597 " %10ld removed as duplicates\n"
1598 " %10ld removed because of the wrong transformation order\n"
1599 " %10ld removed after unsuccessful optimization\n"
1600 " %10ld accepted\n",
1602 StatAccept);
1604 StatAccept = 0;
1605}
1606
1607void find_symmetry_elements(void) {
1609 if (verbose > -1) {
1610 printf("Looking for the inversion center\n");
1611 }
1613 if (verbose > -1) {
1615 printf("Looking for the planes of symmetry\n");
1616 }
1617 find_planes();
1618 if (verbose > -1) {
1620 printf("Looking for infinity axis\n");
1621 }
1623 if (verbose > -1) {
1625 printf("Looking for C2 axes\n");
1626 }
1627 find_c2_axes();
1628 if (verbose > -1) {
1630 printf("Looking for higher axes\n");
1631 }
1633 if (verbose > -1) {
1635 printf("Looking for the improper axes\n");
1636 }
1638 if (verbose > -1) {
1640 }
1641}
1642
1643int compare_axes(const void *a, const void *b) {
1644 SYMMETRY_ELEMENT *axis_a = *(SYMMETRY_ELEMENT **)a;
1645 SYMMETRY_ELEMENT *axis_b = *(SYMMETRY_ELEMENT **)b;
1646 int i, order_a, order_b;
1647
1648 order_a = axis_a->order;
1649 if (order_a == 0)
1650 order_a = 10000;
1651 order_b = axis_b->order;
1652 if (order_b == 0)
1653 order_b = 10000;
1654 if ((i = order_b - order_a) != 0)
1655 return i;
1656 if (axis_a->maxdev > axis_b->maxdev)
1657 return -1;
1658 if (axis_a->maxdev < axis_b->maxdev)
1659 return 1;
1660 return 0;
1661}
1662
1663void sort_symmetry_elements(void) {
1664 if (PlanesCount > 1) {
1665 qsort(Planes, PlanesCount, sizeof(SYMMETRY_ELEMENT *), compare_axes);
1666 }
1667 if (NormalAxesCount > 1) {
1669 compare_axes);
1670 }
1671 if (ImproperAxesCount > 1) {
1673 compare_axes);
1674 }
1675}
1676
1679 report_axes();
1682}
1683
1684void summarize_symmetry_elements(void) {
1685 int i;
1686
1687 NormalAxesCounts = (int *)calloc(MaxAxisOrder + 1, sizeof(int));
1688 ImproperAxesCounts = (int *)calloc(MaxAxisOrder + 1, sizeof(int));
1689 for (i = 0; i < NormalAxesCount; i++)
1690 NormalAxesCounts[NormalAxes[i]->order]++;
1691 for (i = 0; i < ImproperAxesCount; i++)
1693}
1694
1696 int i;
1697 char *symmetry_code =
1698 calloc(1, 10 * (PlanesCount + NormalAxesCount + ImproperAxesCount +
1700 char buf[100];
1702 if (symmetry_code == NULL) {
1703 fprintf(stderr, "Unable to allocate memory for symmetry ID code in "
1704 "report_symmetry_elements_brief()\n");
1705 exit(EXIT_FAILURE);
1706 }
1709 0) {
1710 if (verbose > -1)
1711 printf("Molecule has no symmetry elements\n");
1712 } else {
1713 if (verbose > -1)
1714 printf("Molecule has the following symmetry elements: ");
1715 if (InversionCentersCount > 0)
1716 strcat(symmetry_code, "(i) ");
1717 if (NormalAxesCounts[0] == 1)
1718 strcat(symmetry_code, "(Cinf) ");
1719 if (NormalAxesCounts[0] > 1) {
1720 sprintf(buf, "%d*(Cinf) ", NormalAxesCounts[0]);
1721 strcat(symmetry_code, buf);
1722 }
1723 for (i = MaxAxisOrder; i >= 2; i--) {
1724 if (NormalAxesCounts[i] == 1) {
1725 sprintf(buf, "(C%d) ", i);
1726 strcat(symmetry_code, buf);
1727 }
1728 if (NormalAxesCounts[i] > 1) {
1729 sprintf(buf, "%d*(C%d) ", NormalAxesCounts[i], i);
1730 strcat(symmetry_code, buf);
1731 }
1732 }
1733 for (i = MaxAxisOrder; i >= 2; i--) {
1734 if (ImproperAxesCounts[i] == 1) {
1735 sprintf(buf, "(S%d) ", i);
1736 strcat(symmetry_code, buf);
1737 }
1738 if (ImproperAxesCounts[i] > 1) {
1739 sprintf(buf, "%d*(S%d) ", ImproperAxesCounts[i], i);
1740 strcat(symmetry_code, buf);
1741 }
1742 }
1743 if (PlanesCount == 1)
1744 strcat(symmetry_code, "(sigma) ");
1745 if (PlanesCount > 1) {
1746 sprintf(buf, "%d*(sigma) ", PlanesCount);
1747 strcat(symmetry_code, buf);
1748 }
1749 if (verbose > -1)
1750 printf("%s\n", symmetry_code);
1751 }
1752 SymmetryCode = symmetry_code;
1753}
1754
1755void identify_point_group(int *point_group) {
1756 int i;
1757 int last_matching = -1;
1758 int matching_count = 0;
1759
1760 for (i = 0; i < PointGroupsCount; i++) {
1761 if (strcmp(SymmetryCode, PointGroups[i].symmetry_code) == 0) {
1762 if (PointGroups[i].check() == 1) {
1763 last_matching = i;
1764 matching_count++;
1765 } else {
1766 if (verbose > -2) {
1767 printf("It looks very much like %s, but it is not since %s\n",
1769 }
1770 }
1771 }
1772 }
1773 if (matching_count == 0) {
1774 printf("These symmetry elements match no point group I know of. Sorry.\n");
1775 }
1776 if (matching_count > 1) {
1777 printf("These symmetry elements match more than one group I know of.\n"
1778 "SOMETHING IS VERY WRONG\n");
1779 printf("Matching groups are:\n");
1780 for (i = 0; i < PointGroupsCount; i++) {
1781 if ((strcmp(SymmetryCode, PointGroups[i].symmetry_code) == 0) &&
1782 (PointGroups[i].check() == 1)) {
1783 printf(" %s\n", PointGroups[i].group_name);
1784 }
1785 }
1786 }
1787 if (verbose > 0 && matching_count == 1) {
1788 printf("It seems to be the %s point group\n",
1789 PointGroups[last_matching].group_name);
1790 }
1791
1792 *point_group = last_matching;
1793}
1794
1795/*
1796 * Input/Output
1797 */
1798
1799void FC_FUNC_(symmetries_finite_init,
1800 SYMMETRIES_FINITE_INIT)(const int *natoms, const int *types,
1801 const double *positions,
1802 const int *verbosity, int *point_group) {
1803 int i;
1804
1805 verbose = *verbosity;
1806 AtomsCount = *natoms;
1807 Atoms = calloc(AtomsCount, sizeof(ATOM));
1808 for (i = 0; i < AtomsCount; i++) {
1809 Atoms[i].type = types[i];
1810 Atoms[i].x[0] = positions[3 * i + 0];
1811 Atoms[i].x[1] = positions[3 * i + 1];
1812 Atoms[i].x[2] = positions[3 * i + 2];
1813 if (verbose > 5)
1814 printf("%d %f %f %f\n", Atoms[i].type, Atoms[i].x[0], Atoms[i].x[1],
1815 Atoms[i].x[2]);
1816 }
1817
1818 if (AtomsCount > 0)
1821
1823 if (BadOptimization)
1824 printf("Refinement of some symmetry elements was terminated before "
1825 "convergence was reached.\n"
1826 "Some symmetry elements may remain unidentified.\n");
1827 if (verbose >= 0)
1830 identify_point_group(point_group);
1831}
1832
1833void symmetries_finite_get_group_name(const int *point_group,
1834 char * name) {
1835
1836 strcpy(name, PointGroups[*point_group].group_name);
1837}
1838
1840 const int *point_group, char * elements) {
1841
1842 strncpy(elements, PointGroups[*point_group].symmetry_code, strlen(elements));
1843}
1844
1845void FC_FUNC_(symmetries_finite_end, SYMMETRIES_FINITE_END)() {
1850
1851 free(SymmetryCode);
1852 SymmetryCode = NULL;
1853
1854 free(DistanceFromCenter);
1855 DistanceFromCenter = NULL;
1856
1857 free(NormalAxesCounts);
1858 NormalAxesCounts = NULL;
1859
1860 free(ImproperAxesCounts);
1861 ImproperAxesCounts = NULL;
1862
1863 free(Atoms);
1864 Atoms = NULL;
1865}
void * memcpy(void *__restrict __dest, const void *__restrict __src, size_t __n) __attribute__((__nothrow__
!in the code
Definition: global.h:59
pure real(real64) function center(this)
Center of the filter interval.
real(real64) function s()
ptrdiff_t l
Definition: operate_inc.c:12
ptrdiff_t i
Definition: operate_inc.c:12
ptrdiff_t j
Definition: operate_inc.c:12
for(l=0;l< nri;l++)
Definition: operate_inc.c:15
static double f(double w, void *p)
void(* transform_atom)(struct _SYMMETRY_ELEMENT_ *el, ATOM *from, ATOM *to)
SYMMETRY_ELEMENT * init_inversion_center(void)
char * PointGroupRejectionReason
void set_params(SYMMETRY_ELEMENT *elem, double values[])
SYMMETRY_ELEMENT * init_axis_parameters(double a[3], double b[3], double c[3])
void find_center_of_something(void)
static double OptChangeThreshold
SYMMETRY_ELEMENT * init_c2_axis(int i, int j, double support[3])
struct _SYMMETRY_ELEMENT_ SYMMETRY_ELEMENT
static int verbose
void report_improper_axes(void)
SYMMETRY_ELEMENT * init_mirror_plane(int i, int j)
int establish_pairs(SYMMETRY_ELEMENT *elem)
static SYMMETRY_ELEMENT ** NormalAxes
static long StatPairs
int compare_axes(const void *a, const void *b)
static int PlanesCount
void FC_FUNC_(symmetries_finite_init, SYMMETRIES_FINITE_INIT) const
void destroy_inversion_centers(void)
void sort_symmetry_elements(void)
void find_improper_axes(void)
void find_inversion_centers(void)
double pow2(double x)
void symmetries_finite_get_group_elements(const int *point_group, char *elements)
static double * DistanceFromCenter
double sin(double __x) __attribute__((__nothrow__
int check_transform_order(SYMMETRY_ELEMENT *elem)
static double TolerancePrimary
static int MaxAxisOrder
static int NormalAxesCount
void find_higher_axes(void)
void destroy_symmetry_element(SYMMETRY_ELEMENT *elem)
static long StatOrder
static double ToleranceSame
static long StatAccept
void symmetries_finite_get_group_name(const int *point_group, char *name)
static double MaxOptStep
void destroy_normal_axes()
static int ImproperAxesCount
double asin(double __x) __attribute__((__nothrow__
static char * SymmetryCode
void find_symmetry_elements(void)
double fabs(double __x) __attribute__((__nothrow__
void find_c2_axes(void)
static double GradientStep
void find_infinity_axis(void)
SYMMETRY_ELEMENT * alloc_symmetry_element(void)
static SYMMETRY_ELEMENT * MolecularPlane
void optimize_transformation_params(SYMMETRY_ELEMENT *elem)
static long StatOpt
POINT_GROUP PointGroups[]
static SYMMETRY_ELEMENT ** InversionCenters
void invert_atom(SYMMETRY_ELEMENT *center, ATOM *from, ATOM *to)
static ATOM * Atoms
static double ToleranceFinal
SYMMETRY_ELEMENT * init_improper_axis(int ia, int ib, int ic)
static double CenterOfSomething[3]
double sqrt(double __x) __attribute__((__nothrow__
static double MinOptStep
void report_symmetry_elements_brief()
double eval_optimization_target_function(SYMMETRY_ELEMENT *elem, int *finish)
void summarize_symmetry_elements(void)
void report_planes(void)
SYMMETRY_ELEMENT * init_ultimate_plane(void)
void destroy_improper_axes()
void report_symmetry_elements_verbose(void)
static int AtomsCount
static int BadOptimization
int check_transform_quality(SYMMETRY_ELEMENT *elem)
static int * NormalAxesCounts
static SYMMETRY_ELEMENT ** Planes
void rotate_atom(SYMMETRY_ELEMENT *axis, ATOM *from, ATOM *to)
double cos(double __x) __attribute__((__nothrow__
void report_and_reset_counters(void)
void destroy_planes(void)
FILE * stderr
void identify_point_group(int *point_group)
void report_inversion_centers(void)
static int OptChangeHits
void rotate_reflect_atom(SYMMETRY_ELEMENT *axis, ATOM *from, ATOM *to)
static int MaxOptCycles
void get_params(SYMMETRY_ELEMENT *elem, double values[])
double floor(double __x) __attribute__((__nothrow__
static long StatEarly
int same_transform(SYMMETRY_ELEMENT *a, SYMMETRY_ELEMENT *b)
void find_planes(void)
static long StatTotal
void mirror_atom(SYMMETRY_ELEMENT *plane, ATOM *from, ATOM *to)
static int InversionCentersCount
void report_axes(void)
SYMMETRY_ELEMENT * init_ultimate_axis(void)
static long StatDups
static int * ImproperAxesCounts
SYMMETRY_ELEMENT * init_higher_axis(int ia, int ib, int ic)
static SYMMETRY_ELEMENT ** ImproperAxes
int refine_symmetry_element(SYMMETRY_ELEMENT *elem, int build_table)
real(real64) function values(xx)
Definition: test.F90:2025
static var_type * vars
Definition: varinfo_low.c:2760
void distance(const int iatom, const int jatom, const double coordinates[], double *rr, double *rr2, double *rr6, double *rr7)
Definition: vdw_ts_low.c:2059