1 // $Id: DRAWxtl1.cxx 1114 2011-02-23 20:29:18Z martin $
2 //
3 // module drawxtl1.cxx - part of DRAWxtl V5.5 - the GUI version
4 // Coded using the FLTK 1.1.6 widget set
5 //
6 // This program is copyrighted by Larry W. Finger, Martin Kroeker and Brian Toby
7 // and is distributed under the GNU General Public License - see the
8 // accompanying COPYING file for more details.
9 //
10 // This module contains the following routines:
11 //
12 // add_to_list - add atom position to unit-cell list
13 // add_vert - add vertices (atom positions) to display lists
14 // add_vert_nc - add vertices without checking
15 // analyze_bonds - print the bond distance tables
16 // axeqb - solve matrix equation AX = B for B (A is 3x3 only)
17 // build_box_contents - find all atoms in the display parallelopiped
18 // Calc_Rot - calculate rotation matrices
19 // convert_ellipsoid - handles conversions from Bij or Uij to betaij
20 // Conv_Sym_Mat - converts rotational part of symmetry matrix to a Cartesian rotation
21 // dist - calculates distance between two vertices
22 // dot0_3d - calculated dot product of two 3D vectors
23 // eigen - calculates eigen values and vectors for ellipsoids
24 // find_all_in_box - finds all atoms of a given type in the display box
25 // generate_arrows - generates magnetic vector arrows
26 // generate_bonds - Generate bond descriptions for the display lists
27 // generate_cones - generate lone-pair cone descriptions for display lists
28 // get_atom_id - reworks atom labels
29 // get_input - calls actual routine that reads input, generates lattice metric and does preliminary ellipsoid processing
30 // Locate_Triple - routine to place unit-cell vector triple
31 // make_bmat - calculates lattice metric routines for conversion to cartesian coordinates
32 // modulate_parameters - do the adjustments to positions and occupancies for aperiodic xtals
33 // modulate_uij - do the adjustments to Uij's for aperiodic crystals
34 // polygon_normal_3d - calculates normal vector of a polygon
35 // polygon_solid_angle_3d - calculates projected solid angle of a 3d plane polygon
36 // not_in_slab - check if a point is inside a given parallelepiped
37 // Output_Spheres - adds sphere descriptions to output display lists
38 // plot_vrml_poly - generates polyhedral descriptions to VRML output
39 // print_sym - print symmetry operators
40
41 #include "drawxtl.h"
42 #include "DRAWxtlViewUI.h"
43 #include <stdio.h>
44 #include <string.h>
45 #include <math.h>
46 #include <stdlib.h>
47 #include <FL/glut.H>
48 #if defined(__APPLE__)
49 # include <openGL/glu.h>
50 #else
51 # include <GL/glu.h>
52 #endif
53 #ifdef WIN32
54 #define cbrt(val) pow((val), 1.0 / 3.0)
55 #endif
56 // forward references for routines
57
58 #include "DRAWxtl_proto.h"
59
60 // global variables
61
62 #include "draw_ext.h"
63
64 /* ************************************************************** */
65 /* ************************************************************** */
66
67 void
add_to_list(float xp[3],int no,int nos)68 add_to_list (float xp[3], int no, int nos)
69
70 /* add atom position to unit-cell list if unique */
71 /* no - symmetry operator number */
72 /* nos - signed symmetry operator number to track inversion */
73 {
74 int i, ind;
75
76 for (i = 0; i <= 2; ++i) { /* get coordinates between 0 and 1 */
77 if (xp[i] >= 1.0)
78 xp[i] = xp[i] - 1.0f;
79 if (xp[i] < 0.0)
80 xp[i] = xp[i] + 1.0f;
81 }
82 ind = 0;
83 for (i = 0; i < ncell; ++i) {
84 if ((fabs (drvui->cell_xyz[i][0] - xp[0]) <= 0.00002) &&
85 (fabs (drvui->cell_xyz[i][1] - xp[1]) <= 0.00002) &&
86 (fabs (drvui->cell_xyz[i][2] - xp[2]) <= 0.00002))
87 ind = 1; /*set if same */
88 }
89 if (ind == 0) { /* true for new position */
90 for (i = 0; i <= 2; ++i) /*copy to list */
91 drvui->cell_xyz[ncell][i] = xp[i];
92 drvui->sym_op_signed[ncell] = nos;
93 drvui->sym_op_no[ncell++] = no;
94 }
95 }
96
97 /* ************************************************************** */
98 /* ************************************************************** */
99
100 void
add_vert(float vert[3],int no,int check,int modulate,int sign)101 add_vert (float vert[3], int no, int check, int modulate, int sign)
102
103 /* procedure to add the vertex in 'vert' to the current list if the
104 vertex is contained within the limits of the box, within the limits
105 of a polyhedral center, or needed to complete a molecule. If modulate
106 is true (!= 0), the coordinates in vert will be adjusted for modulation
107 before testing */
108 {
109 float c_vert[3]; /* vertex in cartesian coords */
110
111 float r_vert[3]; /* vertex in rotated coordinates */
112
113 float vert_saved[3]; /* save the unmodulated vetex */
114
115 int i, j; /* loop variables */
116
117 float temp; /* temporary storage */
118
119 int outside; /* temporary for molecular completion */
120
121 int number; /* decoded atom number */
122
123 int poly_center; /* true if this atom is center of a polyhedron */
124
125 int consider_extra;
126
127 float d = 0.0f;
128
129 double occupancy;
130
131 int sym_no;
132
133 if (!check_vert_alloc (NvertM, 0))
134 return; //abort immediately if out of array space
135 number = no / 1000;
136 sym_no = no - 1000 * number;
137 for (i = 0; i < 3; i++)
138 vert_saved[i] = vert[i];
139 if (modulate & (drvui->modulated > 0)) {
140 modulate_parameters (vert, &occupancy, sym_no, number);
141 if (occupancy < 0.01)
142 return;
143 }
144 for (i = 0; i < 3; ++i) { /* convert vertex coordinates to Cartesian */
145 c_vert[i] = 0.0f;
146 for (j = 0; j < 3; ++j)
147 c_vert[i] += (float) drvui->b_mat[i][j] * (vert[j] - origin[j]);
148 }
149 if (number >= 0 && number < drvui->atom_alloc)
150 poly_center = (((drvui->atoms[number].atom_n >> 8) & 255) > 0);
151 else
152 poly_center = 0;
153 consider_extra = (domolcomp != 0 || (drvui->npoly > 1 && !poly_center));
154 outside = 0;
155 if (check == 0) {
156 for (i = 0; i <= 2; i++) { /* check packing range */
157 if ((vert[i] < drvui->frames[drvui->frame_no].cryst_lim[i]) ||
158 (vert[i] > drvui->frames[drvui->frame_no].cryst_lim[i + 3]))
159 outside = 1;
160 }
161 if (outside == 0 && slabmode == 1)
162 outside = not_in_slab (c_vert[0], c_vert[1], c_vert[2]);
163
164 if (((outside == 0) && (fabs (c_vert[0]) > boxlim[0])) ||
165 (fabs (c_vert[1]) > boxlim[1]) || (fabs (c_vert[2]) > boxlim[2]))
166 outside = 1; /* check box limits */
167 if (outside != 0 && !consider_extra)
168 return;
169 }
170 d = 0.0f;
171 if (outside == 1) { /* molecule or polyhedral completion is in progress */
172 float dmin, dmax;
173
174 int polyno; /* number of the polyhedron */
175
176 for (i = 1; i < NvertM; i++) {
177 if (drvui->atoms[drvui->atom_no[i] / 1000].atom_fn != drvui->frame_no)
178 continue;
179
180 dmax = drvui->mol_d * drvui->mol_d;
181 dmin = 0.;
182 d = 0.0f;
183 for (j = 0; j < 3; j++) {
184 temp = (s_vert[3 * i + j] - c_vert[j]);
185 d += temp * temp;
186 }
187 j = drvui->atom_no[i] / 1000;
188 poly_center = (((drvui->atoms[j].atom_n >> 8) & 255) > 0);
189 if (poly_center) {
190 polyno = (drvui->atoms[j].atom_n >> 8) & 255;
191 dmax =
192 drvui->polyhedra[polyno].poly_size *
193 drvui->polyhedra[polyno].poly_size;
194 dmin =
195 drvui->polyhedra[polyno].poly_min * drvui->polyhedra[polyno].poly_min;
196 }
197 if ((domolcomp != 0 || poly_center) && (d <= dmax && d > 0.02)) {
198 outside = 2;
199 if (poly_center && d < dmin)
200 outside = 1;
201 if ((fabs (c_vert[0]) > 10. * boxlim[0])
202 || (fabs (c_vert[1]) > 10. * boxlim[1])
203 || (fabs (c_vert[2]) > 10. * boxlim[2])) {
204 fprintf (drvui->fcns,
205 "mol'completing beyond 10*lattice constant - aborting run\n");
206 outside = 1; /* far out of check box limits - combinatorial explosion */
207 break;
208 }
209 }
210 }
211 }
212 if (outside == 1)
213 return;
214 i = (int) (no / 1000);
215 for (j = 1; j < NvertM; ++j) {
216 if (fabs (s_vert[3 * j] - c_vert[0]) < 0.01
217 && fabs (s_vert[3 * j + 1] - c_vert[1]) < 0.01
218 && fabs (s_vert[3 * j + 2] - c_vert[2]) < 0.01
219 && check_atom_name (drvui->atoms[i].atom_l,
220 drvui->atoms[drvui->atom_no[j] / 1000].atom_l)
221 && (number == i))
222 return; /* eliminate duplicates */
223 }
224 for (i = 0; i <= 2; ++i) { /* add atom to Master Lists */
225 xypos[3 * NvertM + i] = vert[i];
226 xypos_nm[3 * NvertM + i] = vert_saved[i];
227 s_vert[3 * NvertM + i] = c_vert[i];
228 // if (offset[i] > c_vert[i])
229 // offset[i] = c_vert[i];
230 }
231 i = no - 1000 * (no / 1000);
232
233 //if (occupancy != 1.) fprintf (stderr,"occ=%f\n",occupancy);
234 drvui->vert_occ[NvertM] = (float) occupancy;
235 drvui->atom_so[NvertM] = sign;
236 drvui->atom_no[NvertM++] = no;
237 if (!check_vert_alloc (NvertM, 0)) {
238 if (domolcomp != 0) {
239 Error_Box ("Too many vertices for dimensions:\n"
240 " Did you use a 'molcomp' command incorrectly?\n"
241 " If not, please increase parameter MAX_VERTS.");
242 } else {
243 Error_Box
244 ("Too many vertices for dimensions - please increase parameter MAX_VERTS.");
245 }
246 return;
247 }
248 for (i = 0; i <= 2; ++i) { /* calculate position of point after POV rotation */
249 r_vert[i] = 0.0f;
250 for (j = 0; j <= 2; ++j)
251 r_vert[i] += (float) G_Rot[j][i] * c_vert[j];
252 }
253 if (drvui->automation)
254 return; /* if automation in progress, keep old scaling */
255 for (i = 0; i <= 2; ++i) { /* Update Min and Max of output */
256 if (POV_Max[i] < r_vert[i])
257 POV_Max[i] = r_vert[i];
258 if (POV_Min[i] > r_vert[i])
259 POV_Min[i] = r_vert[i];
260 }
261 }
262
263 /* ************************************************************** */
264 /* ************************************************************** */
265
266 void
add_vert_nc(float vert[3])267 add_vert_nc (float vert[3])
268
269 /* procedure to add the vertices in 'vert' to the current working list
270 (NOT Master List ) without checking - NO modulation */
271 {
272
273 int i, j; /* loop variables */
274
275 if (!s_vert) {
276 if ((s_vert = (float *) zalloc ((unsigned) (3 * (long) (drvui->verts_alloc
277 * sizeof (float))))) ==
278 NULL) {
279 Error_Box ("Unable to get initial s_vert allocation");
280 return;
281 }
282 }
283 for (i = 0; i <= 2; ++i) { /* convert vertex coordinates to Cartesian */
284 s_vert[3 * nvert + i] = 0.0f;
285 for (j = 0; j <= 2; ++j)
286 s_vert[3 * nvert + i] += (float) drvui->b_mat[i][j] * (vert[j] - origin[j]);
287 }
288 nvert++;
289 }
290
291 /* ************************************************************** */
292 /* ************************************************************** */
293
294 void
analyze_bonds(void)295 analyze_bonds (void)
296
297 /* routine to print bond distances up to a maximum of 'printdist' input units */
298 {
299 int i, j, k, l, nvert1, *itype;
300
301 float *d, dmin;
302
303 double occ;
304
305 if (!(itype = (int *) zalloc ((unsigned) (drvui->verts_alloc * sizeof (int))))) {
306 Error_Box ("Unable to assign space for bond analysis.");
307 return;
308 }
309 fprintf (drvui->flout, "\n\n Bond Distance Analysis\n");
310 /* Find all atoms of this type in master list */
311 l = nvert = nvert1 = 1;
312 for (j = 0; j < natom; ++j) {
313 if (drvui->atoms[j].atom_fn != drvui->frame_no)
314 continue;
315 l = nvert;
316 find_all_in_box (j);
317 for (k = l; k < nvert; ++k)
318 itype[k] = j; // itype[k] has atom number j
319 }
320 for (i = 0; i < natom; ++i) { /* loop through atom types */
321 if (drvui->atoms[i].atom_fn != drvui->frame_no)
322 continue;
323 /* add the base entry to the beginning of the atom list - picked from vertex list
324 because position originating from x,y,z symmetry operator may be absent due to
325 crenel exclusion */
326 k = -1;
327 for (l = 1; l < nvert; l++) { /* find an entry for atom i */
328 if (itype[l] == i) {
329 for (j = 0; j < 3; ++j) {
330 o_vert[j] = o_vert[3 * l + j];
331 o_vert_nm[j] = o_vert_nm[3 * l + j];
332 s_vert[j] = s_vert[3 * l + j];
333 }
334 k = l;
335 break;
336 }
337 }
338 if (k < 0) {
339 /* There is no entry for atom i in the list due to limits, etc. Add the base one. */
340 for (j = 0; j < 3; j++) {
341 s_vert[j] = 0;
342 o_vert[j] = drvui->atoms[i].atom_xyz[j];
343 o_vert_nm[j] = o_vert[j];
344 for (l = 0; l < 3; l++)
345 s_vert[j] +=
346 (float) drvui->b_mat[l][j] * (drvui->atoms[i].atom_xyz[l] -
347 origin[l]);
348 }
349 }
350 if ((d = (float *) zalloc ((unsigned) (nvert * sizeof (float)))) == NULL) {
351 free (itype);
352 Error_Box ("Unable to get bond distance allocation");
353 return;
354 }
355 for (j = 1; j < nvert; ++j) /* calculate distances */
356 d[j] = dist (0, j);
357 fprintf (drvui->flout,
358 "\nDistances from %c%c%c%c%3d at %8.5f%8.5f%8.5f to\n\n",
359 drvui->atoms[i].atom_l[0], drvui->atoms[i].atom_l[1],
360 drvui->atoms[i].atom_l[2], drvui->atoms[i].atom_l[3],
361 drvui->atoms[i].sv_atom_n, o_vert[0], o_vert[1], o_vert[2]);
362 for (k = nvert1; k < nvert; ++k) {
363 dmin = 999999.9f;
364 for (j = nvert1; j < nvert; ++j) { /* loop through distances */
365 if (d[j] < dmin) {
366 l = j;
367 dmin = d[j];
368 }
369 }
370 if (dmin == 0.0 && itype[l] == i) {
371 d[l] = 999999.9f;
372 } else {
373 if (dmin <= printdist) {
374 fprintf (drvui->flout, " %c%c%c%c%3d at %9.5f%9.5f%9.5f %8.3f",
375 drvui->atoms[itype[l]].atom_l[0],
376 drvui->atoms[itype[l]].atom_l[1],
377 drvui->atoms[itype[l]].atom_l[2],
378 drvui->atoms[itype[l]].atom_l[3],
379 drvui->atoms[itype[l]].sv_atom_n, o_vert[3 * l],
380 o_vert[3 * l + 1], o_vert[3 * l + 2], dmin);
381 d[l] = 999999.9f;
382 if (drvui->atoms[i].atom_ismod != 0 || drvui->atoms[itype[l]].atom_ismod != 0) { /* one or both ends are modulated */
383 float temp_pa[3];
384
385 float vert1[3], vert2[3], vert3[3];
386
387 float dsum = 0.0f, dmax = 0.0f, thisd;
388
389 int m, mx[3], m1, m2, number = 0;
390
391 dmin = 999999.0f;
392 for (m = 0; m < 3; m++) {
393 temp_pa[m] = drvui->phaseshift[m];
394 mx[m] = 50;
395 if (drvui->modulated <= m)
396 mx[m] = 1;
397 }
398 for (m = 0; m < mx[0]; m++) { /* sum over values of initial phaseshift for first modulation */
399 drvui->phaseshift[0] = m * 0.02f;
400 for (m1 = 0; m1 < mx[1]; m1++) {
401 drvui->phaseshift[1] = m1 * 0.02f; /* second modulation */
402 for (m2 = 0; m2 < mx[2]; m2++) {
403 drvui->phaseshift[2] = m2 * 0.02f; /* third modulation */
404 for (j = 0; j < 3; j++) {
405 vert1[j] = o_vert_nm[j];
406 vert2[j] = o_vert_nm[3 * l + j];
407 }
408 modulate_parameters (vert1, &occ, 0, i);
409 modulate_parameters (vert2, &occ, vert_sym_no[l],
410 itype[l]);
411 for (j = 0; j < 3; ++j) { /* convert vertex coordinates to Cartesian */
412 vert3[j] = 0.0f;
413 for (k = 0; k <= 2; ++k)
414 vert3[j] +=
415 (float) drvui->b_mat[j][k] * (vert2[k] -
416 vert1[k]);
417 }
418 thisd =
419 (float) sqrt (vert3[0] * vert3[0] +
420 vert3[1] * vert3[1] +
421 vert3[2] * vert3[2]);
422 dsum += thisd;
423 dmin = min (thisd, dmin);
424 dmax = max (thisd, dmax);
425 number++;
426 }
427 }
428 }
429 dsum /= (float) number; /* get mean distance */
430 for (m = 0; m < 3; m++)
431 drvui->phaseshift[m] = temp_pa[m];
432 fprintf (drvui->flout, "%8.3f %8.3f %8.3f", dsum, dmin, dmax);
433 }
434 fprintf (drvui->flout, "\n");
435 } else
436 break;
437 }
438 }
439 free (d);
440 }
441 free (itype);
442 fprintf (drvui->flout, "\n\n");
443 }
444
445 /* ************************************************************** */
446 /* ************************************************************** */
447
448 void
axeqb(double a1[3][3],double x[3],double b[3])449 axeqb (double a1[3][3], double x[3], double b[3])
450
451 /* solve matrix equation ax = b for x, where a is 3x3 and x and b are
452 vectors */
453 {
454 int i, j;
455
456 double a[3][3], det;
457
458 det = determinant (a1);
459 if (fabs (det) < 1.0e-6)
460 det = 1.0e-6;
461 a[0][0] = (a1[1][1] * a1[2][2] - a1[1][2] * a1[2][1]) / det;
462 a[1][0] = -(a1[1][0] * a1[2][2] - a1[1][2] * a1[2][0]) / det;
463 a[2][0] = (a1[1][0] * a1[2][1] - a1[1][1] * a1[2][0]) / det;
464 a[0][1] = -(a1[0][1] * a1[2][2] - a1[0][2] * a1[2][1]) / det;
465 a[1][1] = (a1[0][0] * a1[2][2] - a1[0][2] * a1[2][0]) / det;
466 a[2][1] = -(a1[0][0] * a1[2][1] - a1[0][1] * a1[2][0]) / det;
467 a[0][2] = (a1[0][1] * a1[1][2] - a1[0][2] * a1[1][1]) / det;
468 a[1][2] = -(a1[0][0] * a1[1][2] - a1[0][2] * a1[1][0]) / det;
469 a[2][2] = (a1[0][0] * a1[1][1] - a1[0][1] * a1[1][0]) / det;
470 for (i = 0; i < 3; i++) {
471 x[i] = 0.0;
472 for (j = 0; j < 3; j++)
473 x[i] += a[i][j] * b[j];
474 }
475 }
476
477 /* ************************************************************** */
478 /* ************************************************************** */
479
480 void
build_box_contents(void)481 build_box_contents (void)
482
483 /* routine to locate all atoms in the display parallelopiped plus
484 any atoms needed to complete molecules or polyhedra */
485 {
486 float vert[3]; /* storage for center */
487
488 int ix[3]; /* array for lattice offset */
489
490 int j, k, saved_nvert;
491
492 int atoms, loop;
493
494 loop = 0;
495 while (1) { // infinite loop here
496 saved_nvert = NvertM; // save number currently in list
497 for (atoms = 0; atoms < natom; atoms++) { // loop through atoms
498 if (drvui->atoms[atoms].atom_fn != drvui->frame_no)
499 continue;
500 for (j = 0; j < 3; j++) { // put the base atom in the range 0-1
501 drvui->atoms[atoms].saved_xyz[j] = drvui->atoms[atoms].atom_xyz[j];
502 if (drvui->atoms[atoms].atom_xyz[j] < 0.0f)
503 drvui->atoms[atoms].atom_xyz[j] += 1.0f;
504 if (drvui->atoms[atoms].atom_xyz[j] > 1.0f)
505 drvui->atoms[atoms].atom_xyz[j] -= 1.0f;
506 }
507 expand_atom (atoms); // generate all atoms of this flavor in unit cell
508 for (j = 0; j < ncell; ++j) { // atoms in cell
509 for (k = 0; k <= 2; ++k)
510 vert[k] = drvui->cell_xyz[j][k];
511 add_vert (vert, 1000 * atoms + drvui->sym_op_no[j], 0, 1, drvui->sym_op_signed[j]); // add point to cell if inside box
512 }
513 }
514 if (saved_nvert == NvertM)
515 break;
516 }
517 while (1) { /* infinite loop here */
518 saved_nvert = NvertM; /* save number currently in list */
519 for (atoms = 0; atoms < natom; atoms++) { /* loop through atoms */
520 if (drvui->atoms[atoms].atom_fn != drvui->frame_no)
521 continue;
522 expand_atom (atoms); /* generate all atoms of this flavor in unit cell */
523 for (j = 0; j < ncell; ++j) { /*atoms in cell */
524 for (k = 0; k <= 2; ++k)
525 vert[k] = drvui->cell_xyz[j][k];
526 add_vert (vert, 1000 * atoms + drvui->sym_op_no[j], 0, 1, drvui->sym_op_signed[j]); /* add point to cell if inside box */
527 }
528 for (ix[0] = (int) drvui->frames[drvui->frame_no].cryst_lim[0] - 1; ix[0] <= (int) drvui->frames[drvui->frame_no].cryst_lim[3] + 1; ++ix[0]) { /* step through cells parallel a */
529 for (ix[1] = (int) drvui->frames[drvui->frame_no].cryst_lim[1] - 1; ix[1] <= (int) drvui->frames[drvui->frame_no].cryst_lim[4] + 1; ++ix[1]) { /* parallel b */
530 for (ix[2] = (int) drvui->frames[drvui->frame_no].cryst_lim[2] - 1; ix[2] <= (int) drvui->frames[drvui->frame_no].cryst_lim[5] + 1; ++ix[2]) { /* parallel c */
531 if (abs (ix[0]) + abs (ix[1]) + abs (ix[2]) != 0) {
532 for (j = 0; j < ncell; ++j) { /*atoms in cell */
533 for (k = 0; k <= 2; ++k)
534 vert[k] = ix[k] + drvui->cell_xyz[j][k];
535 add_vert (vert, 1000 * atoms + drvui->sym_op_no[j], 0, 1, drvui->sym_op_signed[j]); /* add point to cell if inside box */
536 }
537 }
538 }
539 }
540 }
541 }
542 if ((domolcomp == 0 && drvui->npoly == 1) || saved_nvert == NvertM) {
543 for (atoms = 0; atoms < natom; atoms++) { /* loop through atoms */
544 if (drvui->atoms[atoms].atom_fn != drvui->frame_no)
545 continue;
546 for (j = 0; j < 3; j++) { // restore base atom position
547 drvui->atoms[atoms].atom_xyz[j] = drvui->atoms[atoms].saved_xyz[j];
548 }
549 }
550 return;
551 }
552 if (loop == 0)
553 fprintf (drvui->fcns, "\n");
554 loop++;
555 fprintf (drvui->fcns,
556 "Recursive Atom Addition, Loop No. %3d, %4d Atoms in List.\n", loop,
557 NvertM - 1);
558 }
559 } /* end of build_box_contents */
560
561
562 /* ************************************************************** */
563 /* ************************************************************** */
564
565 void
Calc_Rot(float v1[3],float v2[3])566 Calc_Rot (float v1[3], float v2[3])
567
568 /* routine to calculate rotation axes needed to put v1 parallel to Z
569 and the projection of v2 parallel to Y */
570 {
571 float u1[3], u2[3], u3[3], temp;
572
573 double U[3][3];
574
575 int i, j;
576
577 for (i = 0; i < 3; i++)
578 for (j = 0; j < 3; j++)
579 U[i][j] = 0.0f; /* clear U */
580 temp = 0.0f;
581 for (i = 0; i < 3; i++) { /* convert v1 to Cartesian */
582 u1[i] = 0.0f;
583 for (j = 0; j < 3; j++)
584 u1[i] += (float) drvui->b_mat[i][j] * v1[j];
585 temp += u1[i] * u1[i];
586 }
587 if (temp <= 0.0) {
588 u1[2] = 1.0f; /* if null input, make u1 = 0,0,1 */
589 U[2][2] = 1.0f;
590 } else {
591 temp = (float) sqrt (temp);
592 for (i = 0; i < 3; i++) { /* make a unit vector */
593 u1[i] /= temp;
594 U[i][2] = u1[i]; /* put column in U */
595 }
596 }
597 for (i = 0; i < 3; i++) { /* convert v2 to Cartesian */
598 u2[i] = 0.0f;
599 for (j = 0; j < 3; j++)
600 u2[i] += (float) drvui->b_mat[i][j] * v2[j];
601 }
602 u3[0] = u2[1] * u1[2] - u1[1] * u2[2]; /* calculate u2 X u1 */
603 u3[1] = u2[2] * u1[0] - u2[0] * u1[2];
604 u3[2] = u2[0] * u1[1] - u1[0] * u2[1];
605 temp = 0.0f;
606 for (i = 0; i < 3; i++)
607 temp += u3[i] * u3[i];
608 if (temp <= 0.0) {
609 u3[0] = 1.0f; /* if null input, make u3 = 1,0,0 */
610 U[0][0] = 1.0f;
611 } else {
612 temp = (float) sqrt (temp);
613 for (i = 0; i < 3; i++) { /* make a unit vector */
614 u3[i] /= temp;
615 U[i][0] = u3[i]; /* add column to U */
616 }
617 }
618 u2[0] = u1[1] * u3[2] - u3[1] * u1[2]; /* calculate u2 = u1 X u3 */
619 u2[1] = u1[2] * u3[0] - u1[0] * u3[2];
620 u2[2] = u1[0] * u3[1] - u3[0] * u1[1];
621 for (i = 0; i < 3; i++)
622 U[i][1] = u2[i];
623 /* make determinant of U > 0 */
624 if (determinant (U) < 0.0) {
625 for (i = 0; i < 3; i++)
626 U[i][0] = -U[i][0];
627 }
628
629 /* The rotation required to get the desired orientation, which is described
630 by G_Rot, has been set up so G_Rot * U = I. Thus G_Rot = U(inv),
631 U is unitary and U(inv) = U(transpose). */
632
633 for (i = 0; i < 3; i++)
634 for (j = 0; j < 3; j++)
635 G_Rot[i][j] = (float) U[j][i];
636 /* calculate rotations */
637 temp = (float) (G_Rot[1][2] * G_Rot[1][2] + G_Rot[2][2] * G_Rot[2][2]);
638 if (temp > 0.001) {
639 float temp2;
640
641 Yrot = (float) (atan (G_Rot[0][2] / sqrt (temp)) * RAD);
642 temp2 = (float) (G_Rot[1][2] / cos (Yrot / RAD));
643 if (temp2 > 1.0)
644 temp2 = 1.0f;
645 if (temp2 < -1.0)
646 temp2 = -1.0f;
647 Xrot = -(float) asin (temp2) * (float) RAD;
648 temp2 = -(float) (G_Rot[0][1] / cos (Yrot / RAD));
649 if (temp2 > 1.0)
650 temp2 = 1.0f;
651 if (temp2 < -1.0)
652 temp2 = -1.0f;
653 Zrot = (float) asin (temp2) * (float) RAD;
654 } else {
655 Yrot = -(float) asin (G_Rot[0][2]) * (float) RAD;
656 Xrot = 0.0f; /* XXXXXXXXXXX - May not be right */
657 Zrot = 0.0f; /* XXXXXXXXXXX - May not be right */
658 }
659 }
660
661
662 /* ************************************************************** */
663 /* ************************************************************** */
664
665 void
convert_ellipsoid(void)666 convert_ellipsoid (void)
667
668 /* routine to convert ellipsoid coefficients Bij or Uij to betaij form */
669 {
670 int i, j, save;
671
672 for (i = 1; i < drvui->n_ellips; ++i) {
673 if (drvui->ellips[i].ellips_ismod != 0)
674 continue;
675 save = drvui->ellips[i].ell_type / 1000;
676 if (drvui->ellips[i].ell_type % 1000 != 0) { // multiply Uij and Bij by a*(i)*a*(j)/4
677 for (j = 0; j < 3; ++j)
678 drvui->ellips[i].ellips[j] *=
679 drvui->rec_lat_con[j] * drvui->rec_lat_con[j] * 0.25f;
680 drvui->ellips[i].ellips[3] *=
681 drvui->rec_lat_con[0] * drvui->rec_lat_con[1] * 0.25f;
682 drvui->ellips[i].ellips[4] *=
683 drvui->rec_lat_con[0] * drvui->rec_lat_con[2] * 0.25f;
684 drvui->ellips[i].ellips[5] *=
685 drvui->rec_lat_con[1] * drvui->rec_lat_con[2] * 0.25f;
686 if (drvui->ellips[i].ell_type % 1000 == 1) // multiply Uij by 8pi^2
687 for (j = 0; j < 6; ++j)
688 drvui->ellips[i].ellips[j] *= 78.9568f;
689 }
690 drvui->ellips[i].ell_type = 1000 * save;
691 }
692 } /* end of convert_ellipsoid */
693
694 /* ************************************************************** */
695 /* ************************************************************** */
696
697 void
Conv_Sym_Mat(void)698 Conv_Sym_Mat (void)
699
700 /* routine to convert the rotational part of a symmetry matrix to a Cartesian rotation */
701 {
702 int i, j, k, l;
703
704 double mat[3][3];
705
706 float B[3][3], Binv[3][3];
707
708 float tmp[3][3];
709
710 if (drvui->sys == 5) { /* hexagonal systems are special */
711 for (i = 0; i < 3; i++)
712 for (j = 0; j < 3; j++) {
713 B[i][j] = (float) drvui->b_mat[i][j];
714 Binv[i][j] = B[i][j];
715 }
716 matinv (Binv);
717 }
718 for (l = 0; l < drvui->ng; l++) {
719 for (j = 0; j < 3; j++) {
720 for (k = 0; k < 3; k++) {
721 mat[j][k] = drvui->rss[l][j][k];
722 }
723 }
724 if (determinant (mat) < 0.0) { /* if improper rotation, negate all elements */
725 for (j = 0; j < 3; j++) {
726 for (k = 0; k < 3; k++) {
727 mat[j][k] *= -1.0f;
728 }
729 }
730 }
731 for (j = 0; j < 3; j++)
732 for (k = 0; k < 3; k++)
733 drvui->rssC[l][j][k] = (float) mat[j][k]; /* save a copy of the original */
734
735 if (drvui->sys == 5) { /* hexagonal systems are special */
736 float R[3][3];
737
738 for (j = 0; j < 3; j++)
739 for (k = 0; k < 3; k++)
740 R[j][k] = (float) mat[j][k];
741 matmul (B, R, tmp);
742 matmul (tmp, Binv, R); /* The Cartesian matrix is B R Binv */
743 for (j = 0; j < 3; j++)
744 for (k = 0; k < 3; k++)
745 mat[j][k] = R[j][k];
746 }
747 for (j = 0; j < 3; j++) { /* restore matrix */
748 for (k = 0; k < 3; k++) {
749 drvui->rss[l][j][k] = (float) mat[j][k];
750 }
751 }
752 }
753 }
754
755 /* ************************************************************** */
756 /* ************************************************************** */
757
758 float
dist(int j,int k)759 dist (int j, int k)
760
761 /* calculates distance between vertices at s_vert[j] and s_vert[k] */
762 {
763 float d, t;
764
765 int i;
766
767 d = 0.0f;
768 for (i = 0; i <= 2; ++i) {
769 t = s_vert[3 * j + i] - s_vert[3 * k + i];
770 d += t * t;
771 }
772 return ((float) sqrt (d));
773 }
774
775 /* ************************************************************** */
776 /* ************************************************************** */
777
778 float
dot0_3d(float x0,float y0,float z0,float x1,float y1,float z1,float x2,float y2,float z2)779 dot0_3d (float x0, float y0, float z0, float x1, float y1, float z1,
780 float x2, float y2, float z2)
781 {
782 // 3D dot product of (P1-P0) and (P2-P0), John Burkardt
783
784 return (x1 - x0) * (x2 - x0) + (y1 - y0) * (y2 - y0) + (z1 - z0) * (z2 - z0);
785 }
786
787 /* ************************************************************** */
788 /* ************************************************************** */
789
790 int
eigen(float * biso,float beta[6],float valu[3],float vect[3][3])791 eigen (float *biso, float beta[6], float valu[3], float vect[3][3])
792
793 /* routine to calculate eigenvalues 'valu' and vectors 'v' for
794 * array of anisotropic thermal parameters 'beta', 'biso' is the equivalent
795 * isotropic B
796 *
797 * some parts of this routine were copied from ORTEP
798 */
799 {
800 static double errnd = 1.0E-9;
801
802 double b0, b1, b2;
803
804 double a1, a2, phi;
805
806 double a[3][3], b[3][3], w[3][3], u[3];
807
808 float v1[3], v2[3], v3[3];
809
810 double tem, sigma, smax;
811
812 double Vect[3][3];
813
814 int i, j, l, ii, iii, i1, imax = 0;
815
816 int neq;
817
818 static double tpi2 = 2.0 * M_PI * M_PI; /* 2 pi^2 */
819
820 /* put beta into square symmetric matrix b (called M in ORTEP) */
821
822 for (i = 0; i < 3; i++)
823 b[i][i] = beta[i];
824 b[0][1] = b[1][0] = beta[3];
825 b[0][2] = b[2][0] = beta[4];
826 b[1][2] = b[2][1] = beta[5];
827 /* multiply b * ginv to get w */
828 for (i = 0; i < 3; i++) {
829 for (j = 0; j < 3; j++) {
830 w[i][j] = 0.0;
831 for (l = 0; l < 3; l++) {
832 w[i][j] += b[l][i] * drvui->ginv[j][l];
833 }
834 }
835 }
836 sigma = 0.0;
837 for (i = 0; i < 3; i++) {
838 for (j = 0; j < 3; j++) {
839 a[i][j] = w[i][j];
840 sigma += a[i][j] * a[i][j];
841 }
842 }
843 if (sigma <= 0.0)
844 return (0); /* error on null 'w' */
845 sigma = sqrt (sigma);
846 /* get coefficients of third-order characteristic equation
847 * Equation is -L^3 + b2 L^2 + b1 L + b0 = 0
848 */
849 b2 = a[0][0] + a[1][1] + a[2][2];
850 b1 = -a[0][0] * a[1][1] - a[0][0] * a[2][2] - a[1][1] * a[2][2]
851 + a[0][2] * a[2][0] + a[0][1] * a[1][0] + a[1][2] * a[2][1];
852 b0 = a[0][0] * a[1][1] * a[2][2] + a[1][0] * a[2][1] * a[0][2]
853 + a[2][0] * a[1][2] * a[0][1] - a[2][0] * a[0][2] * a[1][1]
854 - a[0][0] * a[2][1] * a[1][2] - a[1][0] * a[0][1] * a[2][2];
855
856 *biso = (float) (4.0 * b2 / 3.0);
857 /* transform to get to form x^3 + a1 x + a2 == 0 */
858 a1 = (-3.0 * b1 - b2 * b2) / 3.0;
859 a2 = (-2.0 * b2 * b2 * b2 - 9.0 * b2 * b1 - 27.0 * b0) / 27.0;
860 tem = (a2 * a2) / 4.0 + (a1 * a1 * a1) / 27.0;
861 if (tem > errnd) {
862 /* two complex roots */
863 fprintf (drvui->flout, "Complex Roots! tem = %f\n", tem);
864 return (0); /* error if complex roots */
865 } else if (tem < -errnd) {
866 /* 3 real and unequal roots */
867 phi = acos (-(a2 / 2.0) / sqrt (-(a1 * a1 * a1) / 27.0));
868 u[0] = 2.0 * sqrt (-a1 / 3.0) * cos (phi / 3.0) + b2 / 3.0;
869 u[1] = 2.0 * sqrt (-a1 / 3.0) * cos (phi / 3.0 + 120.0 / RAD) + b2 / 3.0;
870 u[2] = 2.0 * sqrt (-a1 / 3.0) * cos (phi / 3.0 + 240.0 / RAD) + b2 / 3.0;
871 } else {
872 double sign_a2 = 1.0;
873 /* Workaround brain-dead cbrt routine in Windows */
874 if (a2 < 0.0) {
875 sign_a2 = -1.0;
876 a2 = fabs(a2);
877 }
878 /* 3 real roots with at least 2 being equal */
879 u[0] = -2.0 * sign_a2 * cbrt (a2 / 2.0) + b2 / 3.0;
880 u[1] = sign_a2 * cbrt (a2 / 2.0) + b2 / 3.0;
881 u[2] = u[1];
882 }
883 for (j = 0; j < 2; j++) { /* sort roots in increasing order */
884 if (u[j] > u[j + 1]) {
885 tem = u[j + 1];
886 u[j + 1] = u[j];
887 u[j] = tem;
888 }
889 }
890 for (i = 0; i < 3; i++)
891 if (valu[i] < 0.0)
892 valu[i] = 0.001f;
893 else
894 valu[i] = (float) sqrt (u[i] / tpi2);
895
896 /* count multiple roots */
897 neq = 1;
898 for (j = 0; j < 2; j++)
899 if (fabs (u[j] - u[j + 1]) < errnd)
900 neq++;
901 if (neq > 2) {
902 /* 3 equal roots */
903 for (ii = 0; ii < 3; ii++) {
904 for (i = 0; i < 3; i++)
905 vect[i][ii] = 0.0;
906 vect[ii][ii] = 1.0;
907 }
908 return (1);
909 } else if (neq == 2) {
910 for (ii = 0; ii < 3; ii++) {
911 v1[ii] = v2[ii] = v3[ii] = 0.0;
912 }
913 if (drvui->sys == 5) { // hexagonal
914 /* With 2 equal eigenvalues, we have degenerate vectors */
915 if (fabs (u[0] - u[1]) < errnd) {
916 v3[2] = 1.0f; // 3rd vector is unique
917 v1[0] = 1.0f;
918 } else {
919 v3[0] = 1.0f; // 1st one is
920 v1[2] = 1.0f;
921 }
922 } else if (drvui->sys == 6) { // cubic with atom at xxx
923 if (fabs (u[0] - u[1]) < errnd) {
924 v3[0] = v3[1] = v3[2] = (float) sqrt (3.0) / 3.0f;
925 v1[0] = -(float) sqrt (2.0) / 2.0f;
926 v1[1] = -v1[0];
927 } else {
928 v1[0] = v1[1] = v1[2] = (float) sqrt (3.0) / 3.0f;
929 v3[0] = -(float) sqrt (2.0) / 2.0f;
930 v3[1] = -v3[0];
931 }
932 } else { // tetragonal
933 if (fabs (u[0] - u[1]) < errnd) {
934 v3[2] = 1.0f; // 3rd vector is unique
935 v2[1] = 1.0f;
936 v1[0] = 1.0f;
937 } else {
938 v3[0] = 1.0f; // 1st one is
939 v1[1] = 1.0f;
940 v2[2] = 1.0f;
941 }
942 }
943 vcross (v3, v1, v2); // v1, v2 and v3 are mutually perp. unit vectors
944 vnormalize (v2);
945 for (j = 0; j < 3; j++) {
946 vect[0][j] = v1[j];
947 vect[1][j] = v2[j];
948 vect[2][j] = v3[j];
949 }
950 return (1);
951 }
952 for (ii = 0; ii < 3; ii++) {
953 for (iii = 0; iii < 3; iii++) {
954 for (i = 0; i < 3; i++)
955 a[i][iii] = w[i][iii];
956 a[iii][iii] = w[iii][iii] - u[ii];
957 }
958 smax = 0.0;
959 for (i = 0; i < 3; i++) {
960 i1 = (i < 2) ? i + 1 : 0;
961 b[0][i] = a[1][i] * a[2][i1] - a[2][i] * a[1][i1];
962 b[1][i] = a[2][i] * a[0][i1] - a[0][i] * a[2][i1];
963 b[2][i] = a[0][i] * a[1][i1] - a[1][i] * a[0][i1];
964 tem = b[0][i] * b[0][i] + b[1][i] * b[1][i] + b[2][i] * b[2][i];
965 if (tem > smax) {
966 smax = tem;
967 imax = i;
968 }
969 }
970 if (smax <= errnd)
971 return (0);
972 smax = (float) sqrt (smax);
973 for (i = 0; i < 3; i++)
974 v2[i] = (float) (b[i][imax] / smax);
975 for (i = 0; i < 3; i++) { /* Convert v2 to cartesian system */
976 v1[i] = 0.0;
977 for (j = 0; j < 3; j++)
978 v1[i] += (float) drvui->b_mat[i][j] * v2[j];
979 }
980 vnormalize (v1);
981 for (i = 0; i < 3; i++)
982 Vect[ii][i] = v1[i]; /* unit vector */
983 }
984 if (determinant (Vect) < 0.97) {
985 for (i = 0; i < 3; i++) {
986 v2[i] = (float) Vect[1][i];
987 v3[i] = (float) Vect[2][i];
988 }
989 vcross (v2, v3, v1); /* Set is not orthogonal -the smallest root is probably wrong */
990 for (i = 0; i < 3; i++)
991 Vect[0][i] = v1[i];
992 }
993 for (i = 0; i < 3; i++)
994 for (j = 0; j < 3; j++)
995 vect[i][j] = (float) Vect[i][j];
996 return (1);
997 }
998
999 /* ************************************************************** */
1000 /* ************************************************************** */
1001
1002 void
find_all_in_box(int i)1003 find_all_in_box (int i)
1004
1005 /* routine to locate all atoms of type 'i' in the display box */
1006 /* i - atom number */
1007 {
1008 int n, j, k, l; /* loop indices */
1009
1010 int sym_no; /* symmetry number */
1011
1012 #ifndef WIN32
1013 int* itmp;
1014
1015 float* ftmp;
1016 #endif
1017
1018 for (n = 1; n < NvertM; n++) { /* loop through the Master list */
1019 j = (int) drvui->atom_no[n] / 1000; /* extract atom number */
1020 if (j == i) {
1021 drvui->orig_atom_no[nvert] = i;
1022 sym_no = drvui->atom_no[n] - 1000 * j; /* symmetry operator number */
1023 for (k = 0; k < 3; k++) {
1024 o_vert[3 * nvert + k] = xypos[3 * n + k];
1025 o_vert_nm[3 * nvert + k] = xypos_nm[3 * n + k];
1026 }
1027 for (k = 0; k < 3; k++) {
1028 s_vert[3 * nvert + k] = 0.0f;
1029 for (l = 0; l < 3; l++) { /* calculate cartesian coordinates */
1030 s_vert[3 * nvert + k] +=
1031 (float) drvui->b_mat[k][l] * (o_vert[3 * nvert + l] - origin[l]);
1032 }
1033 }
1034 vert_sym_nos[nvert] = drvui->atom_so[n];
1035 vert_sym_no[nvert++] = sym_no;
1036 if (nvert > drvui->verts_alloc) {
1037 Error_Box
1038 ("Too many vertices for dimensions - please report to developers.");
1039 return;
1040 }
1041 if (nvert > 2 * NvertM) {
1042 #ifdef WIN32
1043 Error_Box ("Overrun of vert_sym_no. Please send 'str' file\n"
1044 "to Larry.Finger@@lwfinger.net.");
1045 exit (0);
1046 #else
1047 itmp = (int *) realloc (vert_sym_no, (nvert + 2) * sizeof (int));
1048 if (itmp)
1049 vert_sym_no = itmp;
1050 else
1051 Error_Box ("Unable to expand storage for vert_sym_no. (Out of memory)\n");
1052 itmp = (int *) realloc (vert_sym_nos, (nvert + 2) * sizeof (int));
1053 if (itmp)
1054 vert_sym_nos = itmp;
1055 else
1056 Error_Box ("Unable to expand storage for vert_sym_nos. (Out of memory)\n");
1057 itmp = (int *) realloc (drvui->orig_atom_no, (nvert + 2) * sizeof (int));
1058 if (itmp)
1059 drvui->orig_atom_no = itmp;
1060 else
1061 Error_Box ("Unable to expand storage for orig_atom_no. (Out of memory)\n");
1062 ftmp = (float *) realloc (o_vert, (6 * (nvert + 2) * sizeof (float)));
1063 if (ftmp)
1064 o_vert = ftmp;
1065 else
1066 Error_Box ("Unable to expand storage for o_vert. (Out of memory)\n");
1067 ftmp = (float *) realloc (o_vert_nm, (6 * (nvert + 2) * sizeof (float)));
1068 if (ftmp)
1069 o_vert_nm = ftmp;
1070 else
1071 Error_Box ("Unable to expand storage for o_vert_nm. (Out of memory)\n");
1072 #endif
1073 }
1074 }
1075 }
1076 } /* end of find_all_in_box */
1077
1078 /* ************************************************************** */
1079 /* ************************************************************** */
1080
1081 void
generate_arrows(void)1082 generate_arrows (void)
1083 // routine to generate arrow objects for magnetic moments
1084 {
1085 int Arrow_Count; // Counter for number output
1086
1087 int i, k, l; // loop counters
1088
1089 float glr, glb, glg;
1090
1091 GLUquadricObj *glu_quadObj;
1092
1093 float base;
1094
1095 float x[3], xp[3], yp[3];
1096
1097 int xi, xj, xk;
1098
1099 float phi, chi;
1100
1101 int savedNvertM = NvertM; // save the end of the vertex list
1102
1103 char col_arrow_v[40];
1104
1105 char col_arrow_p[40];
1106
1107 Arrow_Count = 0;
1108
1109 glu_quadObj = gluNewQuadric ();
1110
1111 for (i = 0; i < drvui->nmag; ++i) { // loop through magnetic arrows
1112 if (drvui->arrows[i].arrow_fn != drvui->frame_no)
1113 continue;
1114 for (k = 0; k < 3; k++) {
1115 x[k] = drvui->arrows[i].mag_xp[k] - drvui->xyzoff[k];
1116 }
1117 for (xi = -3; xi <= 3; xi++) {
1118 xp[0] = (float) xi;
1119 for (xj = -3; xj <= 3; xj++) {
1120 xp[1] = (float) xj;
1121 for (xk = -3; xk <= 3; xk++) {
1122 xp[2] = (float) xk; // xp is a magnetic cell translation
1123 for (k = 0; k < 3; k++) {
1124 yp[k] = x[k]; // yp will be position in nuclear cell
1125 for (l = 0; l < 3; l++) {
1126 yp[k] += drvui->mag_matrix[k][l] * xp[l];
1127 }
1128 }
1129 add_vert (yp, 1000 * (i + drvui->atom_alloc), 0, 0, 0); // add vector to vertex list - NOT modulated??
1130 }
1131 }
1132 }
1133 } // end of loop through magnetic arrows
1134 for (k = savedNvertM; k < NvertM; k++) {
1135 float tmp, tmp1;
1136
1137 i = drvui->atom_no[k] / 1000 - drvui->atom_alloc; // get number of magnetic item
1138 strcpy (col_arrow_v, drvui->arrows[i].col_arrow);
1139 strcpy (col_arrow_p, col_arrow_v);
1140 Transform_VRML_Color (col_arrow_v);
1141 Transform_POV_Color (col_arrow_p);
1142 for (l = 0; l < 3; l++) {
1143 x[l] = drvui->arrows[i].mag_xc[l];
1144 }
1145 double t = sqrt (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
1146
1147 if (t < 0.002)
1148 continue; // skip this arrow if 0 length
1149 for (l = 0; l < 3; l++) {
1150 x[l] /= (float) t;
1151 }
1152 t = sqrt (x[0] * x[0] + x[1] * x[1]);
1153 chi = (float) atan2 (t, (double) x[2]);
1154 if (t > 0.001)
1155 phi = (float) atan2 (x[0] / t, x[1] / t);
1156 else
1157 phi = 0.0f;
1158 glLoadName (200000 + i);
1159 glPushName (k);
1160 glPushMatrix ();
1161 (void) sscanf (col_arrow_v, "%f %f %f", &glr, &glg, &glb);
1162 glColor3f (glr, glg, glb);
1163 glTranslatef (s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1164 glRotatef (-phi * (float) RAD, 0.0, 0.0, 1.0f);
1165 glRotatef (-chi * (float) RAD, 1.0f, 0., 0.);
1166 glTranslatef (0.0, 0.0, -drvui->arrows[i].arrow_length / 2.0f);
1167 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
1168 gluCylinder (glu_quadObj, drvui->arrows[i].arrow_diam,
1169 drvui->arrows[i].arrow_diam, 0.9 * drvui->arrows[i].arrow_length, 10,
1170 1);
1171 glPopMatrix ();
1172 glPushMatrix ();
1173 glTranslatef (s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1174 glRotatef (-phi * (float) RAD, 0.0, 0.0, 1.0f);
1175 glRotatef (-chi * (float) RAD, 1.0f, 0., 0.);
1176 glTranslatef (0.0f, 0.0f, drvui->arrows[i].arrow_length * 0.27f);
1177
1178 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
1179 gluCylinder (glu_quadObj, drvui->arrows[i].arrow_diam * 2.0,
1180 0.01, 0.4 * drvui->arrows[i].arrow_length, 10, 1);
1181 glPopMatrix ();
1182 glPopName ();
1183
1184 // arrow in POV
1185
1186 if (doPOV) {
1187 tmp = 0.5f * drvui->arrows[i].arrow_length;
1188 base = 2.0f * drvui->arrows[i].arrow_diam;
1189 fprintf (drvui->fpoutp, " object{\n");
1190 fprintf (drvui->fpoutp, " union{\n");
1191 fprintf (drvui->fpoutp, " object{\n");
1192 fprintf (drvui->fpoutp,
1193 " cylinder{ < 0,0,%8.5f > , < 0,0,%8.5f >,%8.5f\n",
1194 -0.5f * drvui->arrows[i].arrow_length,
1195 0.27f * drvui->arrows[i].arrow_length, drvui->arrows[i].arrow_diam);
1196 fprintf (drvui->fpoutp, " texture{pigment{color %s}}\n", col_arrow_p);
1197 fprintf (drvui->fpoutp,
1198 " finish{phong %5.2f phong_size %5.2f}}\n }\n",
1199 drvui->Phong_Value, drvui->Phong_Size);
1200 tmp1 = 0.27f * drvui->arrows[i].arrow_length;
1201 tmp = 0.77f * drvui->arrows[i].arrow_length;
1202 fprintf (drvui->fpoutp, " object{\n");
1203 fprintf (drvui->fpoutp,
1204 " cone{<0, 0, %8.5f>, %8.5f, <0, 0, %8.5f>, 0 \n", tmp1, base,
1205 tmp);
1206 fprintf (drvui->fpoutp, " texture{pigment{color %s}}\n", col_arrow_p);
1207 fprintf (drvui->fpoutp, " finish{phong %5.2f phong_size %5.2f}}\n }\n",
1208 drvui->Phong_Value, drvui->Phong_Size);
1209 fprintf (drvui->fpoutp, " rotate <%8.5f,0,0>\n", -chi * RAD);
1210 fprintf (drvui->fpoutp, " rotate <0,0,%8.5f>\n", -phi * RAD);
1211 fprintf (drvui->fpoutp, " translate <%8.5f, %8.5f, %8.5f>\n }}\n",
1212 s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1213 }
1214 if (doAsy) {
1215 fprintf (drvui->fpouta, "draw(pic, (%5.2f,%5.2f,%5.2f)--(%5.2f,%5.2f,%5.2f),\n",
1216 s_vert[3 * k]-0.8f*drvui->arrows[i].mag_xc[0],
1217 s_vert[3 * k + 1]-0.8f*drvui->arrows[i].mag_xc[1],
1218 s_vert[3 * k + 2]-0.8f*drvui->arrows[i].mag_xc[2],
1219 s_vert[3 * k]+1.2f*drvui->arrows[i].mag_xc[0],
1220 s_vert[3 * k + 1]+1.2f*drvui->arrows[i].mag_xc[1],
1221 s_vert[3 * k + 2]+1.2f*drvui->arrows[i].mag_xc[2]);
1222 fprintf (drvui->fpouta, "\tlinewidth(%.2f)+rgb(%.2f,%.2f,%.2f),EndArrow3(DefaultHead3,%.2f),currentlight);\n",
1223 20.f*drvui->arrows[i].arrow_diam,glr,glg,glb,4.f*20.f*drvui->arrows[i].arrow_diam);
1224 }
1225 // arrow in VRML
1226 float d1, d2, d3, alpha, temp;
1227
1228 float cosalp;
1229
1230 d1 = x[1];
1231 d2 = -x[0];
1232 d3 = 0.0f;
1233 temp = (float) sqrt (d1 * d1 + d2 * d2 + d3 * d3);
1234 if (temp < 0.001f) {
1235 d1 = 1.0;
1236 temp = 1.0;
1237 }
1238 d1 /= temp;
1239 d2 /= temp;
1240 d3 /= temp;
1241 cosalp = x[2] / (float) sqrt (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);;
1242 alpha = -(float) acos (cosalp);
1243 if (doVrml) {
1244 if (Vrml2) {
1245 // cylinder for shaft of arrow
1246 fprintf (drvui->fpoutv, " Transform{\n");
1247 fprintf (drvui->fpoutv, " rotation %6.3f %6.3f %6.3f %7.3f\n",
1248 d1, d2, d3, alpha);
1249 fprintf (drvui->fpoutv, " translation %5.3f %5.3f %5.3f\n",
1250 s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1251 fprintf (drvui->fpoutv, " children [\n");
1252 fprintf (drvui->fpoutv, " Transform{\n");
1253 fprintf (drvui->fpoutv, " rotation 1 0 0 1.5708\n");
1254 fprintf (drvui->fpoutv, " translation 0 0 %5.3f\n",
1255 -0.05f * drvui->arrows[i].arrow_length);
1256 fprintf (drvui->fpoutv, " children [\n");
1257 fprintf (drvui->fpoutv, " Shape {\n");
1258 fprintf (drvui->fpoutv, " geometry Cylinder {\n");
1259 fprintf (drvui->fpoutv, " radius %f\n",
1260 drvui->arrows[i].arrow_diam);
1261 fprintf (drvui->fpoutv, " height %f\n",
1262 0.9f * drvui->arrows[i].arrow_length);
1263 fprintf (drvui->fpoutv, " }\n");
1264 fprintf (drvui->fpoutv,
1265 " appearance Appearance { material Material "
1266 "{ diffuseColor %s}}\n", col_arrow_v);
1267 fprintf (drvui->fpoutv, " }\n");
1268 fprintf (drvui->fpoutv, " ]\n");
1269 fprintf (drvui->fpoutv, " }\n");
1270 fprintf (drvui->fpoutv, " ]\n");
1271 fprintf (drvui->fpoutv, " }\n");
1272 // Cone at end of arrow
1273 fprintf (drvui->fpoutv, " Transform{\n");
1274 fprintf (drvui->fpoutv, " rotation %6.3f %6.3f %6.3f %7.3f\n",
1275 d1, d2, d3, alpha);
1276 fprintf (drvui->fpoutv, " translation %5.3f %5.3f %5.3f\n",
1277 s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1278 fprintf (drvui->fpoutv, " children Transform {\n");
1279 fprintf (drvui->fpoutv, " children [\n");
1280 fprintf (drvui->fpoutv, " Transform{\n");
1281 fprintf (drvui->fpoutv, " rotation 1 0 0 1.5708\n");
1282 fprintf (drvui->fpoutv, " translation 0 0 %5.3f\n",
1283 0.5f * drvui->arrows[i].arrow_length);
1284 fprintf (drvui->fpoutv, " children [\n");
1285 fprintf (drvui->fpoutv, " Shape {\n");
1286 fprintf (drvui->fpoutv, " geometry Cone {\n");
1287 fprintf (drvui->fpoutv, " bottomRadius %f\n",
1288 2.0 * drvui->arrows[i].arrow_diam);
1289 fprintf (drvui->fpoutv, " height %f\n",
1290 drvui->arrows[i].arrow_length * 0.4);
1291 fprintf (drvui->fpoutv, " }\n");
1292 fprintf (drvui->fpoutv,
1293 " appearance Appearance { material Material "
1294 "{ diffuseColor %s}}\n", col_arrow_v);
1295 fprintf (drvui->fpoutv, " }\n");
1296 fprintf (drvui->fpoutv, " ]\n");
1297 fprintf (drvui->fpoutv, " }\n");
1298 fprintf (drvui->fpoutv, " ]\n");
1299 fprintf (drvui->fpoutv, " }\n");
1300 fprintf (drvui->fpoutv, " }\n");
1301 } else { // VRML 1
1302 fprintf (drvui->fpoutv, " Separator{\n");
1303 fprintf (drvui->fpoutv, " Transform{\n");
1304 fprintf (drvui->fpoutv, " translation %5.3f %5.3f %5.3f\n",
1305 s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1306 fprintf (drvui->fpoutv, " rotation %6.3f %6.3f %6.3f %7.3f\n }\n",
1307 d1, d2, d3, alpha);
1308 fprintf (drvui->fpoutv, " Material {diffuseColor %s}\n", col_arrow_v);
1309 fprintf (drvui->fpoutv, " Rotation{rotation 1 0 0 1.5708}\n");
1310 fprintf (drvui->fpoutv, " Cylinder {\n");
1311 fprintf (drvui->fpoutv, " parts SIDES\n");
1312 fprintf (drvui->fpoutv, " radius %f\n", drvui->arrows[i].arrow_diam);
1313 fprintf (drvui->fpoutv, " height %f\n",
1314 0.9f * drvui->arrows[i].arrow_length);
1315 fprintf (drvui->fpoutv, " }\n }\n");
1316 fprintf (drvui->fpoutv, " Separator{\n");
1317 fprintf (drvui->fpoutv, " Transform{\n");
1318 fprintf (drvui->fpoutv, " translation %5.3f %5.3f %5.3f\n",
1319 s_vert[3 * k], s_vert[3 * k + 1], s_vert[3 * k + 2]);
1320 fprintf (drvui->fpoutv, " rotation %6.3f %6.3f %6.3f %7.3f\n }\n",
1321 d1, d2, d3, alpha);
1322 fprintf (drvui->fpoutv, " Material {diffuseColor %s}\n", col_arrow_v);
1323 fprintf (drvui->fpoutv, " Transform{\n");
1324 fprintf (drvui->fpoutv, " rotation 1 0 0 1.5708\n");
1325 fprintf (drvui->fpoutv, " translation 0 0 %5.3f\n }\n",
1326 0.5f * drvui->arrows[i].arrow_length);
1327 fprintf (drvui->fpoutv, " Cone {\n");
1328 fprintf (drvui->fpoutv, " parts SIDES\n");
1329 fprintf (drvui->fpoutv, " bottomRadius %f\n",
1330 2.0 * drvui->arrows[i].arrow_diam);
1331 fprintf (drvui->fpoutv, " height %f\n",
1332 drvui->arrows[i].arrow_length * 0.4);
1333 fprintf (drvui->fpoutv, " }\n }\n");
1334 }
1335 }
1336 Arrow_Count++;
1337 }
1338
1339 gluDeleteQuadric (glu_quadObj);
1340
1341 if (Arrow_Count > 0) {
1342 fprintf (drvui->fcns, "%4d magnetic moment arrows.\n", Arrow_Count);
1343 fprintf (drvui->flout, "Generated %4d arrows.\n", Arrow_Count);
1344 }
1345 NvertM = savedNvertM; // restore vertex count
1346 } // end of generate_arrows
1347
1348 /* ************************************************************** */
1349 /* ************************************************************** */
1350
1351 void
generate_bonds(void)1352 generate_bonds (void)
1353
1354 /* routine to generate bond descriptions and add then to edit list */
1355 {
1356 int start, same_type;
1357
1358 int i, j, k; /* loop counters */
1359
1360 int N_Bond; /* Counter for number of bonds of a type */
1361
1362 int l1; /* loop counter for matrix operations */
1363
1364 int nvert1; /* vertex counter */
1365
1366 int Bond_Count; /* Counter for bonds output */
1367
1368 char bnd1[4], bnd2[4]; /* temporary storage of bond labels */
1369
1370 float d; /* bond distance */
1371
1372 float df[3]; /* vector difference */
1373
1374 float beta; /* rotation angle */
1375
1376 float gamma; /* another rotation angle */
1377
1378 float at1[3], at2[3]; /*initial coordinates of bond endpoints */
1379
1380 float clip0, clip1; /* clipping scale factors */
1381
1382 int m, l2, l3, out = 0;
1383 int ell[2] = { 0, 0 }; /* save ellipsoid numbers */
1384
1385 int is_ellipsoid[2] = { 0, 0 }; /* ellipsoid markers - for moving bond endpoint */
1386 int o1, p, p1, q, r; /* loop variables for rotation matrix calculation */
1387
1388 float elrot0[3][3]; /* rotation matrices of ellipsoids */
1389
1390 float Z0[3], Z1[3]; /* components of rms displacement along bond */
1391
1392 int *ellips_id;
1393
1394 int *ellips_num;
1395
1396 int jj;
1397
1398 float factor, rms0, rms1;
1399
1400 float glr, glg, glb;
1401
1402 int numdashes;
1403
1404 float dashes = 5.0f;
1405
1406 float mf;
1407
1408 GLUquadricObj *glu_quadObj;
1409
1410 int omit;
1411
1412 char col_bond_p[40];
1413
1414 char col_bond_v[40];
1415
1416 if (drvui->nbond == 1)
1417 return; /* exit if no work to be done */
1418
1419 Bond_Count = 0;
1420 if (!
1421 (ellips_id =
1422 (int *) zalloc ((unsigned) ((long) (drvui->verts_alloc * sizeof (int)))))) {
1423 Error_Box ("Unable to get ellips_id allocation!");
1424 return;
1425 }
1426 if (!
1427 (ellips_num =
1428 (int *) zalloc ((unsigned) ((long) (drvui->verts_alloc * sizeof (int)))))) {
1429 Error_Box ("Unable to get ellips_num allocation!");
1430 free (ellips_id);
1431 return;
1432 }
1433
1434 glu_quadObj = gluNewQuadric ();
1435
1436 for (i = 1; i < drvui->nbond; ++i) { /* loop through bond types */
1437 if (drvui->bonds[i].bond_fn != drvui->frame_no)
1438 continue; // skip if not in this frame
1439 strcpy (col_bond_p, drvui->bonds[i].col_bond);
1440 strcpy (col_bond_v, col_bond_p);
1441 Transform_VRML_Color (col_bond_v);
1442 Transform_POV_Color (col_bond_p);
1443 nvert = 0; /* empty the vertex list */
1444 for (j = 0; j < natom; j++) { /* loop through atom types */
1445 if (drvui->atoms[j].atom_fn != drvui->frame_no)
1446 continue;
1447 int o = drvui->atoms[j].atom_n;
1448
1449 if (((o & 255) > 0) || (((o >> 24) & 255) > 0)) { // sphere or ellipsoid
1450
1451 o = (drvui->atoms[j].atom_n >> 24) & 255;
1452 if (o > 0) { // ellipsoid
1453 if (check_atom_name
1454 (drvui->bonds[i].bond_l1, drvui->ellips[o].ellips_l)) {
1455 if (drvui->do_ellipsoids == 1 && drvui->El_Cutout == 1 && drvui->ellips[o].ell_type > 100) { /* skip if B iso */
1456 is_ellipsoid[0] = 1; /* or no cutout */
1457 ell[0] = o;
1458 }
1459 o = nvert;
1460 find_all_in_box (j); // get all atoms of type j
1461 for (p = o; p < nvert; p++)
1462 ellips_id[p] = ell[0]; // save ellipsoid id for these entries
1463 ellips_num[p] = j;
1464 }
1465 } else { // sphere
1466 o = drvui->atoms[j].atom_n & 255; // number of sphere
1467 if (check_atom_name
1468 (drvui->bonds[i].bond_l1, drvui->spheres[o].sphere_l)) {
1469 find_all_in_box (j); /* get all atoms of type j */
1470 is_ellipsoid[0] = 0;
1471 }
1472 } /* o > 0 */
1473 }
1474 } /* for (j=0;j<natom;++j) */
1475
1476 bnd1[0] = drvui->bonds[i].bond_l1[0];
1477 bnd1[1] = drvui->bonds[i].bond_l1[1];
1478 bnd1[2] = drvui->bonds[i].bond_l1[2];
1479 bnd1[3] = drvui->bonds[i].bond_l1[3];
1480 nvert1 = nvert;
1481 if (check_atom_name (drvui->bonds[i].bond_l1, drvui->bonds[i].bond_l2)) {
1482 bnd2[0] = drvui->bonds[i].bond_l1[0]; /* atoms are same type */
1483 bnd2[1] = drvui->bonds[i].bond_l1[1];
1484 bnd2[2] = drvui->bonds[i].bond_l1[2];
1485 bnd2[3] = drvui->bonds[i].bond_l1[3];
1486 is_ellipsoid[1] = is_ellipsoid[0];
1487 ell[1] = ell[0];
1488 same_type = 1;
1489 } else {
1490 bnd2[0] = drvui->bonds[i].bond_l2[0];
1491 bnd2[1] = drvui->bonds[i].bond_l2[1];
1492 bnd2[2] = drvui->bonds[i].bond_l2[2];
1493 bnd2[3] = drvui->bonds[i].bond_l2[3];
1494 same_type = 0;
1495 /* atoms i and j are not the same type */
1496 for (j = 0; j < natom; ++j) { /* loop through atom types */
1497 if (drvui->atoms[j].atom_fn != drvui->frame_no)
1498 continue;
1499 if (((drvui->atoms[j].atom_n & 255) > 0) || (((drvui->atoms[j].atom_n >> 24) & 255) > 0)) { // sphere or ellipsoid
1500 int o;
1501
1502 o = (drvui->atoms[j].atom_n >> 24) & 255;
1503 if (o > 0) { // ellipsoid
1504 if (check_atom_name
1505 (drvui->bonds[i].bond_l2, drvui->ellips[o].ellips_l)) {
1506 if (drvui->do_ellipsoids == 1 && drvui->El_Cutout == 1 && drvui->ellips[o].ell_type > 100) { /* skip if B iso */
1507 is_ellipsoid[1] = 1; /* or no cutout */
1508 ell[1] = o;
1509 }
1510 o = nvert;
1511 find_all_in_box (j); /* get all atoms of type j */
1512 for (p = o; p < nvert; p++)
1513 ellips_id[p] = ell[1]; /* save ellipsoid type */
1514 ellips_num[p] = j;
1515 }
1516 } else {
1517 o = drvui->atoms[j].atom_n & 255; // sphere number
1518 if (check_atom_name
1519 (drvui->bonds[i].bond_l2, drvui->spheres[o].sphere_l)) {
1520 find_all_in_box (j); /* get all atoms of type j */
1521 is_ellipsoid[1] = 0;
1522 }
1523 } // o > 0
1524 }
1525 } /* for (j=0;j<natom;++j) */
1526 } /* if atoms are same type */
1527
1528 if (((nvert1 > 0) && (nvert > nvert1)) || (same_type == 1)) {
1529 N_Bond = 0;
1530 for (j = 0; j < nvert1; ++j) { /* loop through atoms of type 1 */
1531 start = nvert1;
1532 if (same_type == 1)
1533 start = j;
1534 for (k = start; k < nvert; ++k) { /* loop on type 2 atoms */
1535 d = dist (j, k); /* get atomic distance */
1536 if ((d >= drvui->bonds[i].bond_min) && (d <= drvui->bonds[i].bond_max)) { /* here if distance within range */
1537 df[0] = df[1] = drvui->Bond_Mult * drvui->bonds[i].bond_size; /* scaling */
1538 df[2] = 0.5f * d;
1539
1540 Bond_Count++; /* update count */
1541 for (p = 0; p < 3; p++) { /* initialize endpoints of bond */
1542 Z0[p] = s_vert[3 * j + p];
1543 Z1[p] = s_vert[3 * k + p];
1544 }
1545
1546 if (drvui->do_ellipsoids == 1 && drvui->El_Cutout == 1) {
1547
1548 /* Elipsoid has octant cutout, calculate rms displacement along bond */
1549
1550 /* calculate multiplier - the 0.9 forces the end of the bond closer to the ellipsoid */
1551 factor =
1552 (float) (drvui->Ellipsoid_Scale * 0.9 /
1553 (PI * d * 1.41421356));
1554
1555 /* Mean-square displacement of a scaled thermal ellipsoid along the direction of a bond is given by: */
1556 /* */
1557 /* D = factor^2 * (X'G'RBR'GX) */
1558 /* */
1559 /* where X is the interatomic vector given in triclinic fractional coordinates, */
1560 /* G is the lattice metric with the ij'th element given by a[i].a[j], */
1561 /* R is the rotational part of the symmetry operation that transform atom i, */
1562 /* B is a square matrix of the betaij thermal factors, */
1563 /* and a prime indicates the transpose of the matrix. */
1564 /* */
1565 if (is_ellipsoid[0] == 1) { /* at one end of bond */
1566 float B[3][3], BP[3][3];
1567
1568 o1 = vert_sym_no[j];
1569 p1 = ellips_id[j]; /* have original ellipsoid number */
1570 jj = ellips_num[j];
1571 if (drvui->ellips[p1].ellips_ismod != 0) { /* apply modulation */
1572 float uij[6];
1573
1574 float vert[3];
1575
1576 vert[0] = o_vert_nm[3 * j];
1577 vert[1] = o_vert_nm[3 * j + 1];
1578 vert[2] = o_vert_nm[3 * j + 2];
1579 if (modulate_uij (vert, p1, jj, o1, uij) == 1) { /* if n.p.d, keep average value */
1580 for (p = 0; p < 3; p++)
1581 B[p][p] = drvui->ellips[p1].ellips[p]; // extract betaij's into square matrix
1582 B[0][1] = B[1][0] = drvui->ellips[p1].ellips[3];
1583 B[0][2] = B[2][0] = drvui->ellips[p1].ellips[4];
1584 B[1][2] = B[2][1] = drvui->ellips[p1].ellips[5];
1585 } else { /* use modulated values */
1586 for (p = 0; p < 3; p++)
1587 B[p][p] = uij[p]; // extract betaij's into square matrix
1588 B[0][1] = B[1][0] = uij[3];
1589 B[0][2] = B[2][0] = uij[4];
1590 B[1][2] = B[2][1] = uij[5];
1591 }
1592 } else { /* not modulated */
1593 for (p = 0; p < 3; p++)
1594 B[p][p] = drvui->ellips[p1].ellips[p]; // extract betaij's into square matrix
1595 B[0][1] = B[1][0] = drvui->ellips[p1].ellips[3];
1596 B[0][2] = B[2][0] = drvui->ellips[p1].ellips[4];
1597 B[1][2] = B[2][1] = drvui->ellips[p1].ellips[5];
1598 }
1599 for (p = 0; p < 3; p++) {
1600 df[p] = 0.0;
1601 for (q = 0; q < 3; ++q) {
1602 elrot0[p][q] = 0.0;
1603 df[p] +=
1604 drvui->ginv[q][p] * (o_vert[3 * k + q] -
1605 o_vert[3 * j + q]);
1606 for (r = 0; r < 3; r++)
1607 elrot0[p][q] +=
1608 B[p][r] * drvui->rssC[o1][q][r];
1609 }
1610 }
1611 for (p = 0; p < 3; p++) {
1612 for (q = 0; q < 3; q++) {
1613 BP[p][q] = 0.0;
1614 for (r = 0; r < 3; r++)
1615 BP[p][q] += (float) drvui->rssC[o1][p][r] * elrot0[r][q]; /* BP is rotated betij matrix */
1616 }
1617 }
1618 for (p = 0; p < 3; p++) {
1619 Z0[p] = 0.0;
1620 for (q = 0; q < 3; q++)
1621 Z0[p] += df[q] * BP[q][p];
1622 }
1623 /* rms0 is the scaled rms displacement value of the ellipsoid in the direction of the bond */
1624 rms0 = factor * (float) sqrt (Z0[0] * df[0] +
1625 Z0[1] * df[1] +
1626 Z0[2] * df[2]);
1627 /* calculate components of end point at intersection of ellipsoid and bond */
1628 for (q = 0; q < 3; q++)
1629 Z0[q] =
1630 s_vert[3 * j + q] + (s_vert[3 * k + q] -
1631 s_vert[3 * j +
1632 q]) * rms0 / d;
1633 }
1634 /* if is_ellipsoid[0] */
1635 if (is_ellipsoid[1] == 1) { /* ellipsoid at other end of bond */
1636 float B[3][3], BP[3][3];
1637
1638 o1 = vert_sym_no[k];
1639 p1 = ellips_id[k]; /* have original ellipsoid number */
1640 jj = ellips_num[k];
1641 if (drvui->ellips[p1].ellips_ismod != 0) { /* apply modulation if needed */
1642 float uij[6];
1643
1644 float vert[3];
1645
1646 vert[0] = o_vert_nm[3 * k];
1647 vert[1] = o_vert_nm[3 * k + 1];
1648 vert[2] = o_vert_nm[3 * k + 2];
1649 if (modulate_uij (vert, p1, jj, o1, uij) == 1) { /* if n.p.d, use average value */
1650 for (p = 0; p < 3; p++)
1651 B[p][p] = drvui->ellips[p1].ellips[p]; /* extract betaij's into square matrix */
1652 B[0][1] = B[1][0] = drvui->ellips[p1].ellips[3];
1653 B[0][2] = B[2][0] = drvui->ellips[p1].ellips[4];
1654 B[1][2] = B[2][1] = drvui->ellips[p1].ellips[5];
1655 } else {
1656 for (p = 0; p < 3; p++)
1657 B[p][p] = uij[p]; // extract betaij's into square matrix
1658 B[0][1] = B[1][0] = uij[3];
1659 B[0][2] = B[2][0] = uij[4];
1660 B[1][2] = B[2][1] = uij[5];
1661 }
1662 } else { /* not modulated */
1663 for (p = 0; p < 3; p++)
1664 B[p][p] = drvui->ellips[p1].ellips[p]; /* extract betaij's into square matrix */
1665 B[0][1] = B[1][0] = drvui->ellips[p1].ellips[3];
1666 B[0][2] = B[2][0] = drvui->ellips[p1].ellips[4];
1667 B[1][2] = B[2][1] = drvui->ellips[p1].ellips[5];
1668 }
1669 for (p = 0; p < 3; p++) {
1670 df[p] = 0.0;
1671 for (q = 0; q < 3; q++) {
1672 elrot0[p][q] = 0.0;
1673 df[p] +=
1674 drvui->ginv[q][p] * (o_vert[3 * k + q] -
1675 o_vert[3 * j + q]);
1676 for (r = 0; r < 3; ++r)
1677 elrot0[p][q] +=
1678 B[p][r] * drvui->rssC[o1][q][r];
1679 }
1680 }
1681 for (p = 0; p < 3; p++) {
1682 for (q = 0; q < 3; q++) {
1683 BP[p][q] = 0.0;
1684 for (r = 0; r < 3; r++)
1685 BP[p][q] += (float) drvui->rssC[o1][p][r] * elrot0[r][q]; /* BP is rotated betij matrix */
1686 }
1687 }
1688 for (p = 0; p < 3; p++) {
1689 Z1[p] = 0.0;
1690 for (q = 0; q < 3; q++)
1691 Z1[p] += df[q] * BP[q][p];
1692 }
1693 /* rms1 is the scaled rms displacement value of the ellipsoid in the direction of the bond */
1694
1695 rms1 = factor * (float) sqrt (Z1[0] * df[0] +
1696 Z1[1] * df[1] +
1697 Z1[2] * df[2]);
1698
1699 /* calculate components of end point at intersection of ellipsoid and bond */
1700 for (q = 0; q < 3; q++)
1701 Z1[q] =
1702 s_vert[3 * k + q] - (s_vert[3 * k + q] -
1703 s_vert[3 * j +
1704 q]) * rms1 / d;
1705 } /* if is_ellipsoid[1] */
1706 }
1707 /* if any segmented ellipsoids present */
1708 for (l1 = 0; l1 < 3; l1++) {
1709 at1[l1] = (float) Z0[l1];
1710 at2[l1] = (float) Z1[l1];
1711 df[l1] = (float) (Z1[l1] - Z0[l1]);
1712 }
1713
1714 if (clipflag != 0) {
1715 out = 0;
1716 l2 = 0;
1717 l3 = 0;
1718 clip0 = 0.0f;
1719 clip1 = 1.0f;
1720 for (l1 = 0; l1 < 3; l1++) {
1721 if (o_vert[3 * j + l1] <
1722 drvui->frames[drvui->frame_no].clip_lim[l1] - 0.01
1723 || o_vert[3 * j + l1] >
1724 drvui->frames[drvui->frame_no].clip_lim[l1 + 3] +
1725 0.01)
1726 l2 = 1;
1727 if (o_vert[3 * k + l1] <
1728 drvui->frames[drvui->frame_no].clip_lim[l1] - 0.01
1729 || o_vert[3 * k + l1] >
1730 drvui->frames[drvui->frame_no].clip_lim[l1 + 3] +
1731 0.01)
1732 l3 = 1;
1733 }
1734 if (l2 && l3)
1735 out = 1; /* both atoms outside */
1736 if (l3)
1737 clip1 = 0.5f;
1738 if (l2)
1739 clip0 = 0.5f;
1740
1741 at2[0] = at1[0] + clip1 * df[0];
1742 at2[1] = at1[1] + clip1 * df[1];
1743 at2[2] = at1[2] + clip1 * df[2];
1744
1745 at1[0] = at1[0] + clip0 * df[0];
1746 at1[1] = at1[1] + clip0 * df[1];
1747 at1[2] = at1[2] + clip0 * df[2];
1748 }
1749 /* if clipflag */
1750 omit = 0;
1751 for (l1 = 0; l1 < Omit->nomits; l1++) {
1752 if (Omit->omit1[l1] == i * 1000
1753 && Omit->omit2[l1] == Bond_Count)
1754 omit = 1;
1755 }
1756
1757 if ((clipflag == 0 || out == 0) && omit == 0) {
1758 if (doPOV) {
1759 fprintf (drvui->fpoutp,
1760 " /* Bond: %c%c %8.5f %8.5f %8.5f TO %c%c %8.5f %8.5f %8.5f*/\n",
1761 bnd1[0], bnd1[1], o_vert[3 * j],
1762 o_vert[3 * j + 1], o_vert[3 * j + 2], bnd2[0],
1763 bnd2[1], o_vert[3 * k], o_vert[3 * k + 1],
1764 o_vert[3 * k + 2]);
1765 if (drvui->bonds[i].bond_style == 0) { /* solid bonds */
1766 fprintf (drvui->fpoutp,
1767 " cylinder{<%8.5f, %8.5f, %8.5f>, <%8.5f, %8.5f, %8.5f>, %8.5f \n",
1768 at1[0], at1[1], at1[2], at2[0], at2[1],
1769 at2[2],
1770 drvui->Bond_Mult *
1771 drvui->bonds[i].bond_size);
1772 fprintf (drvui->fpoutp,
1773 " texture{pigment{color %s }}\n",
1774 col_bond_p);
1775 fprintf (drvui->fpoutp, " }\n");
1776 } else { /* dashed bonds */
1777 numdashes = drvui->bonds[i].bond_style;
1778 dashes = (float) numdashes;
1779
1780 for (m = 0; m < 3; m++)
1781 df[m] = (float) (at2[m] - at1[m]);
1782 fprintf (drvui->fpoutp,
1783 " cylinder{<%8.5f, %8.5f, %8.5f>, <%8.5f, %8.5f, %8.5f>, %8.5f \n",
1784 at1[0], at1[1], at1[2], at1[0] + df[0] / dashes,
1785 at1[1] + df[1] / dashes, at1[2] + df[2] / dashes,
1786 drvui->Bond_Mult *
1787 drvui->bonds[i].bond_size);
1788 fprintf (drvui->fpoutp,
1789 " texture{pigment{color %s }}\n",
1790 col_bond_p);
1791 fprintf (drvui->fpoutp, " }\n");
1792 for (m = 2; m < numdashes - 1; m += 2) {
1793 mf = (float) m;
1794 fprintf (drvui->fpoutp,
1795 " cylinder{<%8.5f, %8.5f, %8.5f>, <%8.5f, %8.5f, %8.5f>, %8.5f \n",
1796 at1[0] + df[0] * mf / dashes,
1797 at1[1] + df[1] * mf / dashes,
1798 at1[2] + df[2] * mf / dashes,
1799 at1[0] + df[0] * (mf + 1.) / dashes,
1800 at1[1] + df[1] * (mf + 1.) / dashes,
1801 at1[2] + df[2] * (mf + 1.) / dashes,
1802 drvui->Bond_Mult *
1803 drvui->bonds[i].bond_size);
1804 fprintf (drvui->fpoutp,
1805 " texture{pigment{color %s }}\n",
1806 col_bond_p);
1807 fprintf (drvui->fpoutp, " }\n");
1808 }
1809 fprintf (drvui->fpoutp,
1810 " cylinder{<%8.5f, %8.5f, %8.5f>, <%8.5f, %8.5f, %8.5f>, %8.5f \n",
1811 at1[0] + df[0] * (dashes - 1.) / dashes,
1812 at1[1] + df[1] * (dashes - 1.) / dashes,
1813 at1[2] + df[2] * (dashes - 1.) / dashes,
1814 at2[0], at2[1], at2[2],
1815 drvui->Bond_Mult *
1816 drvui->bonds[i].bond_size);
1817 fprintf (drvui->fpoutp,
1818 " texture{pigment{color %s }}\n",
1819 col_bond_p);
1820 fprintf (drvui->fpoutp, " }\n");
1821 }
1822 }
1823 if (doAsy) {
1824 (void) sscanf (col_bond_v, "%f %f %f", &glr, &glg, &glb);
1825 fprintf (drvui->fpouta,
1826 " // Bond: %c%c %8.5f %8.5f %8.5f TO %c%c %8.5f %8.5f %8.5f\n",
1827 bnd1[0], bnd1[1], o_vert[3 * j],
1828 o_vert[3 * j + 1], o_vert[3 * j + 2], bnd2[0],
1829 bnd2[1], o_vert[3 * k], o_vert[3 * k + 1],
1830 o_vert[3 * k + 2]);
1831 if (drvui->bonds[i].bond_style == 0) { /* solid bonds */
1832 fprintf (drvui->fpouta," draw(pic, (%8.5f,%8.5f,%8.5f)--(%8.5f,%8.5f,%8.5f),",
1833 at1[0], at1[1], at1[2], at2[0], at2[1], at2[2]);
1834 fprintf (drvui->fpouta, "rgb(%4.2f,%4.2f,%4.2f)+linewidth(%5.2f) );\n",
1835 glr,glg,glb, 25. * drvui->Bond_Mult *
1836 drvui->bonds[i].bond_size);
1837 } else { /* dashed bonds */
1838 numdashes = drvui->bonds[i].bond_style;
1839 dashes = (float) numdashes;
1840
1841 for (m = 0; m < 3; m++)
1842 df[m] = (float) (at2[m] - at1[m]);
1843 fprintf (drvui->fpouta," draw(pic, (%8.5f,%8.5f,%8.5f)--(%8.5f,%8.5f,%8.5f),",
1844 at1[0], at1[1], at1[2], at1[0] + df[0] / dashes,
1845 at1[1] + df[1] / dashes, at1[2] + df[2] / dashes);
1846 fprintf (drvui->fpouta, "rgb(%4.2f,%4.2f,%4.2f)+linewidth(%5.2f)+squarecap );\n",
1847 glr,glg,glb, 25. * drvui->Bond_Mult *
1848 drvui->bonds[i].bond_size);
1849 for (m = 2; m < numdashes - 1; m += 2) {
1850 mf = (float) m;
1851 fprintf (drvui->fpouta," draw(pic, (%8.5f,%8.5f,%8.5f)--(%8.5f,%8.5f,%8.5f),",
1852 at1[0] + df[0] * mf / dashes,
1853 at1[1] + df[1] * mf / dashes,
1854 at1[2] + df[2] * mf / dashes,
1855 at1[0] + df[0] * (mf + 1.) / dashes,
1856 at1[1] + df[1] * (mf + 1.) / dashes,
1857 at1[2] + df[2] * (mf + 1.) / dashes);
1858 fprintf (drvui->fpouta, "rgb(%4.2f,%4.2f,%4.2f)+linewidth(%5.2f)+squarecap );\n",
1859 glr,glg,glb, 25. * drvui->Bond_Mult *
1860 drvui->bonds[i].bond_size);
1861 }
1862 fprintf (drvui->fpouta," draw(pic, (%8.5f,%8.5f,%8.5f)--(%8.5f,%8.5f,%8.5f),",
1863 at1[0] + df[0] * (dashes - 1.) / dashes,
1864 at1[1] + df[1] * (dashes - 1.) / dashes,
1865 at1[2] + df[2] * (dashes - 1.) / dashes,
1866 at2[0], at2[1], at2[2]);
1867 fprintf (drvui->fpouta, "rgb(%4.2f,%4.2f,%4.2f)+linewidth(%5.2f)+squarecap );\n",
1868 glr,glg,glb, 25. * drvui->Bond_Mult *
1869 drvui->bonds[i].bond_size);
1870 }
1871 }
1872 N_Bond++;
1873 if (doVrml && no_comment == 0)
1874 fprintf (drvui->fpoutv,
1875 "# Bond: %c%c %8.5f %8.5f %8.5f TO %c%c %8.5f %8.5f %8.5f\n",
1876 bnd1[0], bnd1[1], o_vert[3 * j],
1877 o_vert[3 * j + 1], o_vert[3 * j + 2], bnd2[0],
1878 bnd2[1], o_vert[3 * k], o_vert[3 * k + 1],
1879 o_vert[3 * k + 2]);
1880 for (m = 0; m < 3; m++)
1881 df[m] = (float) (at2[m] - at1[m]);
1882 d = (float) sqrt (df[0] * df[0] + df[1] * df[1] +
1883 df[2] * df[2]);
1884 /*
1885 calculate rotation matrices needed to rotate bond about Z onto Y-Z plane and then about X so that
1886 bond is parallel to Z. Finally, apply inverse transformation to cylinder, which started parallel
1887 to X. Resulting cylinder will be parallel to vector between atoms.
1888 */
1889 beta = (float) atan2 (df[0], df[1] + 0.0000001f); /* rotation angle about Z (in radians) */
1890 gamma = (float) sqrt (df[0] * df[0] + df[1] * df[1]);
1891 gamma = (float) atan2 (gamma, df[2]); /* Rotation angle about X */
1892 glPushMatrix ();
1893 glLoadName (i * 1000);
1894 glPushName (Bond_Count);
1895 (void) sscanf (col_bond_v, "%f %f %f", &glr, &glg, &glb);
1896 glColor3f (glr, glg, glb);
1897 if (drvui->bonds[i].bond_style == 0) {
1898 glPushMatrix ();
1899 glTranslatef (at1[0], at1[1], at1[2]);
1900 glRotatef (-beta * (float) RAD, 0.0f, 0.0f, 1.);
1901 glRotatef (-gamma * (float) RAD, 1., 0.0f, 0.0f);
1902 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
1903 gluCylinder (glu_quadObj,
1904 drvui->Bond_Mult * drvui->bonds[i].bond_size,
1905 drvui->Bond_Mult * drvui->bonds[i].bond_size,
1906 d, 10, 1);
1907 glPopMatrix ();
1908 } else {
1909 numdashes = drvui->bonds[i].bond_style;
1910 dashes = (float) numdashes;
1911 glPushMatrix ();
1912 glTranslatef (at1[0], at1[1], at1[2]);
1913 glRotatef (-beta * (float) RAD, 0.0f, 0.0f, 1.);
1914 glRotatef (-gamma * (float) RAD, 1., 0.0f, 0.0f);
1915 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
1916 gluCylinder (glu_quadObj,
1917 drvui->Bond_Mult * drvui->bonds[i].bond_size,
1918 drvui->Bond_Mult * drvui->bonds[i].bond_size,
1919 d / dashes, 10, 1);
1920 glPopMatrix ();
1921 for (m = 2; m < numdashes; m += 2) {
1922 mf = (float) m;
1923 glPushMatrix ();
1924 glTranslatef (at1[0] + df[0] * mf / dashes,
1925 at1[1] + df[1] * mf / dashes,
1926 at1[2] + df[2] * mf / dashes);
1927 glRotatef (-beta * (float) RAD, 0.0f, 0.0f, 1.);
1928 glRotatef (-gamma * (float) RAD, 1., 0.0f, 0.0f);
1929 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
1930 gluCylinder (glu_quadObj,
1931 drvui->Bond_Mult *
1932 drvui->bonds[i].bond_size,
1933 drvui->Bond_Mult *
1934 drvui->bonds[i].bond_size, d / dashes,
1935 10, 1);
1936 glPopMatrix ();
1937 }
1938 }
1939 glPopName ();
1940 glPopMatrix ();
1941 if (doVrml) {
1942 if (Vrml2) {
1943 if (drvui->bonds[i].bond_style == 0) { /* solid bond */
1944 fprintf (drvui->fpoutv, " Transform{\n ");
1945 fprintf (drvui->fpoutv,
1946 " translation %5.3f %5.3f %5.3f\n",
1947 (at1[0] + at2[0]) / 2.,
1948 (at1[1] + at2[1]) / 2.,
1949 (at1[2] + at2[2]) / 2.);
1950 fprintf (drvui->fpoutv, " rotation 0 0 1 %f\n",
1951 -beta);
1952 fprintf (drvui->fpoutv,
1953 " children Transform { rotation 1 0 0 %f\n",
1954 1.5708 - gamma);
1955 fprintf (drvui->fpoutv, " children [ Shape {\n");
1956 fprintf (drvui->fpoutv,
1957 " appearance Appearance { material Material { diffuseColor %s} }\n",
1958 col_bond_v);
1959 fprintf (drvui->fpoutv, " geometry Cylinder {");
1960 fprintf (drvui->fpoutv, "radius %f",
1961 drvui->Bond_Mult *
1962 drvui->bonds[i].bond_size);
1963 fprintf (drvui->fpoutv, " height %f", d);
1964 fprintf (drvui->fpoutv, "}}]}}\n");
1965 } else { /* dashed bonds */
1966 numdashes = drvui->bonds[i].bond_style;
1967 dashes = (float) numdashes;
1968 fprintf (drvui->fpoutv, " Transform{\n ");
1969 fprintf (drvui->fpoutv,
1970 " translation %5.3f %5.3f %5.3f\n",
1971 at1[0] + df[0] / dashes,
1972 at1[1] + df[1] / dashes,
1973 at1[2] + df[2] / dashes);
1974 fprintf (drvui->fpoutv, " rotation 0 0 1 %f\n",
1975 -beta);
1976 fprintf (drvui->fpoutv,
1977 " children Transform { rotation 1 0 0 %f\n",
1978 1.5708 - gamma);
1979 fprintf (drvui->fpoutv, " children [ Shape {\n");
1980 fprintf (drvui->fpoutv,
1981 " appearance Appearance { material Material { diffuseColor %s} }\n",
1982 col_bond_v);
1983 fprintf (drvui->fpoutv, " geometry Cylinder {");
1984 fprintf (drvui->fpoutv, "radius %f",
1985 drvui->Bond_Mult *
1986 drvui->bonds[i].bond_size);
1987 fprintf (drvui->fpoutv, " height %f",
1988 d / dashes);
1989 fprintf (drvui->fpoutv, "}}]}}\n");
1990 for (m = 2; m < dashes - 1; m += 2) {
1991 mf = (float) m;
1992 fprintf (drvui->fpoutv, " Transform{\n ");
1993 fprintf (drvui->fpoutv,
1994 " translation %5.3f %5.3f %5.3f\n",
1995 at1[0] +
1996 df[0] * mf / dashes,
1997 at1[1] +
1998 df[1] * mf / dashes,
1999 at1[2] +
2000 df[2] * mf / dashes);
2001 fprintf (drvui->fpoutv,
2002 " rotation 0 0 1 %f\n", -beta);
2003 fprintf (drvui->fpoutv,
2004 " children Transform { rotation 1 0 0 %f\n",
2005 1.5708 - gamma);
2006 fprintf (drvui->fpoutv,
2007 " children [ Shape {\n");
2008 fprintf (drvui->fpoutv,
2009 " appearance Appearance { material Material { diffuseColor %s} }\n",
2010 col_bond_v);
2011 fprintf (drvui->fpoutv,
2012 " geometry Cylinder {");
2013 fprintf (drvui->fpoutv, "radius %f",
2014 drvui->Bond_Mult *
2015 drvui->bonds[i].bond_size);
2016 fprintf (drvui->fpoutv, " height %f",
2017 d / dashes);
2018 fprintf (drvui->fpoutv, "}}]}}\n");
2019 }
2020 }
2021 } else { /* VRML1 */
2022 if (drvui->bonds[i].bond_style == 0) { /* solid bonds */
2023 fprintf (drvui->fpoutv, " Separator{\n");
2024 fprintf (drvui->fpoutv, " Transform{\n ");
2025 fprintf (drvui->fpoutv,
2026 " translation %5.3f %5.3f %5.3f\n }\n",
2027 (at1[0] + at2[0]) / 2.,
2028 (at1[1] + at2[1]) / 2.,
2029 (at1[2] + at2[2]) / 2.);
2030 fprintf (drvui->fpoutv, " Rotation {\n");
2031 fprintf (drvui->fpoutv,
2032 " rotation 0 0 1 %f\n }\n", -beta);
2033 fprintf (drvui->fpoutv, " Rotation {\n");
2034 fprintf (drvui->fpoutv,
2035 " rotation 1 0 0 %f\n }\n",
2036 1.5708 - gamma);
2037 fprintf (drvui->fpoutv,
2038 " Material {\n diffuseColor %s\n",
2039 col_bond_v);
2040 fprintf (drvui->fpoutv, " }\n");
2041 fprintf (drvui->fpoutv, " Cylinder {\n");
2042 fprintf (drvui->fpoutv, " parts SIDES\n");
2043 fprintf (drvui->fpoutv, " radius %f\n",
2044 drvui->Bond_Mult *
2045 drvui->bonds[i].bond_size);
2046 fprintf (drvui->fpoutv, " height %f\n", d);
2047 fprintf (drvui->fpoutv, " }\n }\n");
2048 } else { /* dashed bonds */
2049 fprintf (drvui->fpoutv, " Separator{\n");
2050 fprintf (drvui->fpoutv, " Transform{\n ");
2051 fprintf (drvui->fpoutv,
2052 " translation %5.3f %5.3f %5.3f\n }\n",
2053 at1[0] + df[0] / dashes,
2054 at1[1] + df[1] / dashes,
2055 at1[2] + df[2] / dashes);
2056 fprintf (drvui->fpoutv, " Rotation {\n");
2057 fprintf (drvui->fpoutv,
2058 " rotation 0 0 1 %f\n }\n", -beta);
2059 fprintf (drvui->fpoutv, " Rotation {\n");
2060 fprintf (drvui->fpoutv,
2061 " rotation 1 0 0 %f\n }\n",
2062 1.5708 - gamma);
2063 fprintf (drvui->fpoutv,
2064 " Material {\n diffuseColor %s\n",
2065 col_bond_v);
2066 fprintf (drvui->fpoutv, " }\n");
2067 fprintf (drvui->fpoutv, " Cylinder {\n");
2068 fprintf (drvui->fpoutv, " parts SIDES\n");
2069 fprintf (drvui->fpoutv, " radius %f\n",
2070 drvui->Bond_Mult *
2071 drvui->bonds[i].bond_size);
2072 fprintf (drvui->fpoutv, " height %f\n", d / dashes);
2073 fprintf (drvui->fpoutv, " }\n }\n");
2074 for (m = 2; m < dashes - 1; m += 2) {
2075 mf = (float) m;
2076 fprintf (drvui->fpoutv, " Separator{\n");
2077 fprintf (drvui->fpoutv, " Transform{\n ");
2078 fprintf (drvui->fpoutv,
2079 " translation %5.3f %5.3f %5.3f\n }\n",
2080 at1[0] +
2081 df[0] * mf / dashes,
2082 at1[1] +
2083 df[1] * mf / dashes,
2084 at1[2] +
2085 df[2] * mf / dashes);
2086 fprintf (drvui->fpoutv, " Rotation {\n");
2087 fprintf (drvui->fpoutv,
2088 " rotation 0 0 1 %f\n }\n", -beta);
2089 fprintf (drvui->fpoutv, " Rotation {\n");
2090 fprintf (drvui->fpoutv,
2091 " rotation 1 0 0 %f\n }\n",
2092 1.5708 - gamma);
2093 fprintf (drvui->fpoutv,
2094 " Material {\n diffuseColor %s\n",
2095 col_bond_v);
2096 fprintf (drvui->fpoutv, " }\n");
2097 fprintf (drvui->fpoutv, " Cylinder {\n");
2098 fprintf (drvui->fpoutv, " parts SIDES\n");
2099 fprintf (drvui->fpoutv, " radius %f\n",
2100 drvui->Bond_Mult *
2101 drvui->bonds[i].bond_size);
2102 fprintf (drvui->fpoutv, " height %f\n", d / dashes);
2103 fprintf (drvui->fpoutv, " }\n }\n");
2104 }
2105 }
2106 }
2107 }
2108 }
2109 } /* if distance OK */
2110 } /* for k */
2111 } /* for j */
2112 } /* if nvert */
2113 } /* for i (bond types) */
2114 free (ellips_id);
2115 free (ellips_num);
2116 gluDeleteQuadric (glu_quadObj);
2117 fprintf (drvui->fcns, "%4d bonds.\n", Bond_Count);
2118 fprintf (drvui->flout, "\nGenerated %4d bonds.\n\n", Bond_Count);
2119 } /* end of generate_bonds */
2120
2121 /* ************************************************************** */
2122 /* ************************************************************** */
2123
2124 void
generate_cones(void)2125 generate_cones (void)
2126
2127 /* routine to generate lone pair cone descriptions and add them to edit list */
2128 {
2129 int start, same_type;
2130
2131 int i, j, k, l; /* loop counters */
2132
2133 int N_Bond; /* Counter for number of bonds of a type */
2134
2135 int l1; /* loop counter for matrix operations */
2136
2137 int nvert1; /* vertex counter */
2138
2139 int Bond_Count; /* Counter for bonds output */
2140
2141 char bnd1[4]; /* temporary storage of bond labels */
2142
2143 float d = 0; /* bond distance */
2144
2145 float df[3]; /* vector difference */
2146
2147 float beta; /* rotation angle */
2148
2149 float gamma; /* another rotation angle */
2150
2151 float at1[3], at2[3]; /*initial coordinates of bond endpoints */
2152
2153 float clip0, clip1; /* clipping scale factors */
2154
2155 int m, l2, l3, out = 0;
2156
2157 int p; /* loop variables for rotation matrix calculation */
2158
2159 float Z0[3], Z1[3]; /* components of rms displacement along bond */
2160
2161 float mybonds[20][20]; /*temporary storage for atoms bonded to lonepair site */
2162
2163 int nbnds = 0;
2164
2165 float vecsum[3], vecpro[3];
2166
2167 float glr, glg, glb;
2168
2169 GLUquadricObj *glu_quadObj;
2170
2171 char col_cone_p[40];
2172
2173 char col_cone_v[40];
2174
2175 if (drvui->ncone == 1)
2176 return; /* exit if no work to be done */
2177
2178 Bond_Count = 0;
2179
2180 glu_quadObj = gluNewQuadric ();
2181
2182 for (i = 1; i < drvui->ncone; ++i) { // loop through cone types
2183 if (drvui->cones[i].cone_fn != drvui->frame_no)
2184 continue; // skip if wrong frame
2185 strcpy (col_cone_p, drvui->cones[i].col_cone);
2186 strcpy (col_cone_v, col_cone_p);
2187 Transform_VRML_Color (col_cone_v);
2188 Transform_POV_Color (col_cone_p);
2189 nvert = 0; // empty the vertex list
2190 for (j = 0; j < natom; ++j) { // loop through atom types
2191 if (drvui->atoms[j].atom_fn != drvui->frame_no)
2192 continue;
2193 if (((drvui->atoms[j].atom_n & 255) > 0) || (((drvui->atoms[j].atom_n >> 24) & 255) > 0)) { // sphere or ellipsoid
2194 int o;
2195
2196 o = drvui->atoms[j].atom_n >> 24;
2197 if (o > 0) { // ellipsoid
2198 if (check_atom_name
2199 (drvui->cones[i].cone_l1, drvui->ellips[o].ellips_l))
2200 find_all_in_box (j); /* get all atoms of type j */
2201 } else {
2202 o = drvui->atoms[j].atom_n & 255;
2203 if (check_atom_name
2204 (drvui->cones[i].cone_l1, drvui->spheres[o].sphere_l)) {
2205 find_all_in_box (j); /* get all atoms of type j */
2206 }
2207 } // o > 0
2208 }
2209 } /* for (j=0;j<natom;++j) */
2210
2211 bnd1[0] = drvui->cones[i].cone_l1[0];
2212 bnd1[1] = drvui->cones[i].cone_l1[1];
2213 bnd1[2] = drvui->cones[i].cone_l1[2];
2214 bnd1[3] = drvui->cones[i].cone_l1[3];
2215 nvert1 = nvert;
2216 same_type = 0;
2217
2218 /* atoms i and j are not the same type */
2219 for (j = 0; j < natom; ++j) { /* loop through atom types */
2220 if (drvui->atoms[j].atom_fn != drvui->frame_no)
2221 continue;
2222 if (((drvui->atoms[j].atom_n & 255) > 0) || (((drvui->atoms[j].atom_n >> 24) & 255) > 0)) { // sphere or ellipsoid
2223 find_all_in_box (j); /* get all atoms of type j */
2224 }
2225 } /* for (j=0;j<natom;++j) */
2226
2227 if (((nvert1 > 0) && (nvert > nvert1)) || (same_type == 1)) {
2228 N_Bond = 0;
2229 nbnds = 0;
2230 for (j = 0; j < nvert1; ++j) { /* loop through atoms of type 1 */
2231 nbnds = 0;
2232 start = nvert1;
2233 if (same_type)
2234 start = j;
2235 start = 0;
2236 for (k = start; k < nvert; ++k) { /* loop on type 2 atoms */
2237 d = dist (j, k); /* get atomic distance */
2238 // if (d > 0.2 && d <= drvui->cones[i].cone_height) { /* here if distance within range */
2239 if (d > 0.2 && d <= printdist) { /* here if distance within range */
2240 mybonds[0][nbnds] = s_vert[3 * k];
2241 mybonds[1][nbnds] = s_vert[3 * k + 1];
2242 mybonds[2][nbnds] = s_vert[3 * k + 2];
2243 nbnds++;
2244 } /* if distance ok */
2245 } /* for k */
2246 /*# */
2247 if (nbnds < 4 - abs (drvui->cones[i].numlonepairs)) {
2248 continue;
2249 }
2250 df[0] = df[1] = drvui->Bond_Mult * drvui->bonds[i].bond_size; /* scaling */
2251 df[2] = 0.5f * d;
2252
2253 Bond_Count++; /* update count */
2254 for (p = 0; p < 3; p++) { /* initialize endpoints of bond */
2255 Z0[p] = s_vert[3 * j + p];
2256 }
2257
2258 for (l = 0; l < abs (drvui->cones[i].numlonepairs); l++) {
2259 for (p = 0; p < 3; p++) { /* initialize endpoints of bond */
2260 Z1[p] = 0;
2261
2262 if (drvui->cones[i].numlonepairs == 1) {
2263 for (k = 0; k < nbnds; k++)
2264 Z1[p] -= mybonds[p][k] - s_vert[3 * j + p];
2265 Z1[p] = Z0[p] + Z1[p];
2266 } else if (drvui->cones[i].numlonepairs < 0) {
2267 vecpro[p] = 0.;
2268 vecpro[0] = (s_vert[3 * j + 1] - mybonds[1][0])
2269 * (s_vert[3 * j + 2] - mybonds[2][1])
2270 - (s_vert[3 * j + 1] - mybonds[1][1])
2271 * (s_vert[3 * j + 2] - mybonds[2][0]);
2272 vecpro[1] = (s_vert[3 * j + 2] - mybonds[2][0])
2273 * (s_vert[3 * j] - mybonds[0][1])
2274 - (s_vert[3 * j + 2] - mybonds[2][1])
2275 * (s_vert[3 * j] - mybonds[0][0]);
2276 vecpro[2] = (s_vert[3 * j] - mybonds[0][0])
2277 * (s_vert[3 * j + 1] - mybonds[1][1])
2278 - (s_vert[3 * j] - mybonds[0][1])
2279 * (s_vert[3 * j + 1] - mybonds[1][0]);
2280 if (l > 0)
2281 Z1[p] = Z0[p] + vecpro[p];
2282 else
2283 Z1[p] = Z0[p] - vecpro[p];
2284 } else if (drvui->cones[i].numlonepairs == 2) {
2285 vecsum[p] = 0.;
2286 vecpro[p] = 0.;
2287 for (k = 0; k < nbnds; k++) {
2288 vecsum[p] += mybonds[p][k] - s_vert[3 * j + p];
2289 }
2290
2291 vecpro[0] = (s_vert[3 * j + 1] - mybonds[1][0])
2292 * (s_vert[3 * j + 2] - mybonds[2][1])
2293 - (s_vert[3 * j + 1] - mybonds[1][1])
2294 * (s_vert[3 * j + 2] - mybonds[2][0]);
2295 vecpro[1] = (s_vert[3 * j + 2] - mybonds[2][0])
2296 * (s_vert[3 * j] - mybonds[0][1])
2297 - (s_vert[3 * j + 2] - mybonds[2][1])
2298 * (s_vert[3 * j] - mybonds[0][0]);
2299 vecpro[2] = (s_vert[3 * j] - mybonds[0][0])
2300 * (s_vert[3 * j + 1] - mybonds[1][1])
2301 - (s_vert[3 * j] - mybonds[0][1])
2302 * (s_vert[3 * j + 1] - mybonds[1][0]);
2303 switch (l) {
2304 case 0:
2305 Z1[p] = -vecsum[p] + vecpro[p];
2306 break;
2307 case 1:
2308 Z1[p] = -vecsum[p] - vecpro[p];
2309 break;
2310 }
2311 Z1[p] = Z0[p] + Z1[p];
2312 }
2313 } /* for p */
2314
2315 d = (float) sqrt ((Z0[0] - Z1[0]) * (Z0[0] - Z1[0])
2316 + (Z0[1] - Z1[1]) * (Z0[1] - Z1[1])
2317 + (Z0[2] - Z1[2]) * (Z0[2] - Z1[2]));
2318 d = drvui->cones[i].cone_height / d;
2319 Z1[0] = Z0[0] + d * (Z1[0] - Z0[0]);
2320 Z1[1] = Z0[1] + d * (Z1[1] - Z0[1]);
2321 Z1[2] = Z0[2] + d * (Z1[2] - Z0[2]);
2322
2323 for (l1 = 0; l1 < 3; l1++) {
2324 at1[l1] = (float) Z0[l1];
2325 at2[l1] = (float) Z1[l1];
2326 df[l1] = (float) (Z1[l1] - Z0[l1]);
2327 }
2328
2329 if (clipflag != 0) {
2330 out = 0;
2331 l2 = 0;
2332 l3 = 0;
2333 clip0 = 0.0f;
2334 clip1 = 1.0f;
2335 for (l1 = 0; l1 < 3; l1++) {
2336 if (o_vert[3 * j + l1] <
2337 drvui->frames[drvui->frame_no].clip_lim[l1] - 0.01
2338 || o_vert[3 * j + l1] >
2339 drvui->frames[drvui->frame_no].clip_lim[l1 + 3] + 0.01)
2340 l2 = 1;
2341 if (o_vert[3 * k + l1] <
2342 drvui->frames[drvui->frame_no].clip_lim[l1] - 0.01
2343 || o_vert[3 * k + l1] >
2344 drvui->frames[drvui->frame_no].clip_lim[l1 + 3] + 0.01)
2345 l3 = 1;
2346 }
2347 if (l2 && l3)
2348 out = 1; /* both atoms outside */
2349 if (l2)
2350 out = 1;
2351 if (l3)
2352 clip1 = 0.5f;
2353 if (l2)
2354 clip0 = 0.5f;
2355
2356 at1[0] = at1[0] + clip0 * df[0];
2357 at1[1] = at1[1] + clip0 * df[1];
2358 at1[2] = at1[2] + clip0 * df[2];
2359
2360 at2[0] = at1[0] + clip1 * df[0];
2361 at2[1] = at1[1] + clip1 * df[1];
2362 at2[2] = at1[2] + clip1 * df[2];
2363
2364 }
2365 /* if clipflag */
2366 if (clipflag == 0 || out == 0) {
2367 if (doPOV) {
2368 fprintf (drvui->fpoutp,
2369 " /* LonePair on %c%c %8.5f %8.5f %8.5f */\n",
2370 bnd1[0], bnd1[1], o_vert[3 * j], o_vert[3 * j + 1],
2371 o_vert[3 * j + 2]);
2372 fprintf (drvui->fpoutp,
2373 " cone{<%8.5f, %8.5f, %8.5f>, %8.5f, <%8.5f, %8.5f, %8.5f>, %8.5f \n",
2374 at1[0], at1[1], at1[2],
2375 drvui->Bond_Mult * drvui->cones[i].cone_min, at2[0],
2376 at2[1], at2[2],
2377 drvui->Bond_Mult * drvui->cones[i].cone_max);
2378 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n",
2379 col_cone_p);
2380 fprintf (drvui->fpoutp, " }\n");
2381 fprintf (drvui->fpoutp,
2382 " sphere{<%8.5f, %8.5f, %8.5f> %8.5f \n", at2[0],
2383 at2[1], at2[2],
2384 drvui->Bond_Mult * drvui->cones[i].cone_max);
2385 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n",
2386 col_cone_p);
2387 fprintf (drvui->fpoutp, " }\n");
2388 }
2389 N_Bond++;
2390 if (no_comment == 0 && doVrml)
2391 fprintf (drvui->fpoutv,
2392 "# LonePair on: %c%c %8.5f %8.5f %8.5f\n", bnd1[0],
2393 bnd1[1], o_vert[3 * j], o_vert[3 * j + 1],
2394 o_vert[3 * j + 2]);
2395 for (m = 0; m < 3; m++)
2396 df[m] = (float) (at2[m] - at1[m]);
2397 d = (float) sqrt (df[0] * df[0] + df[1] * df[1] + df[2] * df[2]);
2398 /*
2399 calculate rotation matrices needed to rotate bond about Z onto Y-Z plane and then about X so that
2400 bond is parallel to Z. Finally, apply inverse transformation to cylinder, which started parallel
2401 to X. Resulting cylinder will be parallel to vector between atoms.
2402 */
2403 beta = (float) atan2 (df[0], df[1] + 0.0000001f); /* rotation angle about Z (in radians) */
2404 gamma = (float) sqrt (df[0] * df[0] + df[1] * df[1]);
2405 gamma = (float) atan2 (gamma, df[2]); /* Rotation angle about X */
2406 glPushMatrix ();
2407 sscanf (col_cone_v, "%f %f %f", &glr, &glg, &glb);
2408 glColor3f (glr, glg, glb);
2409 glTranslatef (at1[0], at1[1], at1[2]);
2410 glRotatef (-beta * (float) RAD, 0.0f, 0.0f, 1.0f);
2411 glRotatef (-gamma * (float) RAD, 1.0f, 0.0f, 0.0f);
2412 gluQuadricDrawStyle (glu_quadObj, GLU_FILL);
2413 gluCylinder (glu_quadObj,
2414 drvui->Bond_Mult * drvui->cones[i].cone_min,
2415 drvui->Bond_Mult * drvui->cones[i].cone_max, d, 10,
2416 1);
2417 glPopMatrix ();
2418 glPushMatrix ();
2419 glTranslatef (at2[0], at2[1], at2[2]);
2420 glutSolidSphere (drvui->cones[i].cone_max, 10, 10);
2421 glPopMatrix ();
2422 if (doVrml) {
2423 if (Vrml2) {
2424 fprintf (drvui->fpoutv, " Transform{\n ");
2425 fprintf (drvui->fpoutv,
2426 " translation %5.3f %5.3f %5.3f\n",
2427 (at1[0] + at2[0]) / 2., (at1[1] + at2[1]) / 2.,
2428 (at1[2] + at2[2]) / 2.);
2429 fprintf (drvui->fpoutv, " rotation 0 0 1 %f\n", -beta);
2430 fprintf (drvui->fpoutv,
2431 " children Transform { rotation 0 1 0 %f\n",
2432 3.14);
2433 fprintf (drvui->fpoutv,
2434 " children Transform { rotation 1 0 0 %f\n",
2435 1.5708 + gamma);
2436 fprintf (drvui->fpoutv, " children [ Shape {\n");
2437 fprintf (drvui->fpoutv,
2438 " appearance Appearance { material Material { diffuseColor %s} }\n",
2439 col_cone_v);
2440 fprintf (drvui->fpoutv, " geometry Cone {");
2441 fprintf (drvui->fpoutv, "bottomRadius %f",
2442 drvui->Bond_Mult * drvui->cones[i].cone_max);
2443 fprintf (drvui->fpoutv, " height %f", d);
2444 fprintf (drvui->fpoutv, "}}]}}}\n");
2445 fprintf (drvui->fpoutv, " Transform {");
2446 fprintf (drvui->fpoutv,
2447 " translation %8.5f %8.5f %8.5f\n", at2[0],
2448 at2[1], at2[2]);
2449 fprintf (drvui->fpoutv, " children [\n Shape {");
2450 fprintf (drvui->fpoutv,
2451 "appearance Appearance {material Material { diffuseColor %s} }\n",
2452 col_cone_v);
2453 fprintf (drvui->fpoutv,
2454 " geometry Sphere { radius %8.5f}\n",
2455 drvui->cones[i].cone_max);
2456 fprintf (drvui->fpoutv, "}]}\n");
2457 } else {
2458 fprintf (drvui->fpoutv, " Separator{\n");
2459 fprintf (drvui->fpoutv, " Transform{\n ");
2460 fprintf (drvui->fpoutv,
2461 " translation %5.3f %5.3f %5.3f\n }\n",
2462 (at1[0] + at2[0]) / 2., (at1[1] + at2[1]) / 2.,
2463 (at1[2] + at2[2]) / 2.);
2464 fprintf (drvui->fpoutv, " Rotation {\n");
2465 fprintf (drvui->fpoutv, " rotation 0 0 1 %f\n }\n",
2466 -beta);
2467 fprintf (drvui->fpoutv, " Rotation {\n");
2468 fprintf (drvui->fpoutv, " rotation 0 1 0 %f\n }\n",
2469 3.14);
2470 fprintf (drvui->fpoutv, " Rotation {\n");
2471 fprintf (drvui->fpoutv, " rotation 1 0 0 %f\n }\n",
2472 1.5708 + gamma);
2473 fprintf (drvui->fpoutv,
2474 " Material {\n diffuseColor %s\n",
2475 col_cone_v);
2476 fprintf (drvui->fpoutv, " }\n");
2477 fprintf (drvui->fpoutv, " Cone {\n");
2478 fprintf (drvui->fpoutv, " parts SIDES\n");
2479 fprintf (drvui->fpoutv, " bottomRadius %f\n",
2480 drvui->Bond_Mult * drvui->cones[i].cone_max);
2481 fprintf (drvui->fpoutv, " height %f\n", d);
2482 fprintf (drvui->fpoutv, " }\n }\n");
2483 fprintf (drvui->fpoutv, " Separator {");
2484 fprintf (drvui->fpoutv, " Transform {");
2485 fprintf (drvui->fpoutv,
2486 " translation %8.5f %8.5f %8.5f}\n", at2[0],
2487 at2[1], at2[2]);
2488 fprintf (drvui->fpoutv, " Material {");
2489 fprintf (drvui->fpoutv, " shininess 0.3\n");
2490 fprintf (drvui->fpoutv, " diffuseColor %s}\n",
2491 col_cone_v);
2492 fprintf (drvui->fpoutv, " Sphere { radius %8.5f}\n",
2493 drvui->cones[i].cone_max);
2494 fprintf (drvui->fpoutv, " }\n");
2495 } /* if (Vrml2) */
2496 }
2497 if (doAsy) {
2498 fprintf (drvui->fpouta, "draw (pic, shift(%5.2f,%5.2f,%5.2f)*rotate(%.2f,Z)*rotate(%.2f,X)*\n",
2499 at2[0],at2[1],at2[2],-beta*RAD,-gamma*RAD);
2500 fprintf (drvui->fpouta, "\tscale(%.2f,%.2f,%.2f)*unitcone,rgb(%.2f,%.2f,%.2f));\n",
2501 drvui->Bond_Mult * drvui->cones[i].cone_max,drvui->Bond_Mult * drvui->cones[i].cone_max,-d,glr,glg,glb);
2502 fprintf (drvui->fpouta, "draw (pic, shift(%5.2f,%5.2f,%5.2f)*scale3(%.2f)*unitsphere,rgb(%.2f,%.2f,%.2f));\n",
2503 at2[0],at2[1],at2[2],drvui->cones[i].cone_max,glr,glg,glb);
2504 }
2505 }
2506 /* if not clipped */
2507 } /* for l */
2508 } /* for j */
2509 } /* if nvert */
2510 } /* for i */
2511 gluDeleteQuadric (glu_quadObj);
2512 fprintf (drvui->fcns, "%4d cones.\n", Bond_Count);
2513 fprintf (drvui->flout, "\nGenerated %4d cones.\n\n", Bond_Count);
2514 } /* end of generate_cones */
2515
2516 /* ************************************************************** */
2517 /* ************************************************************** */
2518
2519 void
get_atom_id(void)2520 get_atom_id (void)
2521
2522 /* routine to rework atom numbers to point at polyhedron, sphere
2523 plane, or ellipsoid characteristics list */
2524 {
2525 int i, j; /* local variables */
2526
2527 /* lines commented out with a question mark would be needed for correctness, but
2528 would break the stishovite example that relies on polyhedral completion carrying
2529 over from the first frame into the second */
2530
2531 for (i = 0; i < natom; ++i) {
2532 if (drvui->atoms[i].atom_fn != drvui->frame_no)
2533 continue;
2534 drvui->atoms[i].sv_atom_n = drvui->atoms[i].atom_n; /* save original atom number */
2535 drvui->atoms[i].atom_n = 0; // clear encoded location
2536 for (j = 1; j < drvui->nsphere; ++j) {
2537 //? if (drvui->sphere_fn[j] != drvui->frame_no) continue;
2538 if (check_atom_name (drvui->atoms[i].atom_l, drvui->spheres[j].sphere_l)
2539 && (drvui->atoms[i].sv_atom_n == drvui->spheres[j].sphere_n
2540 || drvui->spheres[j].sphere_n == -1)) {
2541 drvui->atoms[i].atom_n = j; // atom is in sphere list
2542 }
2543 }
2544 for (j = 1; j < drvui->npoly; ++j) {
2545 //? if (drvui->poly_fn[j] != drvui->frame_no) continue;
2546 if (check_atom_name (drvui->atoms[i].atom_l, drvui->polyhedra[j].poly_l)) {
2547 drvui->atoms[i].atom_n += (j << 8); // atom is in poly. list
2548 }
2549 }
2550 for (j = 1; j < drvui->nplane; ++j) {
2551 //? if (drvui->plane_fn[j] != drvui->frame_no) continue;
2552 if (check_atom_name (drvui->atoms[i].atom_l, drvui->planes[j].plane_l)) {
2553 drvui->atoms[i].atom_n += (j << 16); // atom is in plane list
2554 }
2555 }
2556 if (drvui->do_ellipsoids) {
2557 for (j = 1; j < drvui->n_ellips; ++j) {
2558 if (check_atom_name (drvui->atoms[i].atom_l, drvui->ellips[j].ellips_l)
2559 && (drvui->atoms[i].sv_atom_n == drvui->ellips[j].ellips_n)) {
2560 drvui->atoms[i].atom_n += (j << 24); // atom is in ellipsoid list
2561 }
2562 }
2563 }
2564 }
2565 } // end of get_atom_id
2566
2567 /* ************************************************************** */
2568 /* ************************************************************** */
2569
2570 void
get_input(int Quick)2571 get_input (int Quick)
2572
2573 /* Routine to process input file, make lattice metric and do preliminary
2574 ellipsoid processing */
2575 {
2576 int i, j, k;
2577
2578 float biso;
2579
2580 if (!Quick)
2581 fprintf (drvui->flout,
2582 "\n********************* Input: ************************\n");
2583 read_inp (Quick);
2584 if (Quick)
2585 return;
2586 fprintf (drvui->flout, "*****************************************************\n\n");
2587 if (!Labels)
2588 Display_axes = 0; /* If no labels, no vectors */
2589 if (drvui->frame_no == 1 && doVrml) {
2590 if (Vrml2)
2591 fprintf (drvui->fcns, "\nGenerating VRML97 Output.\n");
2592 else
2593 fprintf (drvui->fcns, "\nGenerating VRML1.0 Output.\n");
2594 }
2595 print_sym (); /* print symmetry information */
2596 if (drvui->xyzoff_read)
2597 drvui->origin1_flag = 0;
2598 (void) fclose (drvui->fpin); /* close file */
2599 make_bmat (drvui->sys, drvui->lat_con, drvui->b_mat, drvui->ginv, drvui->rec_lat_con); /* create the lattice metric */
2600 fprintf (drvui->flout, "Lattice metrics\n");
2601 fprintf (drvui->flout, "%8.3f %8.3f %8.3f\n", drvui->b_mat[0][0], drvui->b_mat[0][1],
2602 drvui->b_mat[0][2]);
2603 fprintf (drvui->flout, "%8.3f %8.3f %8.3f\n", drvui->b_mat[1][0], drvui->b_mat[1][1],
2604 drvui->b_mat[1][2]);
2605 fprintf (drvui->flout, "%8.3f %8.3f %8.3f\n", drvui->b_mat[2][0], drvui->b_mat[2][1],
2606 drvui->b_mat[2][2]);
2607 /* set limits of search regions */
2608 if (packflag == 0 && boxflag == 0 && drvui->frame_no == 1) { /* neither given and frame 1 */
2609 for (i = 0; i < 3; i++) {
2610 drvui->frames[drvui->frame_no].cryst_lim[i] = origin[i] - 0.5f; /* set to generate one cell */
2611 drvui->frames[drvui->frame_no].cryst_lim[i + 3] = origin[i] + 0.5f;
2612 }
2613 packflag = 1;
2614 }
2615 if (packflag != 0) { /* pack command given - calculate approximate bounds */
2616 float temp, temp1;
2617
2618 for (i = 0; i < 3; i++) {
2619 temp = origin[i] - drvui->frames[drvui->frame_no].cryst_lim[i];
2620 temp1 = drvui->frames[drvui->frame_no].cryst_lim[i + 3] - origin[i];
2621 if (temp1 > temp)
2622 temp = temp1;
2623 boxlim[i] = (temp + 0.2f) * drvui->lat_con[i] / 0.5f;
2624 }
2625 } else { /* bounds command given - calculate approximate pack limits */
2626 for (i = 0; i < 3; i++) {
2627 drvui->frames[drvui->frame_no].cryst_lim[i] =
2628 origin[i] - boxlim[i] / drvui->lat_con[i] - 0.2f;
2629 drvui->frames[drvui->frame_no].cryst_lim[i + 3] =
2630 boxlim[i] / drvui->lat_con[i] - origin[i] + 0.2f;
2631 }
2632 }
2633 if ((Options & L_OPT) != 0)
2634 Calc_Rot (drvui->lookat_v1, drvui->lookat_v2); // Get rotation angles for "lookat"
2635 if (drvui->n_ellips > 1) {
2636 convert_ellipsoid (); /* change Bij and Uij forms to betaij */
2637 fprintf (drvui->flout, "\nAnalysis of Thermal Ellipsoids:\n");
2638 for (i = 1; i < drvui->n_ellips; i++) {
2639 if (drvui->ellips[i].ellips_ismod != 0) {
2640 fprintf (drvui->flout, " (deferred for %c%c%c%c %d due to modulation)\n",
2641 drvui->ellips[i].ellips_l[0], drvui->ellips[i].ellips_l[1],
2642 drvui->ellips[i].ellips_l[2], drvui->ellips[i].ellips_l[3],
2643 drvui->ellips[i].ellips_n);
2644 continue; // modulation must be applied to unchanged Uij values
2645 }
2646 if (eigen
2647 (&biso, drvui->ellips[i].ellips, drvui->ellips[i].ellips_RMS,
2648 drvui->ellips[i].ellips_EV) == 0) {
2649 fprintf (drvui->flout,
2650 "Ellipsoid for %c%c%c%c %d is non-positive definite and has "
2651 "been removed from the active list.\n",
2652 drvui->ellips[i].ellips_l[0], drvui->ellips[i].ellips_l[1],
2653 drvui->ellips[i].ellips_l[2], drvui->ellips[i].ellips_l[3],
2654 drvui->ellips[i].ellips_n);
2655 drvui->ellips[i].ell_type = -100; // disable this one
2656 drvui->ellips[i].ellips_RMS[0] = (float) sqrt (biso / 78.957); /* convert Biso to mu */
2657 } else {
2658 drvui->ellips[i].ell_type += 1;
2659 fprintf (drvui->flout,
2660 "\n%c%c%c%c %3d Equivalent Isotropic B: %6.2f, U: %6.3f\n",
2661 drvui->ellips[i].ellips_l[0], drvui->ellips[i].ellips_l[1],
2662 drvui->ellips[i].ellips_l[2], drvui->ellips[i].ellips_l[3],
2663 drvui->ellips[i].ellips_n, biso, biso / 78.957);
2664 fprintf (drvui->flout, " RMS Amplitudes: %6.3f %6.3f %6.3f\n",
2665 drvui->ellips[i].ellips_RMS[0], drvui->ellips[i].ellips_RMS[1],
2666 drvui->ellips[i].ellips_RMS[2]);
2667 fprintf (drvui->flout, " Eigen Vectors (by row):\n");
2668 for (j = 0; j < 3; j++) {
2669 fprintf (drvui->flout, " ");
2670 for (k = 0; k < 3; k++)
2671 fprintf (drvui->flout, " %8.5f",
2672 drvui->ellips[i].ellips_EV[j][k]);
2673 fprintf (drvui->flout, "\n");
2674 }
2675 }
2676 }
2677 drvui->Ellipsoid_Scale = P_to_C (drvui->Ellipsoid_Prob);
2678 }
2679 } /* end of get_input */
2680
2681 /* ************************************************************** */
2682 /* ************************************************************** */
2683
2684 void
Locate_Triple(void)2685 Locate_Triple (void)
2686
2687 /* routine to try to place the vector triple where it will be seen in the
2688 POV diagram */
2689 {
2690 int i, j;
2691
2692 float temp, pos[3], cpos[3];
2693
2694 for (i = 0; i < 3; ++i) { // convert origin to Cartesian
2695 pos[i] = 0.0f;
2696 for (j = 0; j < 3; ++j)
2697 pos[i] -= (float) drvui->b_mat[i][j] * origin[j];
2698 }
2699 for (i = 0; i <= 2; ++i) { // calculate position of point after POV rotation
2700 cpos[i] = 0.0f;
2701 for (j = 0; j <= 2; ++j)
2702 cpos[i] += (float) G_Rot[j][i] * pos[j];
2703 }
2704 for (i = 0; i < 3; i++) {
2705 temp = 0.04f * (POV_Max[i] - POV_Min[i]);
2706
2707 if (cpos[i] > 0.)
2708 pos[i] = POV_Max[i] + temp;
2709 else
2710 pos[i] = POV_Min[i] - temp;
2711 }
2712 // calculate coordinates of pos before rotation
2713 for (i = 0; i < 3; i++) {
2714 offset[i] = 0.0f;
2715 for (j = 0; j < 3; j++) {
2716 offset[i] += (float) G_Rot[i][j] * pos[j];
2717 }
2718 }
2719 }
2720
2721 /* ************************************************************** */
2722 /* ************************************************************** */
2723
2724 void
make_bmat(int sys,float lat_con[6],double b_mat[3][3],float ginv[3][3],float rec_lat_con[6])2725 make_bmat (int sys, float lat_con[6], double b_mat[3][3], float ginv[3][3],
2726 float rec_lat_con[6])
2727 {
2728 float snal, snbe, snga, csal, csbe, csga, vol;
2729
2730 int i, j;
2731
2732 for (i = 3; i <= 5; ++i)
2733 if (lat_con[i] <= 0)
2734 lat_con[i] = 90.0f;
2735 if (sys == 5) {
2736 lat_con[1] = lat_con[0];
2737 lat_con[5] = 120.0f;
2738 }
2739 if (sys == 4)
2740 lat_con[1] = lat_con[0];
2741 if (sys == 6) {
2742 lat_con[1] = lat_con[0];
2743 lat_con[2] = lat_con[0];
2744 }
2745 for (i = 0; i <= 2; ++i)
2746 for (j = 0; j <= 2; ++j)
2747 b_mat[j][i] = 0.0f; /* initialize matrix */
2748 b_mat[0][0] = 1.0f;
2749 csal = (float) cos (lat_con[3] / RAD);
2750 snal = (float) sin (lat_con[3] / RAD);
2751 csbe = (float) cos (lat_con[4] / RAD);
2752 snbe = (float) sin (lat_con[4] / RAD);
2753 csga = (float) cos (lat_con[5] / RAD);
2754 snga = (float) sin (lat_con[5] / RAD);
2755 for (i = 0; i < 3; i++)
2756 ginv[i][i] = lat_con[i] * lat_con[i];
2757 ginv[0][1] = ginv[1][0] = lat_con[0] * lat_con[1] * csga;
2758 ginv[0][2] = ginv[2][0] = lat_con[0] * lat_con[2] * csbe;
2759 ginv[1][2] = ginv[2][1] = lat_con[1] * lat_con[2] * csal;
2760 b_mat[0][1] = csga;
2761 b_mat[1][1] = snga;
2762 b_mat[0][2] = csbe;
2763 b_mat[1][2] = (csal - csbe * csga) / snga;
2764 b_mat[2][2] = (float) sqrt (1.0 - b_mat[0][2] * b_mat[0][2] -
2765 b_mat[1][2] * b_mat[1][2]);
2766 for (i = 0; i <= 2; ++i)
2767 for (j = 0; j <= 2; ++j)
2768 b_mat[j][i] *= lat_con[i];
2769 /* calculate cell volume and reciprocal lengths */
2770 vol = (float) (b_mat[0][0] * (b_mat[1][1] * b_mat[2][2] - b_mat[1][2] * b_mat[2][1])
2771 - b_mat[0][1] * (b_mat[1][0] * b_mat[2][2] - b_mat[1][2] * b_mat[2][0])
2772 + b_mat[0][2] * (b_mat[1][0] * b_mat[2][1] -
2773 b_mat[1][1] * b_mat[2][0]));
2774 rec_lat_con[0] = lat_con[1] * lat_con[2] * snal / vol;
2775 rec_lat_con[1] = lat_con[0] * lat_con[2] * snbe / vol;
2776 rec_lat_con[2] = lat_con[1] * lat_con[0] * snga / vol;
2777 rec_lat_con[3] = (csbe * csga - csal) / (snbe * snga);
2778 rec_lat_con[4] = (csal * csga - csbe) / (snal * snga);
2779 rec_lat_con[5] = (csal * csbe - csga) / (snal * snbe);
2780 drvui->subsys_ref_volume = 1.0f / vol; /* save the reciprocal volume */
2781
2782 } /* end of make_bmat */
2783
2784 /* ************************************************************** */
2785 /* ************************************************************** */
2786
2787 void
modulate_parameters(float vert[3],double * occ,int sym_no,int atom_no)2788 modulate_parameters (float vert[3], double *occ, int sym_no, int atom_no)
2789 {
2790 /* routine to generate the occupancy and positional parameters for an atom with
2791 "average" positional parameters in 'vert'
2792 This routine relies very heavily on S. van Smaalen, Z. Kristallogr. 219 (2004) 681-691.
2793 A second reference used for 2 and 3D modulation is Jakubowicz et al., Phys. Rev. B 63 (2001). */
2794
2795 double x4[3] = { 0.0, 0.0, 0.0 }, rx4[3] = {
2796 0.0, 0.0, 0.0};
2797 double arg;
2798 float xadd[3] = { 0.0, 0.0, 0.0 };
2799 float rxadd[3];
2800
2801 float eps_inv[3][3];
2802
2803 int id, i1, i2, i3, j, k, theatom;
2804
2805 int axis;
2806
2807 *occ = 1.0;
2808 if (drvui->modulated <= 0)
2809 return;
2810 for (j = 0; j < 3; j++)
2811 for (k = 0; k < 3; k++)
2812 eps_inv[j][k] = (float) drvui->ss_m[sym_no][j][k];
2813 (void) matinv (eps_inv); /* the 'magical' epsilon^(-1) (van Smaalen eq 14) */
2814 for (k = 0; k < drvui->modulated; k++) {
2815 x4[k] = drvui->phaseshift[k];
2816 for (j = 0; j < 3; j++)
2817 x4[k] += drvui->cell_vec[k][j] * vert[j]; /* contains x4=q1.r, x5=q2.r, x6=q3.r */
2818 }
2819 for (k = 0; k < 3; k++) {
2820 rx4[k] = -drvui->ts_m[sym_no][k];
2821 for (j = 0; j < 3; j++) {
2822 rx4[k] += eps_inv[k][j] * x4[j]; /* rx4 is x4 transformed by epsilon^(-1) */
2823 }
2824 rx4[k] = fmod (rx4[k], 1.0);
2825 if (rx4[k] < 0)
2826 rx4[k] += 1.;
2827 }
2828 if (drvui->atoms[atom_no].occ_ismod & 2) { /* check for crenel occupancy modulation */
2829 float lower = drvui->modulate_x[atom_no].atom_occ_crenel[0];
2830
2831 float upper = drvui->modulate_x[atom_no].atom_occ_crenel[1];
2832
2833 if ((rx4[0] < lower) || (rx4[0] > upper))
2834 *occ = 0.0;
2835 if (((upper - 1.0f) > rx4[0]) || ((lower + 1.0f) < rx4[0]))
2836 *occ = 1.0; /* handle the cases when lower < 0 or upper > 1 */
2837 if (*occ == 0.0)
2838 return; /* done if crenel occ is zero */
2839 }
2840
2841 if (drvui->atoms[atom_no].occ_ismod & 1) { /* check for Fourier occupancy modulation */
2842 *occ = drvui->atoms[atom_no].occupancy;
2843 for (j = 0; j < 3; j++)
2844 rx4[j] *= 2.0 * PI;
2845 for (k = 0; k < drvui->no_site_occ; k++) { /* loop through the fourier mod vectors */
2846 id = drvui->modulate_x[k].atom_occpar_id - 1;
2847 theatom = drvui->modulate_x[k].atom_occpar_atom;
2848 if (theatom != atom_no)
2849 continue;
2850 i1 = drvui->modulate_gbl[id].vector_mult[0]; /* multiplier for first cell mod vector */
2851 i2 = drvui->modulate_gbl[id].vector_mult[1]; /* multiplier for second cell mod vector */
2852 i3 = drvui->modulate_gbl[id].vector_mult[2]; /* multiplier for third cell mod vector */
2853 arg = i1 * rx4[0] + i2 * rx4[1] + i3 * rx4[2];
2854 *occ += (float) (drvui->modulate_x[k].atom_occpar[0] * cos (arg)
2855 + drvui->modulate_x[k].atom_occpar[1] * sin (arg));
2856 }
2857 if (*occ < drvui->atoms[atom_no].min_occ) {
2858 *occ = 0.;
2859 return;
2860 }
2861 }
2862
2863 if (drvui->atoms[atom_no].atom_ismod & 1) { /* Fourier position modulation */
2864 if (!(drvui->atoms[atom_no].occ_ismod & 1)) { //only if not already done in previous step
2865 for (j = 0; j < 3; j++)
2866 rx4[j] *= 2.0 * PI;
2867 }
2868 for (j = 0; j < drvui->no_site_displace; j++) { /* loop through the fourier mod vectors */
2869 id = drvui->modulate_3x[j].atom_modpar_id - 1;
2870 axis = drvui->modulate_3x[j].atom_modpar_axis;
2871 theatom = drvui->modulate_3x[j].atom_modpar_atom;
2872 if (theatom != atom_no)
2873 continue;
2874 i1 = drvui->modulate_gbl[id].vector_mult[0]; /* multiplier for first cell mod vector */
2875 i2 = drvui->modulate_gbl[id].vector_mult[1]; /* multiplier for second cell mod vector */
2876 i3 = drvui->modulate_gbl[id].vector_mult[2]; /* multiplier for third cell mod vector */
2877 arg = i1 * rx4[0] + i2 * rx4[1] + i3 * rx4[2];
2878 xadd[axis] += (float) (drvui->modulate_3x[j].atom_modpar[0] * cos (arg)
2879 + drvui->modulate_3x[j].atom_modpar[1] * sin (arg));
2880 }
2881 for (k = 0; k < 3; k++) {
2882 rxadd[k] = 0.0;
2883 for (i1 = 0; i1 < 3; i1++) {
2884 rxadd[k] += (float) drvui->ss[sym_no][k][i1] * xadd[i1];
2885 }
2886 vert[k] += rxadd[k];
2887 }
2888 }
2889
2890 if (drvui->atoms[atom_no].atom_ismod & 2) { /* Position modulation by a sawtooth function */
2891 if ((drvui->atoms[atom_no].occ_ismod & 1) ^ (drvui->atoms[atom_no].atom_ismod & 1)) { // undo any previous multiplication by 2pi
2892 for (j = 0; j < 3; j++)
2893 rx4[j] /= 2.0 * PI;
2894 }
2895 for (k = 0; k < 3; k++)
2896 xadd[k] = 2.0f * drvui->modulate_x[atom_no].atom_mod_sawtooth[k]
2897 * (((float) rx4[k] - drvui->modulate_x[atom_no].atom_mod_sawtooth[3])
2898 / drvui->modulate_x[atom_no].atom_mod_sawtooth[4]);
2899
2900 for (k = 0; k < 3; k++) {
2901 xadd[k] = 0.0;
2902 for (i1 = 0; i1 < 3; i1++) {
2903 rxadd[k] += (float) drvui->ss[sym_no][k][i1] * xadd[i1];
2904 }
2905 vert[k] += rxadd[k];
2906 }
2907 }
2908 }
2909
2910 /* ************************************************************** */
2911 /* ************************************************************** */
2912
2913 int
modulate_uij(float vert[3],int ellips_no,int atom_no,int sym_no,float uij[])2914 modulate_uij (float vert[3], int ellips_no, int atom_no, int sym_no, float uij[])
2915 {
2916 double x4[3] = { 0.0, 0.0, 0.0 }, rx4[3] = {
2917 0.0, 0.0, 0.0};
2918 double arg;
2919
2920 float eps_inv[3][3];
2921
2922 int id, i1, i2, i3, j, k, theatom;
2923
2924 int term;
2925
2926 float biso;
2927
2928 // float uij[6];
2929
2930 for (j = 0; j < 6; j++)
2931 uij[j] = drvui->ellips[ellips_no].ellips[j]; /* copy unmodulated values */
2932 for (k = 0; k < drvui->modulated; k++) {
2933 x4[k] = drvui->phaseshift[k];
2934 for (j = 0; j < 3; j++)
2935 x4[k] += drvui->cell_vec[k][j] * vert[j]; /* contains x4=q1.r, x5=q2.r, x6=q3.r */
2936 }
2937 for (k = 0; k < 3; k++) {
2938 rx4[k] = drvui->ts_m[sym_no][k];
2939 for (j = 0; j < 3; j++) {
2940 eps_inv[j][k] = (float) drvui->ss_m[sym_no][j][k];
2941 rx4[k] += drvui->ss_m[sym_no][k][j] * x4[j]; /* rx4 is x4 transformed by the _mm sym op */
2942 }
2943 rx4[k] = fmod (rx4[k], 1.0);
2944 }
2945 (void) matinv (eps_inv); /* the 'magical' epsilon^(-1) (van Smaalen eq 14) */
2946
2947 for (j = 0; j < 3; j++)
2948 rx4[j] *= 2.0 * PI;
2949 for (k = 0; k < drvui->no_site_U_terms; k++) { /* loop through the fourier mod vectors */
2950 id = drvui->modulate_3t[k].ellips_modpar_id - 1;
2951 theatom = drvui->modulate_x[k].ellips_modpar_atom;
2952 if (theatom != atom_no)
2953 continue;
2954 term = drvui->modulate_3t[k].ellips_modpar_term;
2955 i1 = drvui->modulate_gbl[id].vector_mult[0]; /* multiplier for first cell mod vector */
2956 i2 = drvui->modulate_gbl[id].vector_mult[1]; /* multiplier for second cell mod vector */
2957 i3 = drvui->modulate_gbl[id].vector_mult[2]; /* multiplier for third cell mod vector */
2958 arg = i1 * rx4[0] + i2 * rx4[1] + i3 * rx4[2];
2959 uij[term] += (float) (drvui->modulate_3t[k].ellips_modpar[0] * cos (arg)
2960 + drvui->modulate_3t[k].ellips_modpar[1] * sin (arg));
2961 }
2962
2963 // FIXME? Does rotation by symmetry need to be done before making the modulation correction?
2964
2965 // convert from Uij to betaij
2966
2967 for (j = 0; j < 3; ++j)
2968 uij[j] *= drvui->rec_lat_con[j] * drvui->rec_lat_con[j] * 0.25f;
2969 uij[3] *= drvui->rec_lat_con[0] * drvui->rec_lat_con[1] * 0.25f;
2970 uij[4] *= drvui->rec_lat_con[0] * drvui->rec_lat_con[2] * 0.25f;
2971 uij[5] *= drvui->rec_lat_con[1] * drvui->rec_lat_con[2] * 0.25f;
2972 // multiply Uij by 8pi^2
2973 for (j = 0; j < 6; ++j)
2974 uij[j] *= 78.9568f;
2975
2976 // calculate eigenvalues and eigenvectors
2977
2978 if (eigen
2979 (&biso, uij, drvui->ellips[ellips_no].ellips_RMS,
2980 drvui->ellips[ellips_no].ellips_EV) == 0) {
2981 fprintf (drvui->flout,
2982 "Ellipsoid for %c%c%c%c %d is non-positive definite and has "
2983 "been removed from the active list.\n",
2984 drvui->ellips[ellips_no].ellips_l[0],
2985 drvui->ellips[ellips_no].ellips_l[1],
2986 drvui->ellips[ellips_no].ellips_l[2],
2987 drvui->ellips[ellips_no].ellips_l[3], drvui->ellips[ellips_no].ellips_n);
2988 return (1);
2989 } else {
2990 fprintf (drvui->flout, "\n%c%c%c%c %3d Equivalent Isotropic B: %6.2f, U: %6.3f\n",
2991 drvui->ellips[ellips_no].ellips_l[0],
2992 drvui->ellips[ellips_no].ellips_l[1],
2993 drvui->ellips[ellips_no].ellips_l[2],
2994 drvui->ellips[ellips_no].ellips_l[3], drvui->ellips[ellips_no].ellips_n,
2995 biso, biso / 78.957);
2996 fprintf (drvui->flout, " RMS Amplitudes: %6.3f %6.3f %6.3f\n",
2997 drvui->ellips[ellips_no].ellips_RMS[0],
2998 drvui->ellips[ellips_no].ellips_RMS[1],
2999 drvui->ellips[ellips_no].ellips_RMS[2]);
3000 fprintf (drvui->flout, " Eigen Vectors (by row):\n");
3001 for (j = 0; j < 3; j++) {
3002 fprintf (drvui->flout, " ");
3003 for (k = 0; k < 3; k++)
3004 fprintf (drvui->flout, " %8.5f",
3005 drvui->ellips[ellips_no].ellips_EV[j][k]);
3006 fprintf (drvui->flout, "\n");
3007 }
3008 }
3009 return 0;
3010
3011 }
3012
3013 /* ************************************************************** */
3014 /* ************************************************************** */
3015
3016 float *
polygon_normal_3d(int n,float v[])3017 polygon_normal_3d (int n, float v[])
3018 /* helper function for not_in_slab(), simplified from John Burkardt's geometry.c */
3019 {
3020 int i;
3021
3022 int j;
3023
3024 float *normal;
3025
3026 float normal_norm;
3027
3028 float *p;
3029
3030 float *v1;
3031
3032 float *v2;
3033
3034 normal = new float[3];
3035
3036 v1 = new float[3];
3037
3038 v2 = new float[3];
3039
3040 p = new float[3];
3041
3042 normal[0] = normal[1] = normal[2] = 0.;
3043
3044 for (i = 0; i < 3; i++) {
3045 v1[i] = v[i + 1 * 3] - v[i + 0 * 3];
3046 }
3047
3048 for (j = 2; j < n; j++) {
3049 for (i = 0; i < 3; i++) {
3050 v2[i] = v[i + j * 3] - v[i + 0 * 3];
3051 }
3052
3053 vcross (v1, v2, p);
3054
3055 for (i = 0; i < 3; i++) {
3056 normal[i] = normal[i] + p[i];
3057 }
3058
3059 for (i = 0; i < 3; i++) {
3060 v1[i] = v2[i];
3061 }
3062
3063 }
3064 normal_norm =
3065 (float) sqrt (normal[0] * normal[0] + normal[1] * normal[1] +
3066 normal[2] * normal[2]);
3067
3068 if (normal_norm != 0.0) {
3069 for (i = 0; i < 3; i++) {
3070 normal[i] = normal[i] / normal_norm;
3071 }
3072 }
3073
3074 delete[]v1;
3075 delete[]v2;
3076
3077 return normal;
3078 }
3079
3080 /* ************************************************************** */
3081 /* ************************************************************** */
3082
3083 float
polygon_solid_angle_3d(int n,float v[],float p[3])3084 polygon_solid_angle_3d (int n, float v[], float p[3])
3085 /* helper function for not_in_slab(), simplified from John Burkardt's geometry.c */
3086 {
3087 float a[3];
3088
3089 double angle;
3090
3091 float area = 0.0;
3092
3093 float b[3];
3094
3095 float c[3];
3096
3097 int j;
3098
3099 int jp1;
3100
3101 float normal1[3];
3102
3103 float normal1_norm;
3104
3105 float normal2[3];
3106
3107 float normal2_norm;
3108
3109 float *plane;
3110
3111 float r1[3];
3112
3113 double s;
3114
3115 float value;
3116
3117 if (n < 3)
3118 return 0.0;
3119
3120 plane = polygon_normal_3d (n, v);
3121
3122 a[0] = v[0 + (n - 1) * 3] - v[0 + 0 * 3];
3123 a[1] = v[1 + (n - 1) * 3] - v[1 + 0 * 3];
3124 a[2] = v[2 + (n - 1) * 3] - v[2 + 0 * 3];
3125
3126 for (j = 0; j < n; j++) {
3127 r1[0] = v[0 + j * 3] - p[0];
3128 r1[1] = v[1 + j * 3] - p[1];
3129 r1[2] = v[2 + j * 3] - p[2];
3130
3131 jp1 = j + 1;
3132 if (jp1 > n - 1)
3133 jp1 = 0;
3134
3135 b[0] = v[0 + jp1 * 3] - v[0 + j * 3];
3136 b[1] = v[1 + jp1 * 3] - v[1 + j * 3];
3137 b[2] = v[2 + jp1 * 3] - v[2 + j * 3];
3138
3139 vcross (a, r1, normal1);
3140
3141 normal1_norm =
3142 (float) sqrt (normal1[0] * normal1[0] + normal1[1] * normal1[1] +
3143 normal1[2] * normal1[2]);
3144
3145 vcross (r1, b, normal2);
3146
3147 normal2_norm =
3148 (float) sqrt (normal2[0] * normal2[0] + normal2[1] * normal2[1] +
3149 normal2[2] * normal2[2]);
3150
3151 s = vdot (normal1, normal2)
3152 / (normal1_norm * normal2_norm);
3153
3154 if (s <= -1.) {
3155 angle = PI;
3156 } else if (s >= 1.) {
3157 angle = 0.;
3158 } else
3159 angle = acos (s);
3160
3161 vcross (a, plane, c);
3162 s = vdot (b, c);
3163
3164 if (0.0 < s)
3165 area += (float) (PI - angle);
3166 else
3167 area += (float) (PI + angle);
3168
3169 a[0] = -b[0];
3170 a[1] = -b[1];
3171 a[2] = -b[2];
3172
3173 }
3174 area = area - (float) (PI * (double) (n - 2));
3175
3176 if (0.0 < vdot (plane, r1))
3177 value = -area;
3178 else
3179 value = area;
3180
3181 delete[]plane;
3182
3183 return value;
3184 }
3185
3186 /* ************************************************************** */
3187 /* ************************************************************** */
3188
3189 int
not_in_slab(float x,float y,float z)3190 not_in_slab (float x, float y, float z)
3191 /* adapted from polyhedron_contains_point_3d (Carvalho et al. Graphics Gems V)
3192 as contained in John Burkardt's geometry.c on http://www.csit.fsu.edu/~burkardt/ */
3193 {
3194 float area;
3195
3196 int face;
3197
3198 int i, k;
3199
3200 int node;
3201
3202 int node_num_face = 4;
3203
3204 float v_face[3 * 4];
3205
3206 float p[3];
3207 static const int face_point[4 * 6] =
3208 { 1, 2, 4, 3, 1, 3, 6, 5, 1, 5, 7, 2, 5, 6, 8, 7, 2, 7, 8, 4, 3, 4, 8, 6 };
3209
3210
3211 p[0] = x;
3212 p[1] = y;
3213 p[2] = z;
3214
3215 area = 0.0;
3216 for (face = 0; face < 6; face++) {
3217
3218 for (k = 0; k < 4; k++) {
3219 node = face_point[k + face * 4];
3220 for (i = 0; i < 3; i++) {
3221 v_face[i + k * 3] = slabv[i + (node - 1) * 3];
3222 }
3223 }
3224 area = area + polygon_solid_angle_3d (node_num_face, v_face, p);
3225 }
3226 if (area < -2. * PI || 2. * PI < area)
3227 return 0;
3228 else
3229 return 1;
3230 }
3231
3232 /* ************************************************************** */
3233 /* ************************************************************** */
3234
3235 void
Output_Spheres(float * radii,int i)3236 Output_Spheres (float *radii, int i)
3237
3238 /* routine to output the spheres with vertices defined in s_vert and o_vert
3239 The radius of the sphere is in radii */
3240 {
3241 float glr, glg, glb;
3242
3243 int k;
3244
3245 int outside = 0;
3246
3247 int l1, l2;
3248
3249 int omit;
3250
3251 char sphere_col_p[40];
3252
3253 char sphere_col_v[40];
3254
3255 strcpy (sphere_col_p, drvui->spheres[i].sphere_col);
3256 strcpy (sphere_col_v, sphere_col_p);
3257 Transform_VRML_Color (sphere_col_v);
3258 Transform_POV_Color (sphere_col_p);
3259
3260 for (k = 0; k < nvert; ++k) {
3261 if (clipflag != 0) {
3262 outside = 0;
3263 for (l1 = 0; l1 < 3; l1++)
3264 if (o_vert[3 * k + l1] <
3265 drvui->frames[drvui->frame_no].clip_lim[l1] - 0.01
3266 || o_vert[3 * k + l1] >
3267 drvui->frames[drvui->frame_no].clip_lim[l1 + 3] + 0.01)
3268 outside = 1;
3269 }
3270
3271 omit = 0;
3272 for (l1 = 0; l1 < Omit->nomits; l1++) {
3273 if (Omit->omit1[l1] == i && Omit->omit2[l1] == k)
3274 omit = 1;
3275 }
3276
3277 if (omit == 0 && (clipflag == 0 || outside == 0)) {
3278 float vert1[3], vert2[3];
3279
3280 for (l1 = 0; l1 < 3; l1++)
3281 vert1[l1] = o_vert[3 * k + l1];
3282 for (l1 = 0; l1 < 3; ++l1) { /* convert vertex coordinates to Cartesian */
3283 vert2[l1] = 0.0f;
3284 for (l2 = 0; l2 < 3; ++l2)
3285 vert2[l1] += (float) drvui->b_mat[l1][l2] * (vert1[l2] - origin[l2]);
3286 }
3287 if (doPOV) {
3288 fprintf (drvui->fpoutp,
3289 " /* Sphere for %c%c%c%c at %8.5f %8.5f %8.5f */ \n",
3290 drvui->spheres[i].sphere_l[0], drvui->spheres[i].sphere_l[1],
3291 drvui->spheres[i].sphere_l[2], drvui->spheres[i].sphere_l[3],
3292 vert1[0], vert1[1], vert1[2]);
3293 fprintf (drvui->fpoutp, " sphere{<%8.5f, %8.5f, %8.5f>, %8.5f\n",
3294 vert2[0], vert2[1], vert2[2], radii[k]);
3295 fprintf (drvui->fpoutp, " texture{pigment{color %s }}\n", sphere_col_p);
3296 fprintf (drvui->fpoutp, " finish{phong %5.2f phong_size %5.2f}\n",
3297 drvui->Phong_Value, drvui->Phong_Size);
3298 fprintf (drvui->fpoutp, " }\n");
3299 }
3300 glPushMatrix ();
3301 glLoadName (i);
3302 glPushName (k);
3303 glTranslatef (vert2[0], vert2[1], vert2[2]);
3304 (void) sscanf (sphere_col_v, "%f %f %f", &glr, &glg, &glb);
3305 glColor3f (glr, glg, glb);
3306 glutSolidSphere (radii[k], 10, 10);
3307 glPopName ();
3308 glPopMatrix ();
3309 if (doAsy) {
3310 fprintf (drvui->fpouta,
3311 " // Sphere for %c%c%c%c at %8.5f %8.5f %8.5f \n",
3312 drvui->spheres[i].sphere_l[0], drvui->spheres[i].sphere_l[1],
3313 drvui->spheres[i].sphere_l[2], drvui->spheres[i].sphere_l[3],
3314 vert1[0], vert1[1], vert1[2]);
3315 fprintf (drvui->fpouta, "draw(pic, shift (%8.5f, %8.5f, %8.5f)*scale3(%.2f)*unitsphere,rgb(%4.2f,%4.2f,%4.2f));\n",
3316 vert2[0], vert2[1], vert2[2], radii[k],glr,glg,glb);
3317 }
3318 if (doVrml) {
3319 if (no_comment == 0)
3320 fprintf (drvui->fpoutv,
3321 "# Sphere for %c%c%c%c at %8.5f %8.5f %8.5f \n",
3322 drvui->spheres[i].sphere_l[0], drvui->spheres[i].sphere_l[1],
3323 drvui->spheres[i].sphere_l[2], drvui->spheres[i].sphere_l[3],
3324 vert1[0], vert1[1], vert1[2]);
3325 if (Vrml2) {
3326 fprintf (drvui->fpoutv, " Transform {\n");
3327 fprintf (drvui->fpoutv, " translation %8.5f %8.5f %8.5f\n",
3328 vert2[0], vert2[1], vert2[2]);
3329 fprintf (drvui->fpoutv, " children [\n Shape {\n");
3330 fprintf (drvui->fpoutv,
3331 " appearance Appearance {material Material { diffuseColor %s} }\n",
3332 sphere_col_v);
3333 fprintf (drvui->fpoutv, " geometry Sphere { radius %8.5f}\n",
3334 radii[k]);
3335 fprintf (drvui->fpoutv, " }\n ]\n }\n");
3336 } else { // VRML1
3337 fprintf (drvui->fpoutv, " Separator {\n");
3338 fprintf (drvui->fpoutv, " Transform {");
3339 fprintf (drvui->fpoutv, " translation %8.5f %8.5f %8.5f}\n", vert2[0],
3340 vert2[1], vert2[2]);
3341 fprintf (drvui->fpoutv, " Material {\n");
3342 fprintf (drvui->fpoutv, " shininess 0.3\n");
3343 fprintf (drvui->fpoutv, " diffuseColor %s\n }\n", sphere_col_v);
3344 fprintf (drvui->fpoutv, " Sphere { radius %8.5f}\n }\n", radii[k]);
3345 }
3346 }
3347 }
3348 }
3349 }
3350
3351 /* ************************************************************** */
3352 /* ************************************************************** */
3353
3354 void
plot_vrml_poly(int polyno)3355 plot_vrml_poly (int polyno)
3356 {
3357 int i, j, k, l, m, n;
3358
3359 float *out;
3360
3361 int *point, equal;
3362
3363 char edgecolor[40];
3364
3365 /* get space for unique vertices and a list to point to them */
3366
3367 if (!(out = (float *) zalloc ((unsigned) ((draw_list + 10) * 3 * sizeof (float))))) {
3368 Error_Box ("Unable to allocate space in plot_vrml_poly.");
3369 return;
3370 }
3371 if (!(point = (int *) zalloc ((unsigned) ((draw_list + 10) * sizeof (int))))) {
3372 free (out);
3373 Error_Box ("Unable to allocate space in plot_vrml_poly.");
3374 return;
3375 }
3376 if (!Vrml2)
3377 fprintf (drvui->fpoutv, "Coordinate3 { point [\n");
3378 n = 0;
3379 for (i = 0; i < draw_list; i++) { /* build list of unique points */
3380 j = poly_list[i];
3381 /* fprintf(fpoutv,"# %d-sided polygon:\n",j); */
3382 for (k = i + 1; k <= i + j; k++) {
3383 l = poly_list[k] - 1;
3384 for (m = 0; m < 3; m++)
3385 out[3 * n + m] = s_vert[3 * l + m];
3386 equal = 0;
3387 for (m = 0; m < n; m++) { /* check if new */
3388 if ((fabs (out[3 * m] - out[3 * n]) < 0.0002)
3389 && (fabs (out[3 * m + 1] - out[3 * n + 1]) < 0.0002)
3390 && (fabs (out[3 * m + 2] - out[3 * n + 2]) < 0.0002)) {
3391 point[k] = m; /* point previously found */
3392 equal = 1;
3393 break;
3394 }
3395 }
3396 if (equal == 0)
3397 point[k] = n++; /* new point - save it by incrementing n */
3398 }
3399 i = i + j;
3400 }
3401 --n;
3402 if (Vrml2)
3403 fprintf (drvui->fpoutv,
3404 "}\n geometry IndexedFaceSet { coord Coordinate{ point [\n");
3405 for (l = 0; l < n; l++) /* output the unique set of coordinates */
3406 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f,\n", out[3 * l], out[3 * l + 1],
3407 out[3 * l + 2]);
3408 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f]\n", out[3 * n], out[3 * n + 1],
3409 out[3 * n + 2]);
3410 if (Vrml2) {
3411 fprintf (drvui->fpoutv, " }\n");
3412 fprintf (drvui->fpoutv, " coordIndex [\n");
3413 } else {
3414 fprintf (drvui->fpoutv, "}\n IndexedFaceSet { coordIndex [\n");
3415 }
3416 l = 0;
3417 for (i = 0; i < draw_list - 1; i++) {
3418 j = poly_list[i];
3419 l++;
3420 for (k = 0; k < j; k++)
3421 fprintf (drvui->fpoutv, "%d,", point[l++]);
3422 if (i + j < draw_list - 2) {
3423 fprintf (drvui->fpoutv, "-1,\n");
3424 } else {
3425 if (Vrml2) {
3426 fprintf (drvui->fpoutv, "-1]\n");
3427 } else {
3428 fprintf (drvui->fpoutv, "-1]\n }\n");
3429 }
3430 }
3431 i = i + j;
3432 }
3433 if (Vrml2) {
3434 fprintf (drvui->fpoutv, " solid FALSE\n");
3435 fprintf (drvui->fpoutv, " convex TRUE\n");
3436 fprintf (drvui->fpoutv, " creaseAngle 1.5708\n");
3437 fprintf (drvui->fpoutv, " }\n }\n");
3438 }
3439 if (edges > 0) {
3440 if (drvui->polyhedra[polyno].poly_rad_edge != 0.)
3441 strncpy (edgecolor, drvui->polyhedra[polyno].poly_col_edge, 39);
3442 else
3443 strncpy (edgecolor, drvui->col_edge, 39);
3444 edgecolor[39] = 0;
3445 Transform_VRML_Color (edgecolor);
3446 if (Vrml2) {
3447 fprintf (drvui->fpoutv, " Shape {\n");
3448 fprintf (drvui->fpoutv,
3449 " geometry IndexedLineSet {\n coord Coordinate{\n point [\n");
3450 for (l = 0; l < n; l++) /* output the unique set of coordinates */
3451 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f,\n", out[3 * l],
3452 out[3 * l + 1], out[3 * l + 2]);
3453 fprintf (drvui->fpoutv, " %5.3f %5.3f %5.3f]\n", out[3 * n],
3454 out[3 * n + 1], out[3 * n + 2]);
3455 fprintf (drvui->fpoutv, " }\n coordIndex [\n ");
3456 } else {
3457 fprintf (drvui->fpoutv, " Material {\n diffuseColor %s \n }\n", edgecolor);
3458 fprintf (drvui->fpoutv, "IndexedLineSet { coordIndex [\n ");
3459 }
3460 l = 0;
3461 for (i = 0; i < draw_list - 1; i++) {
3462 j = poly_list[i];
3463 l++;
3464 for (k = 0; k < j; k++)
3465 fprintf (drvui->fpoutv, "%d,", point[l++]);
3466 if (i + j < draw_list - 2) {
3467 fprintf (drvui->fpoutv, "%d,-1,\n ", point[l - j]);
3468 } else {
3469 if (Vrml2) {
3470 fprintf (drvui->fpoutv, "%d,-1]\n ", point[l - j]);
3471 } else {
3472 fprintf (drvui->fpoutv, "%d,-1]\n }\n }\n", point[l - j]);
3473 }
3474 }
3475 i = i + j;
3476 }
3477 if (Vrml2) {
3478 fprintf (drvui->fpoutv,
3479 " color Color { color [%s]}\n colorIndex[0]\n colorPerVertex FALSE}\n",
3480 edgecolor);
3481 fprintf (drvui->fpoutv, " }\n");
3482 }
3483 } else {
3484 if (!Vrml2)
3485 fprintf (drvui->fpoutv, " }\n");
3486 }
3487 free (out); /* release the space */
3488 free (point);
3489 return;
3490 }
3491
3492 /* ************************************************************** */
3493 /* ************************************************************** */
3494
3495 void
print_sym(void)3496 print_sym (void)
3497
3498 /* print symmetry information */
3499 {
3500 int i, j, k;
3501
3502 switch (drvui->sys) {
3503 case 1:
3504 fprintf (drvui->flout, "Triclinic - ");
3505 break;
3506 case 2:
3507 fprintf (drvui->flout, "Monoclinic - ");
3508 break;
3509 case 3:
3510 fprintf (drvui->flout, "Orthorhombic - ");
3511 break;
3512 case 4:
3513 fprintf (drvui->flout, "Tetragonal - ");
3514 break;
3515 case 5:
3516 fprintf (drvui->flout, "Hexagonal - ");
3517 break;
3518 case 6:
3519 fprintf (drvui->flout, "Cubic - ");
3520 }
3521 if (!drvui->acentric) {
3522 fprintf (drvui->flout, "Centric - ");
3523 } else {
3524 fprintf (drvui->flout, "Acentric - ");
3525 }
3526 switch (drvui->nbr) {
3527 case 1:
3528 fprintf (drvui->flout, "P lattice");
3529 break;
3530 case 2:
3531 fprintf (drvui->flout, "A-centered");
3532 break;
3533 case 3:
3534 fprintf (drvui->flout, "B-centered");
3535 break;
3536 case 4:
3537 fprintf (drvui->flout, "C-centered");
3538 break;
3539 case 5:
3540 fprintf (drvui->flout, "F-centered");
3541 break;
3542 case 6:
3543 fprintf (drvui->flout, "I-centered");
3544 break;
3545 case 7:
3546 fprintf (drvui->flout, "R-centered");
3547 break;
3548 default:
3549 fprintf (drvui->flout, "**** Unknown Lattice ***");
3550 }
3551 fprintf (drvui->flout, " with %d Symmetry Operators\n\n", drvui->ng);
3552 for (i = 1; i <= drvui->ng; ++i) {
3553
3554 for (j = 1; j <= 3; ++j) {
3555 for (k = 1; k <= 3; ++k) {
3556 drvui->rss[i - 1][j - 1][k - 1] = (float) drvui->ss[i - 1][j - 1][k - 1];
3557 if (drvui->ss[i - 1][j - 1][k - 1] != 0) {
3558 if (drvui->ss[i - 1][j - 1][k - 1] > 0) {
3559 fprintf (drvui->flout, "+");
3560 } else {
3561 fprintf (drvui->flout, "-");
3562 }
3563 switch (k) {
3564 case 1:
3565 fprintf (drvui->flout, "x");
3566 break;
3567 case 2:
3568 fprintf (drvui->flout, "y");
3569 break;
3570 case 3:
3571 fprintf (drvui->flout, "z");
3572 }
3573 }
3574 }
3575 if (drvui->ts[i - 1][j - 1] != 0.0) {
3576 if (drvui->ts[i - 1][j - 1] == 0.5)
3577 fprintf (drvui->flout, "+1/2");
3578 if (drvui->ts[i - 1][j - 1] == 0.25)
3579 fprintf (drvui->flout, "+1/4");
3580 if (drvui->ts[i - 1][j - 1] == 0.75)
3581 fprintf (drvui->flout, "+3/4");
3582 if (fabs (drvui->ts[i - 1][j - 1] - 0.33333333) < 0.00001)
3583 fprintf (drvui->flout, "+1/3");
3584 if (fabs (drvui->ts[i - 1][j - 1] - 0.66666667) < 0.00001)
3585 fprintf (drvui->flout, "+2/3");
3586 if (fabs (drvui->ts[i - 1][j - 1] - 0.16666667) < 0.00001)
3587 fprintf (drvui->flout, "+1/6");
3588 if (fabs (drvui->ts[i - 1][j - 1] - 0.83333333) < 0.00001)
3589 fprintf (drvui->flout, "+5/6");
3590 }
3591 if (j != 3)
3592 fprintf (drvui->flout, ",");
3593 }
3594 if (i != drvui->ng) {
3595 if (i % 4 != 0) {
3596 fprintf (drvui->flout, " ; ");
3597 } else {
3598 fprintf (drvui->flout, "\n");
3599 }
3600 }
3601 }
3602 Conv_Sym_Mat ();
3603 fprintf (drvui->flout, "\n");
3604 } /* end of print_sym */
3605