1 /*This file is a part of luscus project*/
2 /*Licensed under the Academic Free License version 3.0*/
3 #include<stdio.h>
4 #include<stdlib.h>
5 #include<string.h>
6 #ifdef WINDOWS
7 #include<strings.h>
8 #endif
9 #include<math.h>
10 #include<sys/types.h>
11 #include<sys/stat.h>
12 #include<unistd.h>
13 #include"luscus.h"
14 #include"gv.h"
15 #include"gv_functions.h"
16 #define MAXCHAR 127
17 #define MIN(x,y) (x<y? x : y)
18 #define MAX(x,y) (x>y? x : y)
19 #define R529 0.52917721067e0
20
21 typedef unsigned char byte_t;
22 MOL *m;
23 MOL *mol;
24 int number_of_elements;
25 ELEM_DATA *e;
26 INPUT_DATA Input_Data;
27 char *input_file_name = NULL;
28 char *input_file_type = NULL; /*19.02.2013 added in order to discriminate file types*/
29 char *current_directory = NULL;
30 int should_bonding_be_determined;
31 static int isBinary;
32 static int isPacked;
33 static int Fmarker = 0;
34 static time_t init_file_time = 0;
35 static long init_file_ntime = 0;
36
37 char file_open_success = 0;
38
39 static int ReadBlock(double*, int, FILE*);
40
41 /* Packing parameters */
42 struct packing_params_t
43 {
44 int range; /* really, half range */
45 int nbytes; /* bytes per packed value */
46
47 double xlimit[4];
48 double xdelta[3];
49 int ylimit[4];
50 int ydelta[3];
51
52 int ileft, iright;
53 };
54
55 /*definition of function*/
56 static int getPackingParams(const char*, struct packing_params_t*);
57 static void unpackData(double*, byte_t*, unsigned int, int, const struct packing_params_t*);
58
deallocate_mol_data(void)59 void deallocate_mol_data(void)
60 {
61 int i;
62 #ifdef EBUG
63 printf("DEALLOCATE MOL DATA\n"); fflush(stdout);
64 #endif
65 deallocate_msrf();
66 for (i = 0; i < n_geometries; i++)
67 deallocate_data(&mol[i]);
68
69 free(mol);
70 n_geometries = 0;
71
72 if (geo_file_position) free(geo_file_position);
73 geo_file_position = NULL;
74 m = NULL;
75
76 remove_all_watched_coords();
77 }
78
deallocate_data(MOL * mp)79 void deallocate_data(MOL *mp)
80 {
81 #ifdef EBUG
82 printf("DEALLOCATE DATA\n"); fflush(stdout);
83 #endif
84
85 deallocate_atoms(mp);
86 deallocate_bonds(mp);
87 deallocate_dipole(mp);
88 deallocate_grids(mp);
89 deallocate_vibrations(mp);
90 deallocate_vectors(mp);
91 deallocate_triangles(mp);
92 deallocate_spheres(mp);
93 deallocate_surfaces(mp);
94 deallocate_cells(mp);
95 deallocate_textboxes(mp);
96 unmark_all();
97 unselect_all();
98 }
99
deallocate_atoms(MOL * mp)100 void deallocate_atoms(MOL *mp)
101 {
102 int i;
103 if (!mp) return;
104
105 if (mp->xyz) free(mp->xyz);
106 mp->xyz = NULL;
107
108 if (mp->charge_m) free(mp->charge_m);
109 mp->charge_m = NULL;
110
111 if (mp->charge_l) free(mp->charge_l);
112 mp->charge_l = NULL;
113
114 if (mp->additional_numeration) free(mp->additional_numeration);
115 mp->additional_numeration = NULL;
116
117 if (mp->symmetry) free(mp->symmetry);
118 mp->symmetry = NULL;
119
120 if (mp->name)
121 {
122 for(i = 0; i < mp->natom; i++)
123 if (mp->name[i])
124 free(mp->name[i]);
125 free(mp->name);
126 }
127 mp->name = NULL;
128
129 if (mp->pixdata)
130 {
131 for(i = 0; i < mp->natom; i++)
132 {
133 if (mp->pixdata[i].pix_index.pixels) free(mp->pixdata[i].pix_index.pixels);
134 if (mp->pixdata[i].pix_addnum.pixels) free(mp->pixdata[i].pix_addnum.pixels);
135 if (mp->pixdata[i].pix_symm.pixels) free(mp->pixdata[i].pix_symm.pixels);
136 if (mp->pixdata[i].pix_name.pixels) free(mp->pixdata[i].pix_name.pixels);
137 if (mp->pixdata[i].pix_charge_m.pixels) free(mp->pixdata[i].pix_charge_m.pixels);
138 if (mp->pixdata[i].pix_charge_l.pixels) free(mp->pixdata[i].pix_charge_l.pixels);
139 }
140 free(mp->pixdata);
141 }
142
143 if (mp->elem)
144 {
145 for(i = 0; i < mp->natom; i++)
146 if (mp->elem[i].name)
147 free(mp->elem[i].name);
148 free(mp->elem);
149 }
150 mp->elem = NULL;
151
152 mp->n_selected = 0;
153 mp->ishow = 0;
154 mp->natom = 0;
155
156 /* if (m->ishow & HAS_ATOMNAMES) m->ishow ^= HAS_ATOMNAMES;
157 if (m->ishow & HAS_MULLIKEN) m->ishow ^= HAS_MULLIKEN;
158 if (m->ishow & HAS_LOPROP) m->ishow ^= HAS_LOPROP;
159 if (m->ishow & HAS_ATOMNUMS) m->ishow ^= HAS_ATOMNUMS;
160 if (m->ishow & HAS_ATOMSYMM) m->ishow ^= HAS_ATOMSYMM;*/
161
162 /*deallocate atom list*/
163 #ifdef EBUG
164 printf("deallocate atom list\n"); fflush(stdout);
165 #endif
166 deallocate_atom_list();
167 }
168
deallocate_bonds(MOL * mp)169 void deallocate_bonds(MOL *mp)
170 {
171 #ifdef EBUG
172 printf("DEALLOCATE BONDS\n"); fflush(stdout);
173 #endif
174 if (mp->bond) free(mp->bond);
175 mp->bond = NULL;
176 mp->nbond = 0;
177 }
178
deallocate_dipole(MOL * mp)179 void deallocate_dipole(MOL *mp)
180 {
181 #ifdef EBUG
182 printf("DEALLOCATE DIPOLE\n"); fflush(stdout);
183 #endif
184 /*no dealocations, just unmark dipole for drawing!*/
185 if (mp->ishow & HAS_DIPOLE) mp->ishow ^= HAS_DIPOLE;
186 }
187
deallocate_grids(MOL * mp)188 void deallocate_grids(MOL *mp)
189 {
190 #ifdef EBUG
191 printf("DEALLOCATE GRIDS\n"); fflush(stdout);
192 #endif
193 if (mp->grid_type) free(mp->grid_type);
194 mp->grid_type = NULL;
195 #ifdef EBUG
196 printf("TU SAM 000\n"); fflush(stdout);
197 #endif
198 if (mp->grid_energy) free(mp->grid_energy);
199 mp->grid_energy = NULL;
200
201 #ifdef EBUG
202 printf("TU SAM 001\n"); fflush(stdout);
203 #endif
204 if (mp->grid_occ) free(mp->grid_occ);
205 mp->grid_occ = NULL;
206 #ifdef EBUG
207 printf("TU SAM 002\n"); fflush(stdout);
208 #endif
209 if (mp->grid_symmetry) free(mp->grid_symmetry);
210 mp->grid_symmetry = NULL;
211 #ifdef EBUG
212 printf("TU SAM 003\n"); fflush(stdout);
213 #endif
214 if (mp->grid_index) free(mp->grid_index);
215 mp->grid_index = NULL;
216 #ifdef EBUG
217 printf("TU SAM 004\n"); fflush(stdout);
218 #endif
219 if (mp->orbital_type) free(mp->orbital_type);
220 mp->orbital_type = NULL;
221 #ifdef EBUG
222 printf("TU SAM 005\n"); fflush(stdout);
223 #endif
224 if (mp->edited) free(mp->edited);
225 mp->edited = NULL;
226 #ifdef EBUG
227 printf("TU SAM 006\n"); fflush(stdout);
228 #endif
229 if (mp->grid.ByteCutOff) free(mp->grid.ByteCutOff);
230 mp->grid.ByteCutOff = NULL;
231 #ifdef EBUG
232 printf("TU SAM 007\n"); fflush(stdout);
233 printf("DELETING GRID.POS POINTER: %ld\n", (long) mp->grid.pos);
234 #endif
235 if (mp->grid.pos) free(mp->grid.pos); /*ERROR!!!*/
236 mp->grid.pos = NULL;
237 #ifdef EBUG
238 printf("TU SAM 008\n"); fflush(stdout);
239 #endif
240 if (mp->grid.title) free(mp->grid.title);
241 mp->grid.title = NULL;
242 #ifdef EBUG
243 printf("TU SAM 009\n"); fflush(stdout);
244 #endif
245 if (mp->grid.titlesArr) free(mp->grid.titlesArr);
246 mp->grid.titlesArr = NULL;
247 #ifdef EBUG
248 printf("TU SAM 010\n"); fflush(stdout);
249 #endif
250 if (mp->grid.FltTitle) free(mp->grid.FltTitle);
251 mp->grid.FltTitle = NULL;
252 #ifdef EBUG
253 printf("TU SAM 011\n"); fflush(stdout);
254 #endif
255 if (mp->grid.iType) free(mp->grid.iType);
256 mp->grid.iType = NULL;
257 #ifdef EBUG
258 printf("TU SAM 012\n"); fflush(stdout);
259 #endif
260 if (mp->grid.dependence) free(mp->grid.dependence);
261 mp->grid.dependence = NULL;
262 #ifdef EBUG
263 printf("TU SAM 013\n"); fflush(stdout);
264 #endif
265 if (mp->grid.values) free(mp->grid.values);
266 mp->grid.values = NULL;
267 #ifdef EBUG
268 printf("TU SAM 014\n"); fflush(stdout);
269 #endif
270
271 mp->ngrids = 0;
272 mp->grid.ngrid = 0;
273 mp->grid.npoints = 0;
274 if (mp->epot) free(mp->epot);
275 mp->epot = NULL;
276 }
277
deallocate_vibrations(MOL * mp)278 void deallocate_vibrations(MOL *mp)
279 {
280 int i;
281 #ifdef EBUG
282 printf("DEALLOCATE VIBRATIONS\n"); fflush(stdout);
283 #endif
284 if (mp->freq) free(mp->freq);
285 mp->freq = NULL;
286 if (mp->ir_intensity) free(mp->ir_intensity);
287 mp->ir_intensity = NULL;
288 if (mp->raman_intensity) free(mp->raman_intensity);
289 mp->raman_intensity = NULL;
290 if (mp->normal_mode)
291 {
292 for(i = 0; i < mp->nvibr; i++) free(mp->normal_mode[i]);
293 free(mp->normal_mode);
294 }
295 mp->normal_mode = NULL;
296 mp->nvibr = 0;
297 if (mp->ishow & HAS_IRINT) mp->ishow ^= HAS_IRINT;
298 if (mp->ishow & HAS_RAMAN) mp->ishow ^= HAS_RAMAN;
299 }
300
deallocate_vectors(MOL * mp)301 void deallocate_vectors(MOL *mp)
302 {
303 #ifdef EBUG
304 printf("DEALLOCATE VECTORS\n"); fflush(stdout);
305 #endif
306 if (mp->vector1) free(mp->vector1);
307 if (mp->vector2) free(mp->vector2);
308 if (mp->radius) free(mp->radius);
309 if (mp->sharpness) free(mp->sharpness);
310 if (mp->vector_color) free(mp->vector_color);
311 mp->vector1 = NULL;
312 mp->vector2 = NULL;
313 mp->radius = NULL;
314 mp->sharpness = NULL;
315 mp->vector_color = NULL;
316 mp->nvector = 0;
317 }
318
deallocate_triangles(MOL * mp)319 void deallocate_triangles(MOL *mp)
320 {
321 #ifdef EBUG
322 printf("DEALLOCATE TRIANGLES\n"); fflush(stdout);
323 #endif
324 if (mp->triangle1) free(mp->triangle1);
325 if (mp->triangle2) free(mp->triangle2);
326 if (mp->triangle3) free(mp->triangle3);
327 if (mp->triangle_color) free(mp->triangle_color);
328
329 mp->triangle1 = NULL;
330 mp->triangle2 = NULL;
331 mp->triangle3 = NULL;
332 mp->triangle_color = NULL;
333 mp->ntriangle = 0;
334 }
335
deallocate_spheres(MOL * mp)336 void deallocate_spheres(MOL *mp)
337 {
338 #ifdef EBUG
339 printf("DEALLOCATE SPHERES\n"); fflush(stdout);
340 #endif
341 if (mp->sphere_center) free(mp->sphere_center);
342 if (mp->sphere_radius) free(mp->sphere_radius);
343 if (mp->sphere_color) free(mp->sphere_color);
344
345 mp->sphere_center = NULL;
346 mp->sphere_radius = NULL;
347 mp->sphere_color = NULL;
348 mp->nsphere = 0;
349 }
350
deallocate_surfaces(MOL * mp)351 void deallocate_surfaces(MOL *mp)
352 {
353 #ifdef EBUG
354 printf("DEALLOCATE SURFACES\n"); fflush(stdout);
355 #endif
356 if (mp->surf1) free(mp->surf1);
357 if (mp->surf2) free(mp->surf2);
358 if (mp->surf3) free(mp->surf3);
359 if (mp->surf_color) free(mp->surf_color);
360
361 mp->surf1 = NULL;
362 mp->surf2 = NULL;
363 mp->surf3 = NULL;
364 mp->surf_color = NULL;
365 mp->nsurf = 0;
366 }
367
deallocate_cells(MOL * mp)368 void deallocate_cells(MOL *mp)
369 {
370 #ifdef EBUG
371 printf("DEALLOCATE CELLS\n"); fflush(stdout);
372 #endif
373 if (mp->cell1) free(mp->cell1);
374 if (mp->cell2) free(mp->cell2);
375 if (mp->cell3) free(mp->cell3);
376 if (mp->cell4) free(mp->cell4);
377 if (mp->cell_color) free(mp->cell_color);
378
379 mp->cell1 = NULL;
380 mp->cell2 = NULL;
381 mp->cell3 = NULL;
382 mp->cell4 = NULL;
383 mp->cell_color = NULL;
384 mp->ncells = 0;
385 }
386
deallocate_textboxes(MOL * mp)387 void deallocate_textboxes(MOL *mp)
388 {
389 int i;
390 #ifdef EBUG
391 printf("DEALLOCATE TEXTBOXES\n"); fflush(stdout);
392 #endif
393 if (!mp) return;
394 for(i = 0; i < mp->ntextboxes; i++)
395 {
396 if (mp->textboxes[i].message) free(mp->textboxes[i].message);
397 if (mp->textboxes[i].font) free(mp->textboxes[i].font);
398 }
399 if (mp->textboxes) free(mp->textboxes);
400 mp->ntextboxes = 0;
401 }
402
new_mol(int number_of_geometries)403 MOL *new_mol(int number_of_geometries)
404 {
405 int i;
406 MOL *mp;
407 #ifdef EBUG
408 printf("NEW MOL\n"); fflush(stdout);
409 #endif
410 n_geometries = number_of_geometries;
411 igeo = 0;
412 ivib = 0;
413 iorb = 0;
414
415 if (number_of_geometries <= 0)
416 {
417 n_geometries = 0;
418 return NULL;
419 }
420 mp = malloc(sizeof(MOL) * number_of_geometries);
421
422 mp->editable = 1;
423
424 for(i = 0; i < number_of_geometries; i++)
425 {
426 m = mp + i;
427
428 m->ishow = 0;
429
430 m->natom = 0;
431 m->xyz = NULL;
432 m->elem = NULL;
433
434 m->charge_m = NULL;
435 m->charge_l = NULL;
436 m->additional_numeration = NULL;
437 m->symmetry = NULL;
438 m->name = NULL;
439
440 m->pixdata = NULL;
441
442 m->n_selected = 0;
443 m->selected[0] = -1;
444 m->selected[1] = -1;
445 m->selected[2] = -1;
446 m->selected[3] = -1;
447 m->n_marked = 0;
448 m->marked = NULL;
449
450 m->nbond = 0;
451 m->bond = NULL;
452
453 m->dipole[0] = 0.0;
454 m->dipole[1] = 0.0;
455 m->dipole[2] = 0.0;
456
457 m->ngrids = 0;
458 m->grid_type = NULL;
459 m->grid_energy = NULL;
460 m->grid_occ = NULL;
461 m->grid_symmetry = NULL;
462 m->grid_index = NULL;
463 m->orbital_type = NULL;
464 m->edited = NULL;
465
466 m->orb_starting_position = 0L;
467 /*long long int*/
468 /* m->orb_starting_position = 0LL;*/
469
470 m->grid.ByteCutOff = NULL;
471 m->grid.pos = NULL;
472 m->grid.title = NULL;
473 m->grid.titlesArr = NULL;
474 m->grid.FltTitle = NULL;
475 m->grid.iType = NULL;
476 m->grid.dependence = NULL;
477 m->grid.values = NULL;
478 m->epot = NULL;
479
480 m->geo_energy = 0.0;
481 m->max_grad = 0.0;
482 m->rms_grad = 0.0;
483
484 m->nvibr = 0;
485 m->freq = NULL;
486 m->ir_intensity = NULL;
487 m->raman_intensity = NULL;
488 m->normal_mode = NULL;
489
490 m->nvector = 0;
491 m->vector1 = NULL;
492 m->vector2 = NULL;
493 m->radius = NULL;
494 m->sharpness = NULL;
495 m->vector_color = NULL;
496
497 m->ntriangle = 0;
498 m->triangle1 = NULL;
499 m->triangle2 = NULL;
500 m->triangle3 = NULL;
501 m->triangle_color = NULL;
502
503 m->nsphere = 0;
504 m->sphere_center = NULL;
505 m->sphere_radius = NULL;
506 m->sphere_color = NULL;
507
508 m->nsurf = 0;
509 m->surf1 = NULL;
510 m->surf2 = NULL;
511 m->surf3 = NULL;
512 m->surf_color = NULL;
513
514 m->ncells = 0;
515 m->cell1 = NULL;
516 m->cell2 = NULL;
517 m->cell3 = NULL;
518 m->cell4 = NULL;
519 m->cell_color = NULL;
520
521 m->ntextboxes = 0;
522 m->textboxes = NULL;
523 }
524 m = mp;
525
526 /*init atom list*/
527 init_atom_list();
528 return mp;
529 }
530
allocate_atoms(MOL * mp,int nnatom)531 void allocate_atoms(MOL *mp, int nnatom)
532 {
533 int i;
534 int old_natom = mp->natom;
535 #ifdef EBUG
536 printf("allocating atoms\n"); fflush(stdout);
537 #endif
538 mp->natom = nnatom;
539 mp->xyz = (XYZ*) realloc(mp->xyz, sizeof(XYZ) * nnatom);
540 mp->charge_m = (double*) realloc(mp->charge_m, sizeof(double) * nnatom);
541 mp->charge_l = (double*) realloc(mp->charge_l, sizeof(double) * nnatom);
542 mp->additional_numeration = (int*) realloc(mp->additional_numeration, sizeof(int) * nnatom);
543 mp->symmetry = (int*) realloc(mp->symmetry, sizeof(int) * nnatom);
544 mp->name = (char**) realloc(mp->name, sizeof(char*) * nnatom);
545 mp->elem = (ELEM_DATA*) realloc(mp->elem, sizeof(ELEM_DATA) * nnatom);
546 mp->pixdata = (PIXDATA*) realloc(mp->pixdata, sizeof(PIXDATA) * nnatom);
547
548 for(i = old_natom; i < nnatom; i++)
549 {
550 mp->elem[i].vdw_rad = 0.0;
551 mp->elem[i].bond_rad = 0.0;
552 mp->elem[i].valency = 0;
553 mp->elem[i].name = NULL;
554
555 mp->charge_m[i] = 0.0;
556 mp->charge_l[i] = 0.0;
557 mp->additional_numeration[i] = 0;
558 mp->symmetry[i] = 0;
559 mp->name[i] = NULL;
560
561 mp->pixdata[i].pix_index.width = 0;
562 mp->pixdata[i].pix_index.height = 0;
563 mp->pixdata[i].pix_index.pixels = NULL;
564 mp->pixdata[i].pix_addnum.width = 0;
565 mp->pixdata[i].pix_addnum.height = 0;
566 mp->pixdata[i].pix_addnum.pixels = NULL;
567 mp->pixdata[i].pix_symm.width = 0;
568 mp->pixdata[i].pix_symm.height = 0;
569 mp->pixdata[i].pix_symm.pixels = NULL;
570 mp->pixdata[i].pix_name.width = 0;
571 mp->pixdata[i].pix_name.height = 0;
572 mp->pixdata[i].pix_name.pixels = NULL;
573 mp->pixdata[i].pix_charge_m.width = 0;
574 mp->pixdata[i].pix_charge_m.height = 0;
575 mp->pixdata[i].pix_charge_m.pixels = NULL;
576 mp->pixdata[i].pix_charge_l.width = 0;
577 mp->pixdata[i].pix_charge_l.height = 0;
578 mp->pixdata[i].pix_charge_l.pixels = NULL;
579 }
580 #ifdef EBUG
581 printf("allocating atoms done\n"); fflush(stdout);
582 #endif
583 }
584
allocate_bonds(MOL * mp,int nnbonds)585 void allocate_bonds(MOL *mp, int nnbonds)
586 {
587 mp->nbond = nnbonds;
588 mp->bond = (BOND*) realloc(mp->bond, sizeof(BOND) * (nnbonds + 1));
589 }
590
allocate_grids(MOL * mp,int nngrid)591 void allocate_grids(MOL *mp, int nngrid)
592 {
593 mp->grid_type = (int*) realloc(mp->grid_type, sizeof(int) * nngrid);
594 mp->grid_energy = (double*) realloc(mp->grid_energy, sizeof(double) * nngrid);
595 mp->grid_occ = (double*) realloc(mp->grid_occ, sizeof(double) * nngrid);
596 mp->grid_symmetry = (int*) realloc(mp->grid_symmetry, sizeof(int) * nngrid);
597 mp->grid_index = (int*) realloc(mp->grid_index, sizeof(int) * nngrid);
598 mp->orbital_type = (int*) realloc(mp->orbital_type, sizeof(int) * nngrid);
599 mp->edited = (int*) realloc(mp->edited, sizeof(int) * nngrid);
600
601 mp->grid.ByteCutOff = (int*) malloc(sizeof(int) * mp->grid.npoints);
602
603 mp->grid.pos = (long*) malloc(sizeof(long) * (unsigned long) mp->grid.nBlock * (unsigned long) mp->grid.ngrid);
604 #ifdef EBUG
605 printf("ALLOCATED GRID.POS: %ld\n", (long) mp->grid.pos); fflush(stdout);
606 #endif
607 /* mp->grid.pos = (long long int*) malloc(sizeof(long long int) * (unsigned int) mp->grid.nBlock * mp->grid.ngrid);*/
608
609 /* mp->grid.title = NULL;*/
610 mp->grid.titlesArr = (char**) malloc(sizeof(char*) * nngrid);
611 mp->grid.FltTitle = (int*) malloc(sizeof(int) * nngrid);
612 mp->grid.iType = (char*) malloc(sizeof(char) * nngrid);
613 mp->grid.dependence = (int*) malloc(sizeof(int) * nngrid);
614 mp->grid.values = (double*) malloc(sizeof(double) * mp->grid.npoints);
615 }
616
allocate_vibs(MOL * mp,int nnvibs)617 void allocate_vibs(MOL *mp, int nnvibs)
618 {
619 mp->nvibr = nnvibs;
620 mp->freq = (double*) realloc(mp->freq, sizeof(double) * nnvibs);
621 mp->ir_intensity = (double*) realloc(mp->ir_intensity, sizeof(double) * nnvibs);
622 mp->raman_intensity = (double*) realloc(mp->raman_intensity, sizeof(double) * nnvibs);
623 mp->normal_mode = (XYZ**) realloc(mp->normal_mode, sizeof(XYZ*) * nnvibs);
624 mp->normal_mode[nnvibs-1] = (XYZ*) malloc(sizeof(XYZ) * mp->natom);
625 }
626
allocate_vectors(MOL * mp,int nnvec)627 void allocate_vectors(MOL *mp, int nnvec)
628 {
629 mp->nvector = nnvec;
630 mp->vector1 = (XYZ*) realloc(mp->vector1, sizeof(XYZ) * nnvec);
631 mp->vector2 = (XYZ*) realloc(mp->vector2, sizeof(XYZ) * nnvec);
632 mp->radius = (double*) realloc(mp->radius, sizeof(double) * nnvec);
633 mp->sharpness = (double*) realloc(mp->sharpness, sizeof(double) * nnvec);
634 mp->vector_color = (color_t*) realloc(mp->vector_color, sizeof(color_t) * nnvec);
635 }
636
allocate_triangles(MOL * mp,int nntri)637 void allocate_triangles(MOL *mp, int nntri)
638 {
639 mp->ntriangle = nntri;
640 mp->triangle1 = (XYZ*) realloc(mp->triangle1, sizeof(XYZ) * nntri);
641 mp->triangle2 = (XYZ*) realloc(mp->triangle2, sizeof(XYZ) * nntri);
642 mp->triangle3 = (XYZ*) realloc(mp->triangle3, sizeof(XYZ) * nntri);
643 mp->triangle_color = (color_t*) realloc(mp->triangle_color, sizeof(color_t) * nntri);
644 }
645
allocate_spheres(MOL * mp,int nnspher)646 void allocate_spheres(MOL *mp, int nnspher)
647 {
648 mp->nsphere = nnspher;
649 mp->sphere_center = (XYZ*) realloc(mp->sphere_center, sizeof(XYZ) * nnspher);
650 mp->sphere_radius = (double*) realloc(mp->sphere_radius, sizeof(double) * nnspher);
651 mp->sphere_color = (color_t*) realloc(mp->sphere_color, sizeof(color_t) * nnspher);
652 }
653
allocate_surfaces(MOL * mp,int nnsurf)654 void allocate_surfaces(MOL *mp, int nnsurf)
655 {
656 mp->nsurf = nnsurf;
657 mp->surf1 = (XYZ*) realloc(mp->surf1, sizeof(XYZ) * nnsurf);
658 mp->surf2 = (XYZ*) realloc(mp->surf2, sizeof(XYZ) * nnsurf);
659 mp->surf3 = (XYZ*) realloc(mp->surf3, sizeof(XYZ) * nnsurf);
660 mp->surf_color = (color_t*) realloc(mp->surf_color, sizeof(color_t) * nnsurf);
661 }
662
allocate_cells(MOL * mp,int nncells)663 void allocate_cells(MOL *mp, int nncells)
664 {
665 mp->ncells = nncells;
666 mp->cell1 = (XYZ*) realloc(mp->cell1, sizeof(XYZ) * nncells);
667 mp->cell2 = (XYZ*) realloc(mp->cell2, sizeof(XYZ) * nncells);
668 mp->cell3 = (XYZ*) realloc(mp->cell3, sizeof(XYZ) * nncells);
669 mp->cell4 = (XYZ*) realloc(mp->cell4, sizeof(XYZ) * nncells);
670 mp->cell_color = (color_t*) realloc(mp->cell_color, sizeof(color_t) * nncells);
671 }
672
allocate_textbox(MOL * mp)673 void allocate_textbox(MOL *mp)
674 {
675 mp->ntextboxes++;
676 mp->textboxes = (TEXTBOX_T*) realloc(mp->textboxes, sizeof(TEXTBOX_T) * mp->ntextboxes);
677 mp->textboxes[mp->ntextboxes-1].message = NULL;
678 mp->textboxes[mp->ntextboxes-1].font = NULL;
679 mp->textboxes[mp->ntextboxes-1].pixtext.pixels = NULL;
680 }
681
delete_last_textbox(MOL * mp)682 void delete_last_textbox(MOL *mp)
683 {
684 if (!mp) return;
685 if (mp->ntextboxes)
686 {
687 if (mp->textboxes[mp->ntextboxes-1].message) free(mp->textboxes[mp->ntextboxes-1].message);
688 if (mp->textboxes[mp->ntextboxes-1].font) free(mp->textboxes[mp->ntextboxes-1].font);
689 if (mp->textboxes[mp->ntextboxes-1].pixtext.pixels) free(mp->textboxes[mp->ntextboxes-1].pixtext.pixels);
690 mp->ntextboxes--;
691 mp->textboxes = (TEXTBOX_T*) realloc(mp->textboxes, sizeof(TEXTBOX_T) * mp->ntextboxes);
692 }
693 }
694
delete_ith_textbox(MOL * mp,int it)695 void delete_ith_textbox(MOL *mp, int it)
696 {
697 int i;
698 if (it >= 0 && it < mp->ntextboxes)
699 {
700 if (mp->textboxes[it].message) free(mp->textboxes[it].message);
701 if (mp->textboxes[it].font) free(mp->textboxes[it].font);
702 if (mp->textboxes[it].pixtext.pixels) free(mp->textboxes[it].pixtext.pixels);
703 for(i = it; i < mp->ntextboxes-1; i++)
704 mp->textboxes[i+1] = mp->textboxes[i];
705 mp->ntextboxes--;
706 mp->textboxes = (TEXTBOX_T*) realloc(mp->textboxes, sizeof(TEXTBOX_T) * mp->ntextboxes);
707 }
708 }
709
read_line(FILE * in)710 char *read_line(FILE *in)
711 {
712 #ifdef LINUX
713 size_t nchar = 0;
714 char *line = NULL;
715 ssize_t res;
716
717 if (feof(in))
718 {
719 line = (char*) malloc(sizeof(char));
720 line[0] = 0;
721 return line;
722 }
723 res = getline(&line, &nchar, in);
724 if (res < 0)
725 {
726 free(line);
727 line = (char*) malloc(sizeof(char));
728 line[0] = 0;
729 }
730 return line;
731 #endif
732 #ifdef WINDOWS
733 int i = 1;
734 char c;
735 char *str = NULL;
736 str = (char*) malloc(sizeof(char) * i);
737 str[i-1] = 0;
738 do
739 {
740 c = fgetc(in);
741 if (c == 10 || c == 0 || feof(in)) return str;
742 else
743 {
744 str = (char*) realloc(str, sizeof(char) * ++i);
745 str[i-2] = c;
746 str[i-1] = 0;
747 }
748 }
749 while(c != 10 && c != 0 && !feof(in));
750
751 return str;
752 #endif
753 }
754
get_file_exist(char * filename)755 char get_file_exist(char* filename)
756 {
757 struct stat file_info;
758 int rc;
759 rc=stat(filename, &file_info);
760 if (rc == 0) return 1;
761 else return 0;
762 }
763
get_file_time(char * filename,long * file_ntime)764 time_t get_file_time(char* filename, long *file_ntime)
765 {
766 struct stat file_info;
767 int rc;
768
769 rc=stat(filename, &file_info);
770
771 *file_ntime = file_info.st_mtim.tv_nsec;
772
773 return file_info.st_mtime;
774 }
775
open_file(char * filename,char * file_description)776 void open_file(char *filename, char *file_description)
777 {
778 int i;
779 char *ext = strrchr(filename, '.'); /*POSSIBLE BUG!!!*/
780 char *command;
781 char *newfile;
782
783 #ifdef EBUG
784 printf("opening file |%s| file description |%s|\n", filename, file_description); fflush(stdout);
785 #endif
786
787 delete_all_orbitals_from_list();
788
789 if (input_file_name) free(input_file_name);
790
791 /*if (input_file_type) free(input_file_type);*/
792 input_file_name = strdup(filename);
793
794 if (!input_file_type && file_description) input_file_type = strdup(file_description);
795
796 if (current_directory) free(current_directory);
797
798 current_directory = strdup(filename);
799
800 #ifdef WINDOWS
801 for(i = strlen(current_directory)-1; i > 0 && current_directory[i] != '\\'; i--);
802 #else
803 for(i = strlen(current_directory)-1; i > 0 && current_directory[i] != '/'; i--);
804 #endif
805 if (i != 0) current_directory[i] = 0;
806
807 init_file_time = get_file_time(input_file_name, &init_file_ntime);
808 #ifdef EBUG
809 printf("INITIAL FILE TIME = %ld : %ld\n", init_file_time, init_file_ntime);
810 #endif
811
812 if (ext == NULL)
813 {
814 make_warning("ERROR: Can't open file: unknown filetype!");
815 return;
816 }
817
818 ext += 1;
819
820 /*19.02.2013; now determine right plugin for converting file*/
821 if (file_description) /*determination based on user-selected file type*/
822 {
823 if (strcmp(file_description, LUSCUS_FILE_DES) == 0)
824 {
825 /* if (m)
826 if (m->ngrids)
827 delete_all_orbitals_from_list();*/
828
829 open_gv_file(input_file_name);
830 if (!file_open_success)
831 {
832 init_file_time = 0;
833 close_file();
834 /*make false m structure*/
835 return;
836 }
837 Do_center();
838 luscus_gtk_show_or_hide_widgets();
839 luscus_gtk_menubar_show_or_hide_widgets();
840 luscus_gtk_pop_message_from_statusbar1();
841 luscus_gtk_push_message_to_statusbar1(input_file_name);
842 luscus_gtk_update_3Dobject_info();
843 file_open_success = 0;
844 return;
845 }
846 else
847 {
848 for (i = 0; i < n_input_types; i++)
849 {
850 if (strcmp(input_filetypes[i].description, file_description) == 0) break;
851 }
852 }
853 }
854 else if (strcmp(ext, LUSCUS_FILE_EXT) == 0) /* special option for LUS files !*/
855 {
856 /* if (m)
857 if (m->ngrids)
858 delete_all_orbitals_from_list();*/
859
860 open_gv_file(input_file_name);
861
862 if (!file_open_success)
863 {
864 init_file_time = 0;
865 close_file();
866 /*make false m structure*/
867 return;
868 }
869 Do_center();
870 luscus_gtk_show_or_hide_widgets();
871 luscus_gtk_menubar_show_or_hide_widgets();
872 luscus_gtk_pop_message_from_statusbar1();
873 luscus_gtk_push_message_to_statusbar1(input_file_name);
874 luscus_gtk_update_3Dobject_info();
875 file_open_success = 0;
876 return;
877 }
878 else /*determination based on file extension*/
879 {
880 for (i = 0; i < n_input_types; i++)
881 {
882 if (strcmp(input_filetypes[i].extension, ext) == 0) break;
883 }
884 }
885
886 if (i < n_input_types && input_filetypes[i].libpath)
887 {
888 command = (char*) malloc(sizeof(char) * (strlen(input_filetypes[i].libpath) + strlen(input_filetypes[i].forward) + strlen(filename) + 7));
889 command[0] = 0;
890 /*#ifdef WINDOWS
891 strcat(command, "\"");
892 #endif*/
893 strcat(command, input_filetypes[i].libpath);
894 #ifdef WINDOWS
895 strcat(command, "\\");
896 #else
897 strcat(command, "/");
898 #endif
899 strcat(command, input_filetypes[i].forward);
900 /*#ifdef WINDOWS
901 strcat(command, "\" \"");
902 strcat(command, filename);
903 strcat(command, "\"");
904 #else*/
905 strcat(command, " ");
906 strcat(command, filename);
907 /*#endif*/
908 if (system(command) < 0)
909 {
910 free(command); /*converting file to GV failed*/
911 command = (char*) malloc(sizeof(char) * (strlen(filename) + 27));
912 command[0] = 0;
913 strcat(command,"ERROR: Can't process file: ");
914 strcat(command, filename);
915 make_warning(command);
916 free(command);
917 return;
918 }
919 free(command);
920 }
921 else
922 {
923 /* 08.04.2013: this else diseabled This prevents opening a lus. file!*/
924 /* 03.01.2014: this else enabled again*/
925 make_warning("ERROR: Can't open file: unknown filetype!");
926 return;
927 }
928
929 newfile = make_luscus_file_name(filename);
930
931 /* if (m)
932 if (m->ngrids)
933 delete_all_orbitals_from_list();*/
934
935 open_gv_file(newfile);
936
937 free(newfile);
938
939 if (!file_open_success)
940 {
941 init_file_time = 0;
942 close_file();
943 /*make false m structure*/
944 return;
945 }
946
947 Do_center();
948 luscus_gtk_show_or_hide_widgets();
949 luscus_gtk_menubar_show_or_hide_widgets();
950 luscus_gtk_pop_message_from_statusbar1();
951 luscus_gtk_push_message_to_statusbar1(input_file_name);
952 luscus_gtk_update_3Dobject_info();
953 file_open_success = 0;
954 }
955
open_current_luscus_file(void)956 FILE *open_current_luscus_file(void)
957 {
958 FILE *fp;
959 char *luscus_filename;
960 luscus_filename = make_luscus_file_name(input_file_name);
961 fp = fopen(luscus_filename, "rb");
962 free(luscus_filename);
963 return fp;
964 }
965
open_gv_file(char * filename)966 void open_gv_file(char *filename)
967 {
968 char *msg;
969 FILE *in = NULL;
970
971 in = fopen(filename, "rb");
972 if (in == NULL)
973 {
974 msg = (char*) malloc(sizeof(char) * (strlen(filename) + 25));
975 msg[0] = 0;
976 strcat(msg,"ERROR: Can't open file: ");
977 strcat(msg, filename);
978 make_warning(msg);
979 free(msg);
980 return;
981 }
982
983 remove_backup();
984 parse_gv_file(in);
985 append_backup();
986 fclose(in);
987 if (m->ngrids) do_load_grid();
988 file_open_success = 1;
989 luscus_gtk_update_geo_info();
990 }
991
parse_gv_file(FILE * in)992 void parse_gv_file(FILE *in)
993 {
994 int ngeoms;
995 deallocate_mol_data();
996 ngeoms = read_geo_positions(in);
997 mol = new_mol(ngeoms);
998
999 read_all_sections(in);
1000 draw_all_pixdata();
1001 /* luscus_gtk_update_geo_info();*/
1002 }
1003
read_geo_positions(FILE * in)1004 int read_geo_positions(FILE *in)
1005 {
1006 int i, j;
1007 long lastpos;
1008 char *line;
1009 int ngeo;
1010 int sr45 = 0;
1011 ngeo = 1; /*file beginning one geometry is present*/
1012 geo_file_position = malloc(sizeof(long));
1013 geo_file_position[0] = ftell(in);
1014
1015 line = read_line(in);
1016 while(!feof(in) && strcasestr(line, "<break>") == 0)
1017 {
1018 free(line);
1019 line = read_line(in);
1020 if (feof(in))
1021 {
1022 if (sr45) ngeo--;
1023 break;
1024 }
1025
1026 if (strcasestr(line, "<grid>"))
1027 {
1028 char *ptr;
1029 long orboff;
1030
1031 for(i = 0; i < 8; i++)
1032 {
1033 free(line);
1034 line = read_line(in);
1035 }
1036 ptr = strchr(line, '=') + 1;
1037 orboff = atol(ptr);
1038
1039 fseek(in, orboff, SEEK_CUR);
1040 }
1041 else if (strcasestr(line, "<break>"))
1042 {
1043 if (sr45) ngeo--;
1044 break;
1045 }
1046 else if (strcasestr(line, "<end>"))
1047 {
1048 ngeo++;
1049 sr45 = 1;
1050 geo_file_position = realloc(geo_file_position, sizeof(long) * ngeo);
1051 geo_file_position[ngeo-1] = ftell(in);
1052 }
1053 else sr45 = 0;
1054 }
1055 if (line) free(line);
1056
1057 /*This detects the "<end>" at the end of the file*/
1058 lastpos = ftell(in);
1059 geo_file_position = realloc(geo_file_position, sizeof(long) * ngeo);
1060
1061 /* rewind(in);*/ /*Do I need this?*/
1062 return ngeo;
1063 }
1064
read_all_sections(FILE * in)1065 void read_all_sections(FILE *in)
1066 {
1067 int i;
1068 for (i = 0; i < n_geometries; i++)
1069 {
1070 should_bonding_be_determined = 1;
1071 read_section(in, i, mol+i);
1072 }
1073 insert_all_atoms_into_list();
1074 }
1075
read_all_grids_from_all_sections(int * m_ngrids)1076 void read_all_grids_from_all_sections(int *m_ngrids)
1077 {
1078 int m_igrid = 0;
1079 FILE *in;
1080 char *line;
1081
1082 if (!input_file_name)
1083 {
1084 luscus_gtk_pop_message_from_statusbar1();
1085 /* luscus_gtk_push_message_to_statusbar1("ERROR: Can't find file that contains grid data!");*/
1086 return;
1087 }
1088 in = open_current_luscus_file();
1089 if (in == NULL)
1090 {
1091 luscus_gtk_push_message_to_statusbar1("ERROR: Can't open file that contains grid data!");
1092 return;
1093 }
1094
1095 while(!feof(in))
1096 {
1097 line = read_line(in);
1098 if (strcasestr(line, "<end>")) m_igrid++;
1099 if (strcasestr(line, "<grid>"))
1100 if (m_ngrids[m_igrid]) read_grid_block(in, mol + m_igrid);
1101 free(line);
1102 }
1103
1104 fclose(in);
1105 }
1106
read_section(FILE * in,int isec,MOL * mp)1107 void read_section(FILE *in, int isec, MOL *mp)
1108 {
1109 int next = 1;
1110 char *line;
1111
1112 if (isec >= n_geometries) return;
1113 fseek(in, geo_file_position[isec], 0);
1114 read_coord_block(in, mp);
1115 while(next)
1116 {
1117 if (feof(in)) next = 0;
1118 if (next)
1119 {
1120 line = read_line(in);
1121 if (strcasestr(line, "<atom>")) read_additional_atom_block(in, mp);
1122 else if (strcasestr(line, "<energy>")) read_energy_block(in, mp);
1123 else if (strcasestr(line, "<rms_grad>")) read_rms_grad_block(in, mp);
1124 else if (strcasestr(line, "<max_grad>")) read_max_grad_block(in, mp);
1125 else if (strcasestr(line, "<bond>")) read_bond_block(in, mp);
1126 else if (strcasestr(line, "<dipole>")) read_dipole_block(in, mp);
1127 else if (strcasestr(line, "<vector>")) read_vector_block(in, mp);
1128 else if (strcasestr(line, "<triangle>")) read_triangle_block(in, mp);
1129 else if (strcasestr(line, "<sphere>")) read_sphere_block(in, mp);
1130 else if (strcasestr(line, "<surface>")) read_surface_block(in, mp);
1131 else if (strcasestr(line, "<cell>")) read_cell_block(in, mp);
1132 else if (strcasestr(line, "<textbox>")) read_textbox_block(in, mp);
1133 else if (strcasestr(line, "<grid>")) read_grid_block(in, mp);
1134 else if (strcasestr(line, "<vibration>")) read_vibration_block(in, mp);
1135 else if (strcasestr(line, "<editable>")) read_editable_block(in, mp); /*19.02.2013 new block!*/
1136 else if (strcasestr(line, "<sleep>")) read_sleep_block(in);
1137 else if (strcasestr(line, "<write>")) read_write_block(in);
1138 else if (strcasestr(line, "<end>")) next = 0;
1139 free(line);
1140 }
1141 }
1142 if (should_bonding_be_determined) determine_bonding(mp);
1143
1144 set_sizes();
1145 }
1146
read_coord_block(FILE * in,MOL * mp)1147 void read_coord_block(FILE *in, MOL *mp)
1148 {
1149 int i;
1150 char *line;
1151 int natom;
1152 int nread = 0;
1153 char *tmp;
1154
1155 line = read_line(in);
1156 if (line)
1157 {
1158 natom = atoi(line);
1159 allocate_atoms(mp, natom);
1160 }
1161 else
1162 {
1163 deallocate_mol_data();
1164 make_warning("ERROR: Can't read file!");
1165 return;
1166 }
1167 free(line);
1168 line = read_line(in); free(line); /*comment line*/
1169 for(i = 0; i < natom; i++)
1170 {
1171 line = read_line(in);
1172 if (line)
1173 {
1174 if (line[0])
1175 {
1176 mp->elem[i].name = get_one_word(get_ptr_ith_word(line,1));
1177 set_element_data(mp, i);
1178
1179 tmp = get_ptr_ith_word(line, 2);
1180 if (tmp) mp->xyz[i][0] = atof(tmp);
1181 else mp->xyz[i][0] = 0.0;
1182
1183 tmp = get_ptr_ith_word(line, 3);
1184 if (tmp) mp->xyz[i][1] = atof(tmp);
1185 else mp->xyz[i][1] = 0.0;
1186
1187 tmp = get_ptr_ith_word(line, 4);
1188 if (tmp) mp->xyz[i][2] = atof(tmp);
1189 else mp->xyz[i][2] = 0.0;
1190 nread++;
1191 }
1192 free(line);
1193 }
1194 }
1195 if (nread != natom)
1196 {
1197 make_warning("WARNING: Broken file; some atoms are missing!\n");
1198 /*In the case of broken file!!!*/
1199 allocate_atoms(mp, nread);
1200 }
1201
1202 }
1203
read_additional_atom_block(FILE * in,MOL * mp)1204 void read_additional_atom_block(FILE *in, MOL *mp)
1205 {
1206 char *line;
1207 int next = 1;
1208 int iatom;
1209
1210 if (mp->natom <= 0) return;
1211
1212 iatom = 0;
1213
1214 while(next)
1215 {
1216 line = read_line(in);
1217 if (feof(in)) next = 0;
1218 if (iatom == mp->natom) next = 0; /*Check that file doesn't contain exess of atoms*/
1219 if (next)
1220 {
1221 if (strcasestr(line, "</atom>")) next = 0;
1222 if (strcasestr(line,"name"))
1223 {
1224 mp->name[iatom] = my_read_str_value((char*) strcasestr(line,"name"));
1225 if (!(mp->ishow & HAS_ATOMNAMES)) mp->ishow ^= HAS_ATOMNAMES;
1226 }
1227 if ((char*) strcasestr(line,"mulliken_charge"))
1228 {
1229 mp->charge_m[iatom] = my_read_dble_value((char*) strcasestr(line,"mulliken_charge"));
1230 if (!(mp->ishow & HAS_MULLIKEN)) mp->ishow ^= HAS_MULLIKEN;
1231 }
1232 if (strcasestr(line,"loprop_charge"))
1233 {
1234 mp->charge_l[iatom] = my_read_dble_value((char*) strcasestr(line,"loprop_charge"));
1235 if (!(mp->ishow & HAS_LOPROP)) mp->ishow ^= HAS_LOPROP;
1236 }
1237 if (strcasestr(line,"number"))
1238 {
1239 mp->additional_numeration[iatom] = my_read_int_value((char*) strcasestr(line,"number"));
1240 if (!(mp->ishow & HAS_ATOMNUMS)) mp->ishow ^= HAS_ATOMNUMS;
1241 }
1242 if (strcasestr(line,"symmetry"))
1243 {
1244 mp->symmetry[iatom] = my_read_int_value((char*) strcasestr(line,"symmetry"));
1245 if (!(mp->ishow & HAS_ATOMSYMM)) mp->ishow ^= HAS_ATOMSYMM;
1246 }
1247 /*added 08.04.2013 atom colors*/
1248 if (strcasestr(line,"red"))
1249 {
1250 mp->elem[iatom].color[0] = my_read_dble_value((char*) strcasestr(line,"red"));
1251 if (!(mp->ishow & HAS_COLOUR)) mp->ishow ^= HAS_COLOUR;
1252 }
1253 if (strcasestr(line,"green"))
1254 {
1255 mp->elem[iatom].color[1] = my_read_dble_value((char*) strcasestr(line,"green"));
1256 if (!(mp->ishow & HAS_COLOUR)) mp->ishow ^= HAS_COLOUR;
1257 }
1258 if (strcasestr(line,"blue"))
1259 {
1260 mp->elem[iatom].color[2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1261 if (!(mp->ishow & HAS_COLOUR)) mp->ishow ^= HAS_COLOUR;
1262 }
1263 iatom++;
1264 }
1265 free(line);
1266 }
1267
1268 }
1269
read_energy_block(FILE * in,MOL * mp)1270 void read_energy_block(FILE *in, MOL *mp)
1271 {
1272 char *line;
1273 line = read_line(in);
1274 mp->geo_energy = atof(line);
1275 if (!(mp->ishow & HAS_ENERGY)) mp->ishow ^= HAS_ENERGY;
1276 free(line);
1277 }
1278
read_rms_grad_block(FILE * in,MOL * mp)1279 void read_rms_grad_block(FILE *in, MOL *mp)
1280 {
1281 char *line;
1282 line = read_line(in);
1283 mp->max_grad = atof(line);
1284 if (!(mp->ishow & HAS_RMS_GRAD)) mp->ishow ^= HAS_RMS_GRAD;
1285 free(line);
1286 }
1287
read_max_grad_block(FILE * in,MOL * mp)1288 void read_max_grad_block(FILE *in, MOL *mp)
1289 {
1290 char *line;
1291 line = read_line(in);
1292 mp->rms_grad = atof(line);
1293 if (!(mp->ishow & HAS_MAX_GRAD)) mp->ishow ^= HAS_MAX_GRAD;
1294 free(line);
1295 }
1296
read_bond_block(FILE * in,MOL * mp)1297 void read_bond_block(FILE *in, MOL *mp)
1298 {
1299 char *line;
1300 int next = 1;
1301 int at1, at2, bo;
1302
1303 line = read_line(in);
1304 if (!line_is_empty(line))
1305 if (my_read_int_value((char*) strcasestr(line,"automatic")) == 0) should_bonding_be_determined = 0;
1306 free(line);
1307
1308 while(next)
1309 {
1310 line = read_line(in);
1311 if (feof(in)) next = 0;
1312
1313 if (next)
1314 {
1315 if (strcasestr(line, "</bond>")) next = 0;
1316 else
1317 {
1318 sscanf(line, "%d %d %d", &at1, &at2, &bo);
1319 if (at1 > 0 && at1 <= mp->natom)
1320 if (at2 > 0 && at2 <= mp->natom)
1321 add_bond(mp, at1 - 1, at2 - 1, bo);
1322 }
1323 }
1324 free(line);
1325 }
1326
1327 }
1328
determine_bonding(MOL * mp)1329 void determine_bonding(MOL *mp)
1330 {
1331 int i, j;
1332 int nbond = 0;
1333 int nne1, nne2;
1334 #ifdef EBUG
1335 printf("entering function: determine bonding!\n"); fflush(stdout);
1336 #endif
1337 if (mp->natom > MAX_ATOMS_FOR_BONDING) return;
1338 /*connect atoms with single bonds*/
1339 for(i = 0; i < mp->natom - 1; i++)
1340 for(j = i + 1; j < mp->natom; j++)
1341 {
1342 #ifndef NO_RESTRICT_BOND
1343 if (fabs(mp->xyz[i][0] - mp->xyz[j][0]) < 3.0 &&
1344 fabs(mp->xyz[i][1] - mp->xyz[j][1]) < 3.0 &&
1345 fabs(mp->xyz[i][2] - mp->xyz[j][2]) < 3.0)
1346 {
1347 #endif
1348 /* should be filter that does not check atoms that are some distance apart!*/
1349 if (find_bond(mp, i, j) == -1)
1350 {
1351 if (atom_distance(mp, i, j) < 1.15 *(mp->elem[i].bond_rad + mp->elem[j].bond_rad))
1352 {
1353 nne1 = get_number_of_bonds_on_atom(i, mp);
1354 nne2 = get_number_of_bonds_on_atom(j, mp);
1355 if (nne1 < mp->elem[i].valency && nne2 < mp->elem[j].valency)
1356 {
1357 nbond++;
1358 add_bond(mp, i, j, SINGLE_BOND);
1359 }
1360 }
1361 }
1362 #ifndef NO_RESTRICT_BOND
1363 }
1364 #endif
1365 }
1366
1367 #ifdef EBUG
1368 printf("bond order one = %d\n", Input_Data.bond_order_one);
1369 #endif
1370 if (Input_Data.bond_order_one) return;
1371
1372 /*find unsatisfied atoms and increase bond order*/
1373 for(i = 0; i < mp->nbond; i++)
1374 if (mp->bond[i].bond_type == SINGLE_BOND)
1375 {
1376 nne1 = get_number_of_bonds_on_atom(mp->bond[i].iat1, mp);
1377 nne2 = get_number_of_bonds_on_atom(mp->bond[i].iat2, mp);
1378 if (nne1 < m->elem[mp->bond[i].iat1].valency && nne2 < m->elem[mp->bond[i].iat2].valency)
1379 mp->bond[i].bond_type = DOUBLE_BOND;
1380 }
1381
1382 /*search for triple bonds*/
1383 for(i = 0; i < mp->nbond; i++)
1384 if (mp->bond[i].bond_type == DOUBLE_BOND)
1385 {
1386 nne1 = get_number_of_bonds_on_atom(mp->bond[i].iat1, mp);
1387 nne2 = get_number_of_bonds_on_atom(mp->bond[i].iat2, mp);
1388 if (nne1 < m->elem[mp->bond[i].iat1].valency && nne2 < m->elem[mp->bond[i].iat2].valency)
1389 mp->bond[i].bond_type = TRIPLE_BOND;
1390 }
1391
1392 /*search for conjugations*/
1393 /* for(i = 0; i < mp->nbond; i++)
1394 if (mp->bond[i].bond_type == SINGLE_BOND)
1395 {
1396 nne1 = get_number_of_bond_types_on_atom(mp, mp->bond[i].iat1, DOUBLE_BOND);
1397 nne2 = get_number_of_bond_types_on_atom(mp, mp->bond[i].iat2, DOUBLE_BOND);
1398 if (nne1 == 1 && nne2 == 1)
1399 {
1400 mp->bond[i].bond_type = S_P_BOND;
1401 for(j = 0; j < mp->nbond; j++)
1402 if (mp->bond[j].bond_type == DOUBLE_BOND)
1403 {
1404 if (mp->bond[j].iat1 == mp->bond[i].iat1 || mp->bond[j].iat2 == mp->bond[i].iat1)
1405 mp->bond[j].bond_type = S_P_BOND;
1406 if (mp->bond[j].iat1 == mp->bond[i].iat2 || mp->bond[j].iat2 == mp->bond[i].iat2)
1407 mp->bond[j].bond_type = S_P_BOND;
1408 }
1409 }
1410 }*/
1411
1412 #ifdef EBUG
1413 printf("leaving function: determine bonding!\n"); fflush(stdout);
1414 #endif
1415 }
1416
get_number_of_bonds_on_atom(int iatom,MOL * mp)1417 int get_number_of_bonds_on_atom(int iatom, MOL *mp)
1418 {
1419 int i;
1420 int ibond = 0;
1421 for(i = 0; i < mp->nbond; i++)
1422 if (mp->bond[i].iat1 == iatom || mp->bond[i].iat2 == iatom)
1423 if (mp->bond[i].bond_type == SINGLE_BOND) ibond++;
1424 else if (mp->bond[i].bond_type == DOUBLE_BOND) ibond+=2;
1425 else if (mp->bond[i].bond_type == TRIPLE_BOND) ibond+=3;
1426
1427 return ibond;
1428 }
1429
1430 /*void adapt_bonding(MOL *mp)
1431 {
1432 float *bondor;
1433 float *atomval;
1434 int i, j, k;
1435
1436 bondor = (float*) malloc(mp->nbond * sizeof(float));
1437 atomval = (float*) malloc(mp->natom * sizeof(float));
1438
1439 /-*initialize values*-/
1440 for (i = 0; i < mp->nbond; i++)
1441 {
1442 switch(mp->bond[i].bond_type)
1443 {
1444 case NO_BOND: bondor[i] = 0.F;
1445 case SINGLE_BOND: bondor[i] = 1.F;
1446 case DOUBLE_BOND: bondor[i] = 2.F;
1447 case TRIPLE_BOND: bondor[i] = 3.F;
1448 case PARTIAL_BOND: bondor[i] = 0.5F;
1449 case S_P_BOND: bondor[i] = 1.5F;
1450 case LINE_BOND: bondor[i] = 0.F;
1451 }
1452 printf("bondors = %f\n", bondor[i]);
1453 }
1454
1455 for(i = 0; i < mp->natom; i++) /-*number of iterations = number of atoms; could be changed*-/
1456 {
1457 /-*1. calculate current valencies*-/
1458 for(j = 0; j < mp->natom; j++)
1459 for(k = 0; k < mp->nbond; k++)
1460 if (j == mp->bond[k].iat1 || j == mp->bond[k].iat2)
1461 atomval[j] += (float) bondor[k];
1462
1463 /-*2. correct bond orders*-/
1464 for(j = 0; j < mp->natom; j++)
1465 for(k = 0; k < mp->nbond; k++)
1466 {
1467 if (j == mp->bond[k].iat1 || j == mp->bond[k].iat2)
1468 bondor[k] += ((float) mp->elem[j].valency - atomval[j]) / (float) (2* mp->elem[j].valency);
1469 printf("bondor[%d] = %f\n", k, bondor[k]);
1470 }
1471 putchar(10);
1472 }
1473
1474 /-*transform corrected bond orders into the bond types*-/
1475 for(i = 0; i < mp->nbond; i++)
1476 {
1477 if (bondor[i] < 0.25F) mp->bond[i].bond_type = NO_BOND;
1478 else if (bondor[i] < 0.5F) mp->bond[i].bond_type = PARTIAL_BOND;
1479 else if (bondor[i] < 1.25F) mp->bond[i].bond_type = SINGLE_BOND;
1480 else if (bondor[i] < 1.75F) mp->bond[i].bond_type = S_P_BOND;
1481 else if (bondor[i] < 2.5F) mp->bond[i].bond_type = DOUBLE_BOND;
1482 else mp->bond[i].bond_type = TRIPLE_BOND;
1483 }
1484
1485 free(bondor);
1486 free(atomval);
1487 }*/
1488
1489 /*int get_number_of_bond_types_on_atom(MOL *mp, int iatom, int bond_type)
1490 {
1491 int i;
1492 int nt = 0;
1493 for(i = 0; i < mp->nbond; i++)
1494 if (mp->bond[i].iat1 == iatom || mp->bond[i].iat2 == iatom)
1495 if (mp->bond[i].bond_type == bond_type) nt++;
1496 return nt;
1497 }*/
1498
read_dipole_block(FILE * in,MOL * mp)1499 void read_dipole_block(FILE *in, MOL *mp)
1500 {
1501 char *line;
1502
1503 line = read_line(in);
1504 sscanf(line, "%lf %lf %lf", &mp->dipole[0], &mp->dipole[1], &mp->dipole[2]);
1505
1506 if (!(mp->ishow & HAS_DIPOLE)) mp->ishow ^= HAS_DIPOLE;
1507
1508 free(line);
1509 }
1510
read_vector_block(FILE * in,MOL * mp)1511 void read_vector_block(FILE *in, MOL *mp)
1512 {
1513 char *line;
1514 int nvec = mp->nvector + 1;
1515
1516 allocate_vectors(mp, nvec);
1517 vector_set_default_values(mp, nvec-1);
1518 line = read_line(in);
1519
1520 if (strcasestr(line,"red"))
1521 mp->vector_color[nvec-1][0] = my_read_dble_value((char*) strcasestr(line,"red"));
1522 if (strcasestr(line,"green"))
1523 mp->vector_color[nvec-1][1] = my_read_dble_value((char*) strcasestr(line,"green"));
1524 if (strcasestr(line,"blue"))
1525 mp->vector_color[nvec-1][2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1526 if (strcasestr(line,"transparency"))
1527 mp->vector_color[nvec-1][3] = my_read_dble_value((char*) strcasestr(line,"transparency"));
1528 if (strcasestr(line,"radius"))
1529 mp->radius[nvec-1] = my_read_dble_value((char*) strcasestr(line,"radius"));
1530 if (strcasestr(line,"sharpness"))
1531 mp->sharpness[nvec-1] = my_read_dble_value((char*) strcasestr(line,"sharpness"));
1532
1533 free(line);
1534
1535 line = read_line(in);
1536 sscanf(line, "%lf %lf %lf", &mp->vector1[nvec-1][0], &mp->vector1[nvec-1][1], &mp->vector1[nvec-1][2]);
1537 free(line);
1538
1539 line = read_line(in);
1540 sscanf(line, "%lf %lf %lf", &mp->vector2[nvec-1][0], &mp->vector2[nvec-1][1], &mp->vector2[nvec-1][2]);
1541 free(line);
1542 }
1543
vector_set_default_values(MOL * mp,int ivec)1544 void vector_set_default_values(MOL *mp, int ivec)
1545 {
1546 int i;
1547 if (ivec >= mp->nvector) return;
1548 for (i = 0; i < 4; i++)
1549 mp->vector_color[ivec][i] = Input_Data.extracolor[i];
1550 /*default vector radius and sharpness here*/
1551 mp->radius[ivec] = 0.05F;
1552 mp->sharpness[ivec] = 0.1F;
1553 }
1554
read_triangle_block(FILE * in,MOL * mp)1555 void read_triangle_block(FILE *in, MOL *mp)
1556 {
1557 char *line;
1558 int ntriang = mp->ntriangle + 1;
1559
1560 allocate_triangles(mp, ntriang);
1561 triangle_set_defeault_values(mp, ntriang-1);
1562 line = read_line(in);
1563 if (strcasestr(line,"red"))
1564 mp->triangle_color[ntriang-1][0] = my_read_dble_value((char*) strcasestr(line,"red"));
1565 if (strcasestr(line,"green"))
1566 mp->triangle_color[ntriang-1][1] = my_read_dble_value((char*) strcasestr(line,"green"));
1567 if (strcasestr(line,"blue"))
1568 mp->triangle_color[ntriang-1][2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1569 if (strcasestr(line,"transparency"))
1570 mp->triangle_color[ntriang-1][3] = my_read_dble_value((char*) strcasestr(line,"transparency"));
1571 free(line);
1572
1573 line = read_line(in);
1574 sscanf(line, "%lf %lf %lf", &mp->triangle1[ntriang-1][0], &mp->triangle1[ntriang-1][1], &mp->triangle1[ntriang-1][2]);
1575 free(line);
1576
1577 line = read_line(in);
1578 sscanf(line, "%lf %lf %lf", &mp->triangle2[ntriang-1][0], &mp->triangle2[ntriang-1][1], &mp->triangle2[ntriang-1][2]);
1579 free(line);
1580
1581 line = read_line(in);
1582 sscanf(line, "%lf %lf %lf", &mp->triangle3[ntriang-1][0], &mp->triangle3[ntriang-1][1], &mp->triangle3[ntriang-1][2]);
1583 free(line);
1584 }
1585
triangle_set_defeault_values(MOL * mp,int itriang)1586 void triangle_set_defeault_values(MOL *mp, int itriang)
1587 {
1588 int i;
1589 if (itriang >= mp->ntriangle) return;
1590 for(i = 0; i < 4; i++)
1591 mp->triangle_color[itriang][i] = Input_Data.extracolor[i];
1592 }
1593
read_sphere_block(FILE * in,MOL * mp)1594 void read_sphere_block(FILE *in, MOL *mp)
1595 {
1596 char *line;
1597 int nspheres = mp->nsphere + 1;
1598
1599 allocate_spheres(mp, nspheres);
1600 sphere_set_defeault_values(mp, nspheres-1);
1601
1602 line = read_line(in);
1603 if (strcasestr(line,"red"))
1604 mp->sphere_color[nspheres-1][0] = my_read_dble_value((char*) strcasestr(line,"red"));
1605 if (strcasestr(line,"green"))
1606 mp->sphere_color[nspheres-1][1] = my_read_dble_value((char*) strcasestr(line,"green"));
1607 if (strcasestr(line,"blue"))
1608 mp->sphere_color[nspheres-1][2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1609 if (strcasestr(line,"transparency"))
1610 mp->sphere_color[nspheres-1][3] = my_read_dble_value((char*) strcasestr(line,"transparency"));
1611 if (strcasestr(line,"radius"))
1612 mp->sphere_radius[nspheres-1] = my_read_dble_value((char*) strcasestr(line,"radius"));
1613
1614 free(line);
1615
1616 line = read_line(in);
1617 sscanf(line, "%lf %lf %lf", &mp->sphere_center[nspheres-1][0], &mp->sphere_center[nspheres-1][1], &mp->sphere_center[nspheres-1][2]);
1618 free(line);
1619 }
1620
sphere_set_defeault_values(MOL * mp,int isphere)1621 void sphere_set_defeault_values(MOL *mp, int isphere)
1622 {
1623 int i;
1624 if (isphere >= mp->nsphere) return;
1625 for(i = 0; i < 4; i++)
1626 mp->sphere_color[isphere][i] = Input_Data.extracolor[i];
1627 mp->sphere_radius[isphere] = 2.0;
1628 }
1629
read_surface_block(FILE * in,MOL * mp)1630 void read_surface_block(FILE *in, MOL *mp)
1631 {
1632 char *line;
1633 int nsurf = mp->nsurf+1;
1634
1635 allocate_surfaces(mp, nsurf);
1636 surface_set_default_values(mp, nsurf-1);
1637 line = read_line(in);
1638
1639 if (strcasestr(line,"red"))
1640 mp->surf_color[nsurf-1][0] = my_read_dble_value((char*) strcasestr(line,"red"));
1641 if (strcasestr(line,"green"))
1642 mp->surf_color[nsurf-1][1] = my_read_dble_value((char*) strcasestr(line,"green"));
1643 if (strcasestr(line,"blue"))
1644 mp->surf_color[nsurf-1][2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1645 if (strcasestr(line,"transparency"))
1646 mp->surf_color[nsurf-1][3] = my_read_dble_value((char*) strcasestr(line,"transparency"));
1647
1648
1649 free(line);
1650
1651 line = read_line(in);
1652 sscanf(line, "%lf %lf %lf", &mp->surf1[nsurf-1][0], &mp->surf1[nsurf-1][1], &mp->surf1[nsurf-1][2]);
1653 free(line);
1654
1655 line = read_line(in);
1656 sscanf(line, "%lf %lf %lf", &mp->surf2[nsurf-1][0], &mp->surf2[nsurf-1][1], &mp->surf2[nsurf-1][2]);
1657 free(line);
1658
1659 line = read_line(in);
1660 sscanf(line, "%lf %lf %lf", &mp->surf3[nsurf-1][0], &mp->surf3[nsurf-1][1], &mp->surf3[nsurf-1][2]);
1661 free(line);
1662 }
1663
surface_set_default_values(MOL * mp,int isurf)1664 void surface_set_default_values(MOL *mp, int isurf)
1665 {
1666 int i;
1667 if (isurf >= mp->nsurf) return;
1668 for(i = 0; i < 4; i++)
1669 mp->surf_color[isurf][i] = Input_Data.extracolor[i];
1670 }
1671
read_cell_block(FILE * in,MOL * mp)1672 void read_cell_block(FILE *in, MOL *mp)
1673 {
1674 char *line;
1675 int next = 1;
1676 int i;
1677 int ncell;
1678
1679 ncell = mp->ncells + 1;
1680 allocate_cells(mp, ncell);
1681 cell_set_default_values(mp, ncell-1);
1682
1683 line = read_line(in);
1684
1685 if (strcasestr(line,"red"))
1686 mp->cell_color[ncell-1][0] = my_read_dble_value((char*) strcasestr(line,"red"));
1687 if (strcasestr(line,"green"))
1688 mp->cell_color[ncell-1][1] = my_read_dble_value((char*) strcasestr(line,"green"));
1689 if (strcasestr(line,"blue"))
1690 mp->cell_color[ncell-1][2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1691 if (strcasestr(line,"transparency"))
1692 mp->cell_color[ncell-1][3] = my_read_dble_value((char*) strcasestr(line,"transparency"));
1693
1694 free(line);
1695
1696 line = read_line(in);
1697 sscanf(line, "%lf %lf %lf", &mp->cell1[ncell-1][0], &mp->cell1[ncell-1][1], &mp->cell1[ncell-1][2]);
1698 free(line);
1699
1700 line = read_line(in);
1701 sscanf(line, "%lf %lf %lf", &mp->cell2[ncell-1][0], &mp->cell2[ncell-1][1], &mp->cell2[ncell-1][2]);
1702 free(line);
1703
1704 line = read_line(in);
1705 sscanf(line, "%lf %lf %lf", &mp->cell3[ncell-1][0], &mp->cell3[ncell-1][1], &mp->cell3[ncell-1][2]);
1706 free(line);
1707
1708 line = read_line(in);
1709 sscanf(line, "%lf %lf %lf", &mp->cell4[ncell-1][0], &mp->cell4[ncell-1][1], &mp->cell4[ncell-1][2]);
1710 free(line);
1711 }
1712
cell_set_default_values(MOL * mp,int icell)1713 void cell_set_default_values(MOL *mp, int icell)
1714 {
1715 int i;
1716 if (icell >= mp->ncells) return;
1717 for(i = 0; i < 4; i++)
1718 mp->cell_color[icell][i] = Input_Data.extracolor[i];
1719 }
1720
read_textbox_block(FILE * in,MOL * mp)1721 void read_textbox_block(FILE *in, MOL *mp)
1722 {
1723 char *line;
1724 int itb;
1725
1726 allocate_textbox(mp);
1727 itb = mp->ntextboxes-1;
1728 textbox_set_dfault_values(mp, itb);
1729
1730 line = read_line(in);
1731
1732 if (strcasestr(line,"red"))
1733 mp->textboxes[itb].color[0] = my_read_dble_value((char*) strcasestr(line,"red"));
1734 if (strcasestr(line,"green"))
1735 mp->textboxes[itb].color[1] = my_read_dble_value((char*) strcasestr(line,"green"));
1736 if (strcasestr(line,"blue"))
1737 mp->textboxes[itb].color[2] = my_read_dble_value((char*) strcasestr(line,"blue"));
1738 if (strcasestr(line,"x"))
1739 mp->textboxes[itb].coord_x = my_read_int_value((char*) strcasestr(line,"x"));
1740 if (strcasestr(line,"y"))
1741 mp->textboxes[itb].coord_y = my_read_int_value((char*) strcasestr(line,"y"));
1742
1743 free(line);
1744
1745 line = read_line(in);
1746 if (line[strlen(line)-1] == 10) line[strlen(line)-1] = 0;
1747 mp->textboxes[itb].font = strdup(line);
1748 free(line);
1749
1750 line = read_line(in);
1751 if (line[strlen(line)-1] == 10) line[strlen(line)-1] = 0;
1752 mp->textboxes[itb].message = strdup(line);
1753 free(line);
1754 }
1755
textbox_set_dfault_values(MOL * mp,int itb)1756 void textbox_set_dfault_values(MOL* mp, int itb)
1757 {
1758 int i;
1759
1760 if (itb >= mp->ntextboxes) return;
1761 for(i = 0; i < 3; i++)
1762 mp->textboxes[itb].color[i] = Input_Data.label_color[i];
1763 mp->textboxes[itb].coord_x = 1;
1764 mp->textboxes[itb].coord_x = 2;
1765 mp->textboxes[itb].font = NULL;
1766 mp->textboxes[itb].message = NULL;
1767 }
1768
read_grid_block(FILE * in,MOL * mp)1769 void read_grid_block(FILE *in, MOL *mp)
1770 {
1771 int i, j, ii;
1772 char *line;
1773 char *word;
1774 char *ptr;
1775 int next = 1;
1776 int i0, i1, i2;
1777 int ie0, ie1, ie2;
1778 int icount;
1779 int ishow, iaia;
1780 double pp[3];
1781 double myR529 = 1.0;
1782
1783 long int rec0, rec1;
1784 char buf[512];
1785
1786
1787
1788 line = read_line(in);
1789 /*failsafe: if the line is empty; return*/
1790 if (strlen(line) < 2) return;
1791
1792 deallocate_grids(mp);
1793
1794 #ifdef EBUG
1795 printf("Start reading grid section\n");
1796 #endif
1797
1798 if (strcasestr(line,"N_OF_MO"))
1799 mp->grid.nmo = my_read_int_value((char*) strcasestr(line,"N_OF_MO"));
1800 if (strcasestr(line,"N_of_Grids"))
1801 mp->ngrids = mp->grid.ngrid = my_read_int_value((char*) strcasestr(line,"N_of_Grids"));
1802 if (strcasestr(line,"N_of_Points"))
1803 mp->grid.npoints = my_read_int_value((char*) strcasestr(line,"N_of_Points"));
1804 if (strcasestr(line,"Block_Size"))
1805 mp->grid.block_size = my_read_int_value((char*) strcasestr(line,"Block_Size"));
1806 if (strcasestr(line,"N_Blocks"))
1807 mp->grid.nBlock = my_read_int_value((char*) strcasestr(line,"N_Blocks"));
1808 if (strcasestr(line,"Is_cutoff"))
1809 mp->grid.isCutOff = my_read_int_value((char*) strcasestr(line,"Is_cutoff"));
1810 if (strcasestr(line,"CutOff"))
1811 mp->grid.CutOff = my_read_dble_value((char*) strcasestr(line,"CutOff"));
1812 if (strcasestr(line,"N_P"))
1813 mp->grid.nPointsCutOff = my_read_int_value((char*) strcasestr(line,"N_P"));
1814 if (strcasestr(line,"BOHR")) myR529 = R529;
1815 free(line);
1816
1817 #ifdef EBUG
1818 printf("N_OF_MO = %d\n", mp->grid.nmo);
1819 printf("N_of_Grids = %d\n", mp->ngrids);
1820 printf("N_of_Points = %d\n", mp->grid.npoints);
1821 printf("Block_Size = %d\n", mp->grid.block_size);
1822 printf("N_Blocks = %d\n", mp->grid.nBlock);
1823 printf("Is_cutoff = %d\n", mp->grid.isCutOff);
1824 printf("CutOff = %f\n", mp->grid.CutOff);
1825 printf("N_P = %d\n", mp->grid.nPointsCutOff); fflush(stdout);
1826 printf("Read header line\n");
1827 #endif
1828
1829 allocate_grids(mp, mp->grid.ngrid);
1830
1831 allocate_msrf(mp->grid.ngrid);
1832 line = read_line(in);
1833 ptr = strchr(line, '=') + 1;
1834 sscanf(ptr, "%d %d %d %d %d %d %d", &(mp->grid.iniIndex[0]),
1835 &(mp->grid.iniIndex[1]), &(mp->grid.iniIndex[2]),
1836 &(mp->grid.iniIndex[3]), &(mp->grid.iniIndex[4]),
1837 &(mp->grid.iniIndex[5]), &(mp->grid.iniIndex[6]));
1838 free(line);
1839
1840 line = read_line(in);
1841 ptr = strchr(line, '=') + 1;
1842 sscanf (ptr, "%d %d %d", &(mp->grid.npt[0]), &(mp->grid.npt[1]), &(mp->grid.npt[2]));
1843
1844 free(line);
1845
1846 if (mp->grid.npoints != mp->grid.npt[0] * mp->grid.npt[1] * mp->grid.npt[2])
1847 {
1848 deallocate_grids(mp);
1849 make_warning("ERROR: Ill-defined number of grids");
1850 /*read lined untill </GRID>*/
1851 return;
1852 }
1853 if (mp->grid.ngrid == 0)
1854 {
1855 deallocate_grids(mp);
1856 make_warning("ERROR: No data in grid section");
1857 return;
1858 }
1859
1860 line = read_line(in);
1861 ptr = strchr(line, '=') + 1;
1862 sscanf (ptr, "%lf %lf %lf", &(mp->grid.origin[0]), &(mp->grid.origin[1]), &(mp->grid.origin[2]));
1863 free(line);
1864
1865 line = read_line(in);
1866 ptr = strchr(line, '=') + 1;
1867 sscanf (ptr, "%lf %lf %lf", &(mp->grid.axisvec1[0]),
1868 &(mp->grid.axisvec1[1]),
1869 &(mp->grid.axisvec1[2]));
1870 free(line);
1871
1872 line = read_line(in);
1873 ptr = strchr(line, '=') + 1;
1874 sscanf (ptr, "%lf %lf %lf", &(mp->grid.axisvec2[0]),
1875 &(mp->grid.axisvec2[1]),
1876 &(mp->grid.axisvec2[2]));
1877 free(line);
1878
1879 line = read_line(in);
1880 ptr = strchr(line, '=') + 1;
1881 sscanf (ptr, "%lf %lf %lf", &(mp->grid.axisvec3[0]),
1882 &(mp->grid.axisvec3[1]),
1883 &(mp->grid.axisvec3[2]));
1884 free(line);
1885
1886 for(i = 0; i < 3; i++)
1887 {
1888 mp->grid.origin[i] *= R529;
1889 mp->grid.axisvec1[i] *= R529;
1890 mp->grid.axisvec2[i] *= R529;
1891 mp->grid.axisvec3[i] *= R529;
1892 }
1893
1894 if(mp->grid.isCutOff)
1895 {
1896 ie0=mp->grid.npt[0];
1897 ie1=mp->grid.npt[1];
1898 ie2=mp->grid.npt[2];
1899
1900 icount=0;
1901
1902 for(i0 = 0; i0 < ie0; i0++)
1903 for(i1 = 0; i1 < ie1; i1++)
1904 for(i2=0; i2 < ie2; i2++)
1905 {
1906 pp[0]=(mp->grid.origin[0]+mp->grid.axisvec1[0]*i0/(ie0-1)+mp->grid.axisvec2[0]*i1/(ie1-1)+mp->grid.axisvec3[0]*i2/(ie2-1));
1907 pp[1]=(mp->grid.origin[1]+mp->grid.axisvec1[1]*i0/(ie0-1)+mp->grid.axisvec2[1]*i1/(ie1-1)+mp->grid.axisvec3[1]*i2/(ie2-1));
1908 pp[2]=(mp->grid.origin[2]+mp->grid.axisvec1[2]*i0/(ie0-1)+mp->grid.axisvec2[2]*i1/(ie1-1)+mp->grid.axisvec3[2]*i2/(ie2-1));
1909 ishow=0;
1910 for(i = 0; i < mp->natom; i++)
1911 {
1912 iaia=0;
1913 if(fabs(mp->xyz[i][0]-pp[0]) < mp->grid.CutOff*myR529) iaia++;
1914 if(fabs(mp->xyz[i][1]-pp[1]) < mp->grid.CutOff*myR529) iaia++;
1915 if(fabs(mp->xyz[i][2]-pp[2]) < mp->grid.CutOff*myR529) iaia++;
1916
1917 if(iaia==3) ishow=1;
1918 }
1919 if(ishow==1) mp->grid.ByteCutOff[icount]=1;
1920 else mp->grid.ByteCutOff[icount]=0;
1921 icount++;
1922 }
1923 }
1924
1925 /*file positions*/
1926 #ifdef EBUG
1927 printf("Read up to Axis\n"); fflush(stdout);
1928 #endif
1929
1930 ii = 0;
1931 for (i = 0; i < mp->ngrids; i++)
1932 {
1933 #ifdef EBUG
1934 printf("ACCESSING GRID.POS POINTER: %ld\n", (long) mp->grid.pos); fflush(stdout);
1935 #endif
1936
1937 /*THIS MIGHT BE AN ERROR nBlock might be 0!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
1938
1939 /* if(!mp->grid.npoints || !mp->grid.nBlock) return; */
1940 if(mp->grid.nBlock)
1941 {
1942 for(j = 0; j < mp->grid.nBlock-1; j++)
1943 {
1944 mp->grid.pos[i * mp->grid.nBlock + j] = ((long)mp->ngrids * (long) j + (long) i) * 8L * (long) mp->grid.block_size;
1945 /*long long int*/
1946 /* mp->grid.pos[i * mp->grid.nBlock + j] = (long long int) mp->ngrids * (long long int) j + (long long int) i;
1947 mp->grid.pos[i * mp->grid.nBlock + j] *= 8LL;
1948 mp->grid.pos[i * mp->grid.nBlock + j] *= (long long int) mp->grid.block_size;*/
1949
1950
1951 #ifdef EBUG
1952 printf("orbital position igrid:%d iblock:%d %ld = %ld\n", i, j, mp->grid.pos[i * mp->grid.nBlock + j], (long) (mp->ngrids * j + i) * 8 * mp->grid.block_size); fflush(stdout);
1953 /* printf("orbital position igrid:%d iblock:%d %lld = %lld = %ld\n", i, j, ((long long int) mp->ngrids * (long long int) j + (long long int) i) * 8LL * mp->grid.block_size,
1954 mp->grid.pos[i * mp->grid.nBlock + j],
1955 (long) (mp->ngrids * j + i) * 8 * mp->grid.block_size); fflush(stdout);*/
1956 #endif
1957 }
1958
1959 mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1] = (long) (mp->grid.nBlock-1) * (long) mp->grid.block_size;
1960 mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1] *= (long) (mp->ngrids - i);
1961 mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1] += (long) mp->grid.npoints * (long) i;
1962 mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1] *= 8LL;
1963
1964 /*long long int*/
1965 /* {
1966 long long int tmp1, tmp2;
1967 tmp1 = (long long int) mp->ngrids * (long long int) (mp->grid.nBlock-1);
1968 tmp1 *= 8LL;
1969 tmp1 *= (long long int) mp->grid.block_size;
1970
1971 tmp2 = (long long int) (mp->grid.nBlock-1);
1972 tmp2 *= (long long int) mp->grid.block_size;
1973 tmp2 = (long long int) mp->grid.npoints - tmp2;
1974 tmp2 *= 8LL;
1975 tmp2 *= (long long int) i;
1976 mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1] = tmp1 + tmp2;
1977 }*/
1978
1979 #ifdef EBUG
1980 printf("orbital position igrid:%d iblock:%d %ld = %ld\n", i, mp->grid.nBlock-1, mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1], (long) (mp->ngrids * (mp->grid.nBlock-1)) * 8 * mp->grid.block_size +
1981 (long) (mp->grid.npoints - (mp->grid.nBlock-1) * mp->grid.block_size) * 8 * i); fflush(stdout);
1982 /*long long int*/
1983 /* printf("orbital position igrid:%d iblock:%d %lld = %ld\n", i, mp->grid.nBlock-1, mp->grid.pos[i * mp->grid.nBlock + mp->grid.nBlock-1], (long) (mp->ngrids * (mp->grid.nBlock-1)) * 8 * mp->grid.block_size +
1984 (long) (mp->grid.npoints - (mp->grid.nBlock-1) * mp->grid.block_size) * 8 * i); fflush(stdout);*/
1985 #endif
1986 }
1987 }
1988
1989 /*
1990 {
1991
1992 {
1993 mp->grid.pos[i * mp->grid.nBlock + j] = 8 * i * mp->grid.npoints + 8 * j * mp->grid.block_size;
1994 }
1995 }*/
1996
1997 /* mp->orb_ending_position = (long) mp->ngrids * (long) mp->grid.nBloc * (long) 8 * (long) mp->grid.block_size;*/
1998 /*orbital index position*/
1999
2000
2001 line = read_line(in);
2002 ptr = strchr(line, '=') + 1;
2003
2004 mp->orb_ending_position = atol(ptr);
2005 /*long long int*/
2006 /* mp->orb_ending_position = atoll(ptr);*/
2007 free(line);
2008
2009 #ifdef EBUG
2010 printf("End-of-orbitals location\n");
2011 #endif
2012
2013 /*orbital data*/
2014 for(i = 0; i < mp->ngrids; i++)
2015 {
2016 line = read_line(in);
2017
2018 mp->edited[i] = 0;
2019 if (strcasestr(line,"GridName"))
2020 {
2021 mp->grid.titlesArr[i] = my_read_str_value((char*) strcasestr(line,"GridName"));
2022 if (strcasestr(line, "Orbital"))
2023 {
2024 mp->grid_type[i] = ORBITAL;
2025 if (strcasestr(line, "sym"))
2026 mp->grid_symmetry[i] = my_read_int_value((char*) strcasestr(line,"sym"));
2027 else
2028 mp->grid_symmetry[i] = 0;
2029
2030 if (strcasestr(line, "index"))
2031 mp->grid_index[i] = my_read_int_value((char*) strcasestr(line,"index"));
2032 else
2033 mp->grid_index[i] = 0; /*some internal numbering should be put here*/
2034
2035 if (strcasestr(line, "Energ"))
2036 mp->grid_energy[i] = my_read_dble_value((char*) strcasestr(line,"energ"));
2037 else
2038 mp->grid_energy[i] = 0.0;
2039
2040 if (strcasestr(line, "OCc"))
2041 mp->grid_occ[i] = my_read_dble_value((char*) strcasestr(line,"occ"));
2042 else
2043 mp->grid_occ[i] = 0.0;
2044
2045 if (strcasestr(line, "type"))
2046 {
2047 for(ptr = (char*) strcasestr(line, "type"); ptr[0] != '='; ptr++);
2048 ptr++;
2049 for(; ptr[0] == 32; ptr++);
2050 switch (ptr[0])
2051 {
2052 case 'u': mp->orbital_type[i] = ORBITAL_TYPE_U; break;
2053 case 'f': mp->orbital_type[i] = ORBITAL_TYPE_F; break;
2054 case 'i': mp->orbital_type[i] = ORBITAL_TYPE_I; break;
2055 case '1': mp->orbital_type[i] = ORBITAL_TYPE_1; break;
2056 case '2': mp->orbital_type[i] = ORBITAL_TYPE_2; break;
2057 case '3': mp->orbital_type[i] = ORBITAL_TYPE_3; break;
2058 case 's': mp->orbital_type[i] = ORBITAL_TYPE_S; break;
2059 case 'd': mp->orbital_type[i] = ORBITAL_TYPE_D; break;
2060 default: mp->orbital_type[i] = ORBITAL_TYPE_U;
2061 }
2062 }
2063 else mp->orbital_type[i] = ORBITAL_TYPE_U;
2064 }
2065 else if (strcasestr(line, "Electrostatic"))
2066 {
2067 mp->ishow ^= HAS_E_POT;
2068 mp->grid_type[i] = E_POTEN;
2069
2070 mp->orbital_type[i] = ORBITAL_TYPE_U;
2071 mp->grid_symmetry[i] = 0;
2072 mp->grid_index[i] = 0;
2073 mp->grid_energy[i] = 0.0;
2074 mp->grid_occ[i] = 0.0;
2075 }
2076 else
2077 {
2078 mp->grid_type[i] = CUSTOM;
2079
2080 mp->orbital_type[i] = ORBITAL_TYPE_U;
2081 mp->grid_type[i] = CUSTOM;
2082 mp->grid_symmetry[i] = 0;
2083 mp->grid_index[i] = 0;
2084 mp->grid_energy[i] = 0.0;
2085 mp->grid_occ[i] = 0.0;
2086 }
2087 }
2088 /* else if (strcasestr(line, "Density")) mp->grid_type[i] = DENSITY;
2089 else if (strcasestr(line, "E_potential")) mp->grid_type[i] = E_POTEN;
2090 else mp->grid_type[i] = CUSTOM;*/
2091 free(line);
2092 }
2093
2094 /* line = read_line(in);
2095 free(line);*/
2096 #ifdef EBUG
2097 printf("Reading <DENSITY>\n"); fflush(stdout);
2098 #endif
2099
2100 while(next)
2101 {
2102 line = read_line(in);
2103 if (strcasestr(line, "density")) next = 0;
2104 if (strcasestr(line, "/grid") || feof(in))
2105 {
2106 make_warning("Could not find <DENSITY> section in luscus file. You might be attemting to open the old version of luscus file. Blame Valera Veryazov for changing luscus files without backward compatibility! Use old version of luscus or update your molcas.");
2107 deallocate_grids(mp);
2108 m->ngrids=0;
2109 free(line);
2110 return;
2111 }
2112 free(line);
2113 }
2114 next = 1;
2115
2116 iorb = mp->ngrids - 1;
2117
2118 /* find where <DENSITY> ends */
2119
2120 mp->orb_starting_position = ftell(in);
2121
2122 #ifdef EBUG
2123 printf("Position Set!\n"); fflush(stdout);
2124
2125 for(i = 0; i < mp->ngrids; i++)
2126 {
2127 printf("orb #%d type = %d E = %f occ = %f symm = %d index = %d edit = %d type = %d\n", i,
2128 mp->grid_type[i], mp->grid_energy[i], mp->grid_occ[i], mp->grid_symmetry[i],
2129 mp->grid_index[i], mp->edited[i], mp->orbital_type[i]);
2130 }
2131 #endif
2132
2133 /* if(sizeof(long int) == 4)
2134 {
2135 if (fread(&rec0, sizeof(int), 1, in)!=1) {fprintf(stderr, "READ ERROR\n"); fclose(in); return;}
2136 if (fread(buf,(unsigned int)rec0,1,in)!=1) {fprintf(stderr, "READ ERROR\n"); fclose(in); return;}
2137 if (fread(&rec1, sizeof(int), 1, in)!=1) {fprintf(stderr, "READ ERROR\n"); fclose(in); return;}
2138 if (rec0!=rec1) Fmarker=1;
2139 }*/
2140
2141 /*rewind to the end of the grid section*/
2142
2143 /* fseek(in, mp->grid.pos[(mp->ngrids - 1) * (mp->grid.nBlock - 1)], SEEK_CUR);*/
2144
2145 for(i = 0; i < mp->ngrids; i++)
2146 if (mp->grid_type[i] == E_POTEN)
2147 read_epoten_grid(in, mp, i);
2148
2149 fseek(in, mp->orb_ending_position, SEEK_CUR);
2150
2151 #ifdef EBUG
2152 printf("Jump to the Orbitals!\n");
2153 #endif
2154
2155 while(next)
2156 {
2157 line = read_line(in);
2158 if (strcasestr(line, "/grid")) next = 0;
2159 free(line);
2160 if (feof(in)) next = 0;
2161 }
2162 #ifdef EBUG
2163 printf("Return from read grid!\n");
2164 #endif
2165
2166 return;
2167
2168 }
2169
read_epoten_grid(FILE * in,MOL * mp,int i_grid)2170 void read_epoten_grid(FILE *in, MOL *mp, int i_grid)
2171 {
2172 int i;
2173 int iblock;
2174 int np;
2175 double *readbuffer;
2176 int success;
2177 int itek = 0;
2178 #ifdef EBUG
2179 FILE *tst;
2180 #endif
2181
2182 /*
2183 char *line;
2184 char word[50];
2185 int next = 1;
2186 int i, j;
2187 int x, y, z;
2188 long staring_position;
2189 double bestiso =0.1;
2190 double g;
2191 */
2192
2193 if (in == NULL)
2194 {
2195 make_warning("ERROR in reading Epot: Can't open file containig grid data!");
2196 /*register the existence of epot somewhere if possible*/
2197 return;
2198 }
2199
2200 isBinary = 1;
2201 isPacked = 0;
2202
2203 if (i_grid > m->ngrids || i_grid < 0)
2204 {
2205 make_warning("ERROR in reading Epot: Grid number out of range!");
2206 /*register the existence of epot somewhere if possible*/
2207 return;
2208 }
2209
2210 /*allocating data for readbuffer*/
2211 readbuffer = (double *) calloc ((unsigned int) m->grid.block_size, sizeof (double));
2212 /*allocating data for epot*/
2213 if (m->epot == NULL)
2214 m->epot = (double*) malloc(m->grid.npoints * sizeof(double));
2215
2216 iblock = 0;
2217 /* x = y = z = 0;*/
2218 for (i = 0; i < m->grid.npoints; i++) /* loop thro' file */
2219 {
2220 if (i == m->grid.block_size * iblock)
2221 {
2222 #ifdef EBUG
2223 printf("ACCESSING GRID.POS POINTER: %ld\n", (long) m->grid.pos); fflush(stdout);
2224 #endif
2225 fseek(in, m->orb_starting_position + m->grid.pos[i_grid * m->grid.nBlock + iblock], SEEK_SET);
2226 iblock++;
2227 np = m->grid.block_size;
2228 if (iblock == m->grid.nBlock)
2229 np = m->grid.npoints - m->grid.block_size * (m->grid.nBlock - 1);
2230 if(m->grid.isCutOff)
2231 {
2232 /* calculate number of uncutted data in this block */
2233 int ic = 0, ii;
2234 for(ii = 0; ii < np; ii++)
2235 if(m->grid.ByteCutOff[ii+i] == 1) ic++;
2236 success = ReadBlock(readbuffer, ic, in);
2237 }
2238 else
2239 {
2240 success = ReadBlock(readbuffer, np*sizeof(double), in);
2241 }
2242 itek = 0;
2243 }
2244 if(m->grid.isCutOff && m->grid.ByteCutOff[i]==0)
2245 {
2246 m->epot[i]=0.0;
2247 itek--;
2248 }
2249 else
2250 m->epot[i] = readbuffer[itek];
2251 itek++;
2252 /* For first point, set the min&max values;
2253 for subsequent points, conditionally set. */
2254 } /* endof loop */
2255
2256 free(readbuffer);
2257
2258 #ifdef EBUG
2259 tst = fopen("epot.txt", "w");
2260
2261 for(i = 0; i < m->grid.npoints; i++)
2262 {
2263 fprintf(tst, "%d %f\n", i, m->epot[i]);
2264 }
2265
2266 fclose(tst);
2267 #endif
2268
2269 }
2270
read_vibration_block(FILE * in,MOL * mp)2271 void read_vibration_block(FILE *in, MOL *mp)
2272 {
2273 int i;
2274 char *line;
2275 int next = 1;
2276 int nvibs = mp->nvibr + 1;
2277 allocate_vibs(mp, nvibs);
2278 line = read_line(in);
2279
2280 mp->freq[nvibs-1] = 0.0;
2281 mp->ir_intensity[nvibs-1] = 0.0;
2282 mp->raman_intensity[nvibs-1] = 0.0;
2283
2284 if (strcasestr(line,"freq"))
2285 mp->freq[nvibs-1] = my_read_dble_value((char*) strcasestr(line,"freq"));
2286 if (strcasestr(line,"ir_int"))
2287 {
2288 mp->ir_intensity[nvibs-1] = my_read_dble_value((char*) strcasestr(line,"ir_int"));
2289 if (!(mp->ishow & HAS_IRINT)) mp->ishow ^= HAS_IRINT;
2290 }
2291 if (strcasestr(line,"raman_int"))
2292 {
2293 mp->raman_intensity[nvibs-1] = my_read_dble_value((char*) strcasestr(line,"raman_int"));
2294 if (!(mp->ishow & HAS_RAMAN)) mp->ishow ^= HAS_RAMAN;
2295 }
2296
2297 free(line);
2298
2299 for(i = 0; i < mp->natom; i++)
2300 {
2301 line = read_line(in);
2302 sscanf(line, "%lf %lf %lf", &mp->normal_mode[nvibs-1][i][0], &mp->normal_mode[nvibs-1][i][1], &mp->normal_mode[nvibs-1][i][2]);
2303 free(line);
2304 }
2305 }
2306
read_editable_block(FILE * in,MOL * mp)2307 void read_editable_block(FILE *in, MOL *mp)
2308 {
2309 char *line;
2310 line = read_line(in);
2311 if (strcasestr(line,"yes")) mp->editable=1;
2312 else if (strcasestr(line,"no")) mp->editable=0;
2313 free(line);
2314 }
2315
read_sleep_block(FILE * in)2316 void read_sleep_block(FILE *in)
2317 {
2318 char *line;
2319 int seconds;
2320
2321 line = read_line(in);
2322 sscanf(line, "%d", &seconds);
2323 free(line);
2324 /* redraw();
2325 do_sleep(seconds);
2326 redraw();*/ /*BUG reading file should be redesigned if this option should take effect!*/
2327 }
2328
read_write_block(FILE * in)2329 void read_write_block(FILE *in)
2330 {
2331 char *line = NULL;
2332 do
2333 {
2334 line = read_line(in);
2335 if (line == NULL) break;
2336 if (line[strlen(line)-1] == 10) line[strlen(line)-1] = 0;
2337 luscus_gtk_pop_message_from_statusbar2();
2338 luscus_gtk_push_message_to_statusbar2(line);
2339 free(line);
2340 } while(!strcasestr(line,"write") && feof(in));
2341 }
2342
close_file(void)2343 void close_file(void)
2344 {
2345 #ifdef EBUG
2346 printf("CLOSING FILE\n");
2347 #endif
2348 delete_all_orbitals_from_list();
2349 deallocate_mol_data();
2350 mol = new_mol(1); /*Allocate one geometry; user can draw new molecule!*/
2351 luscus_gtk_show_or_hide_widgets();
2352 free(input_file_name);
2353 input_file_name = strdup("new_file.lus");
2354 remove_backup();
2355 init_file_time = 0;
2356 }
2357
check_if_file_changed(void)2358 void check_if_file_changed(void)
2359 {
2360 char tmppath[512];
2361 time_t checked_file_time;
2362 long checked_file_ntime;
2363 #ifdef EBUG
2364 printf("CHECKING IF FILE %s CHANGED\n", input_file_name);
2365 #endif
2366 if (!init_file_time) return;
2367 checked_file_time = get_file_time(input_file_name, &checked_file_ntime);
2368 #ifdef EBUG
2369 printf("INIT FILE TIME: %ld : %ld CHECKED FILE TIME: %ld : %ld\n", init_file_time, init_file_ntime, checked_file_time, checked_file_ntime);
2370 #endif
2371 if (init_file_time != checked_file_time || init_file_ntime != checked_file_ntime)
2372 {
2373 #ifdef EBUG
2374 printf("THERE IS A CHANGE!\n");
2375 #endif
2376 strncpy(tmppath, input_file_name, 511);
2377 if (input_file_type) open_file(tmppath, input_file_type);
2378 else
2379 {
2380 close_file();
2381 open_file(tmppath, NULL);
2382 }
2383 rerender_3d();
2384 redraw();
2385 }
2386 }
2387
get_input_filename(void)2388 char *get_input_filename(void)
2389 {
2390 if (input_file_name) return input_file_name;
2391 else return "newfile.lus";
2392 }
2393
get_current_directory(void)2394 char *get_current_directory(void)
2395 {
2396 return current_directory;
2397 }
2398
nextfname(void)2399 char *nextfname(void)
2400 {
2401 int n;
2402 int i;
2403 char *new_filename;
2404
2405 if (input_file_name == NULL) return NULL;
2406
2407 n = strlen(input_file_name);
2408 for (n = strlen(input_file_name); n > 0 && input_file_name[n] != '.'; n--);
2409 if (n == 0) n = strlen(input_file_name);
2410
2411 new_filename = (char*) malloc(sizeof(char) * (n+5));
2412 new_filename[0] = 0;
2413 for(i = 0; i < n; i++) new_filename[i] = input_file_name[i];
2414 new_filename[n] = '_';
2415 new_filename[n+1] = 'n';
2416 new_filename[n+2] = 'e';
2417 new_filename[n+3] = 'w';
2418 new_filename[n+4] = 0;
2419 return new_filename;
2420 }
2421
make_luscus_file_name(char * original_filename)2422 char *make_luscus_file_name(char* original_filename)
2423 {
2424 int i;
2425 int extlen;
2426 char *ext;
2427 char *luscus_filename;
2428
2429 for(i = strlen(original_filename); i > 0 && original_filename[i] != '.'; i--, ext = original_filename+i);
2430 if (strlen(original_filename) - i > MAX_EXT_LENGTH) ext = NULL; /*extensions greater than MAX_EXT_LENGTH do not make sense!*/
2431
2432 if (ext == NULL)
2433 {
2434 luscus_filename = (char*) malloc(sizeof(char) * (strlen(original_filename) + strlen(input_filetypes[n_input_types-1].extension) + 2));
2435 extlen = 0;
2436
2437 luscus_filename[0] = 0;
2438 for(i = 0; i <= strlen(original_filename) - extlen; i++) luscus_filename[i] = original_filename[i];
2439 luscus_filename[i] = 0;
2440 strcat(luscus_filename, ".");
2441 strcat(luscus_filename, input_filetypes[0].extension);
2442 }
2443 else
2444 {
2445 luscus_filename = (char*) malloc(sizeof(char) * (strlen(original_filename) - strlen(ext) + strlen(input_filetypes[n_input_types-1].extension) + 2));
2446 extlen = strlen(ext);
2447
2448 luscus_filename[0] = 0;
2449 for(i = 0; i <= strlen(original_filename) - extlen; i++) luscus_filename[i] = original_filename[i];
2450 luscus_filename[i] = 0;
2451 strcat(luscus_filename, input_filetypes[0].extension);
2452 }
2453
2454 return luscus_filename;
2455 }
2456
get_new_section(void)2457 void get_new_section(void)
2458 {
2459 m = mol + igeo;
2460
2461 /*populate atom-list with new data*/
2462 deallocate_atom_list();
2463 insert_all_atoms_into_list();
2464 }
2465
set_element_data(MOL * mp,int iat)2466 void set_element_data(MOL *mp, int iat)
2467 {
2468 int i;
2469 for(i = 0; i < number_of_elements; i++)
2470 if (strcasecmp(mp->elem[iat].name, e[i].name) == 0) break;
2471
2472 mp->elem[iat].color[0] = e[i].color[0];
2473 mp->elem[iat].color[1] = e[i].color[1];
2474 mp->elem[iat].color[2] = e[i].color[2];
2475 mp->elem[iat].vdw_rad = e[i].vdw_rad;
2476 mp->elem[iat].bond_rad = e[i].bond_rad;
2477 mp->elem[iat].valency = e[i].valency;
2478 mp->elem[iat].periodic_pos_x = e[i].periodic_pos_x;
2479 mp->elem[iat].periodic_pos_y = e[i].periodic_pos_y;
2480 }
2481
atom_distance(MOL * mp,int iat1,int iat2)2482 double atom_distance(MOL *mp, int iat1, int iat2)
2483 {
2484 if (m == NULL) return 0;
2485 return sqrt((mp->xyz[iat1][0] - mp->xyz[iat2][0]) * (mp->xyz[iat1][0] - mp->xyz[iat2][0]) +
2486 (mp->xyz[iat1][1] - mp->xyz[iat2][1]) * (mp->xyz[iat1][1] - mp->xyz[iat2][1]) +
2487 (mp->xyz[iat1][2] - mp->xyz[iat2][2]) * (mp->xyz[iat1][2] - mp->xyz[iat2][2]));
2488 }
2489
add_bond(MOL * mp,int iat1,int iat2,int bond_order)2490 void add_bond(MOL *mp, int iat1, int iat2, int bond_order)
2491 {
2492 int i;
2493 int nbond;
2494 for(i = 0; i < mp->nbond; i++)
2495 if ((mp->bond[i].iat1 == iat1 && mp->bond[i].iat2 == iat2) ||
2496 (mp->bond[i].iat1 == iat2 && mp->bond[i].iat2 == iat1))
2497 {
2498 mp->bond[i].bond_type = bond_order;
2499 return;
2500 }
2501
2502 /*bond is not recorded yet; make new record*/
2503 nbond = mp->nbond+1;
2504 allocate_bonds(mp, nbond);
2505 mp->bond[mp->nbond-1].iat1 = iat1;
2506 mp->bond[mp->nbond-1].iat2 = iat2;
2507 mp->bond[mp->nbond-1].bond_type = bond_order;
2508 }
2509
do_load_grid(void)2510 void do_load_grid(void)
2511 {
2512 read_grid_from_file();
2513 make_surfaces();
2514 }
2515
read_grid_from_file(void)2516 void read_grid_from_file(void)
2517 {
2518 FILE *in;
2519 char *line;
2520 char *luscus_filename;
2521 char word[50];
2522 int next = 1;
2523 int i, j;
2524 int np;
2525 int x, y, z;
2526 int itek = 0;
2527 int iblock;
2528 int success;
2529 long staring_position;
2530 double *readbuffer;
2531 double bestiso =0.1;
2532 double g;
2533
2534 if (!input_file_name)
2535 {
2536 make_warning("ERROR: File containing grid data is undefined!");
2537 deallocate_grids(m);
2538 return;
2539 }
2540
2541 luscus_filename = make_luscus_file_name(input_file_name);
2542
2543 in = fopen(luscus_filename, "rb");
2544 free(luscus_filename);
2545 if (in == NULL)
2546 {
2547 make_warning("ERROR: Can't open file containig grid data!");
2548 deallocate_grids(m);
2549 return;
2550 }
2551
2552 isBinary = 1;
2553 isPacked = 0;
2554
2555 if (iorb > m->ngrids || iorb < 0)
2556 {
2557 make_warning("ERROR: Grid number out of range!");
2558 deallocate_grids(m);
2559 return;
2560 }
2561
2562 /*allocating data for readbuffer*/
2563
2564 readbuffer = (double *) calloc ((unsigned int) m->grid.block_size, sizeof (double));
2565
2566 iblock = 0;
2567 x = y = z = 0;
2568 for (i = 0; i < m->grid.npoints; i++, x++) /* loop thro' file */
2569 {
2570 if (i == m->grid.block_size * iblock)
2571 {
2572 /*#ifdef EBUG
2573 printf("seeking position %ld\n", m->orb_starting_position + m->grid.pos[iorb * m->grid.nBlock + iblock]);
2574 #endif*/
2575 #ifdef EBUG
2576 printf("ACCESSING GRID.POS POINTER: %ld\n", (long) m->grid.pos); fflush(stdout);
2577 #endif
2578 fseek(in, m->orb_starting_position + m->grid.pos[iorb * m->grid.nBlock + iblock], SEEK_SET);
2579 iblock++;
2580 np = m->grid.block_size;
2581 if (iblock == m->grid.nBlock)
2582 np = m->grid.npoints - m->grid.block_size * (m->grid.nBlock - 1);
2583 if(m->grid.isCutOff)
2584 {
2585 /* calculate number of uncutted data in this block */
2586 int ic = 0, ii;
2587 for(ii = 0; ii < np; ii++)
2588 if(m->grid.ByteCutOff[ii+i] == 1) ic++;
2589 success = ReadBlock(readbuffer, ic, in);
2590 }
2591 else
2592 {
2593 success = ReadBlock(readbuffer, np*8, in);
2594 }
2595 itek = 0;
2596 }
2597 if (x == m->grid.npt[0])
2598 {
2599 x = 0;
2600 y++;
2601 if (y == m->grid.npt[1])
2602 {
2603 y = 0;
2604 z++;
2605 }
2606 }
2607 if(m->grid.isCutOff && m->grid.ByteCutOff[i]==0)
2608 {
2609 m->grid.values[i]=0.0;
2610 itek--;
2611 }
2612 else
2613 m->grid.values[i] = readbuffer[itek];
2614 itek++;
2615 /* For first point, set the min&max values;
2616 for subsequent points, conditionally set. */
2617 if (i == 0)
2618 {
2619 m->grid.minval = 10000.0;
2620 m->grid.maxval = -10000.0;
2621 }
2622 else
2623 {
2624 /* let's cut off huge values */
2625 if (m->grid.values[i] > 2) m->grid.values[i] = 2;
2626 m->grid.minval = MIN (m->grid.minval, m->grid.values[i]);
2627 m->grid.maxval = MAX (m->grid.maxval, m->grid.values[i]);
2628 }
2629 } /* endof loop */
2630
2631 {
2632 int iguess;
2633 double cutoff[21];
2634 double gguess[21];
2635 double mysmall;
2636 double mybig;
2637 double searchval;
2638 int isd=0;
2639 int ptotal;
2640 if (m->grid_type[iorb] == DENSITY) isd=1;
2641 /* if (strcasestr(m->grid.title,"Densit")) isd=1; */
2642 g = MAX(fabs(m->grid.maxval), fabs(m->grid.minval)) / 20;
2643
2644 ptotal = m->grid.npt[0] * m->grid.npt[1] * m->grid.npt[2];
2645 for (iguess = 0; iguess < 21; iguess++)
2646 {
2647 gguess[iguess] = 0;
2648 cutoff[iguess] = (iguess-.001)*g;
2649 if(iguess==0) cutoff[0]=0;
2650 for (i = 0; i < ptotal; i++)
2651 if(fabs(m->grid.values[i])>cutoff[iguess])
2652 {
2653 if(isd) gguess[iguess]+=m->grid.values[i];
2654 else gguess[iguess]+=m->grid.values[i]*m->grid.values[i];
2655 }
2656 }
2657
2658 for(i = 1; i < 21; i++)
2659 {
2660 searchval=gguess[0]*(i/20.);
2661 for (iguess = 0; iguess < 20; iguess++)
2662 {
2663 if(gguess[iguess]>=searchval && searchval>gguess[iguess+1])
2664 {
2665 mybig = gguess[iguess];
2666 mysmall = gguess[iguess+1];
2667 if(fabs(mybig-mysmall) < 1e-10) bestiso = .5*(cutoff[iguess]+cutoff[iguess+1]) ;
2668 else bestiso = ((searchval-mysmall)/(mybig-mysmall)) * (cutoff[iguess+1]-cutoff[iguess]) + cutoff[iguess];
2669 }
2670 }
2671 if(bestiso > m->grid.maxval) bestiso = -bestiso ;
2672 m->grid.isod[i]=bestiso;
2673 }
2674 }
2675
2676 m->grid.guess = m->grid.isod[14];
2677 sprintf(word, "%7.0le", m->grid.guess);
2678 sscanf(word, "%lf", &m->grid.guess);
2679 /* close the file and fly away */
2680
2681 free(readbuffer);
2682 fclose(in);
2683
2684 #ifdef NOCOMPILE
2685 #ifdef EBUG
2686 in = fopen("test.txt", "w");
2687
2688 for(i = 0; i < m->grid.npoints; i++)
2689 {
2690 fprintf(in, "%d %f\n", i, m->grid.values[i]);
2691 }
2692
2693 fclose(in);
2694 #endif
2695 #endif
2696 }
2697
2698 /***
2699 Read or skip block. When reading, unpack values if necessary.
2700
2701 In: fd -- file descriptor
2702 buf -- pointer to buffer to read into
2703 np -- expected number of values in block (buffer length, in doubles)
2704 Returns: success flag
2705 Notes:
2706 1) If pointer is NULL, does fake read (i.e. skips the block)
2707 2) After successfull reading/skipping the file pointer is positioned
2708 just after the block (presumably at the beginning of the next one)
2709 ***/
ReadBlock(double * buf,int np,FILE * fp)2710 static int ReadBlock(double* buf, int np, FILE* fp)
2711 {
2712 int j;
2713 char line[MAXCHAR];
2714 struct packing_params_t pk_prm;
2715 byte_t* packed_data;
2716
2717 if (fread(buf, (size_t) np, 1, fp) !=1) return -1;
2718
2719 return np;
2720 }
2721
2722 #ifdef NOCOMPILE
2723 /** Read block/record from unformatted Fortran file.
2724 Returns:
2725 -1 : on failure
2726 number of bytes
2727 read into buffer : on success
2728
2729 If buf==NULL, fakes read.
2730 ***/
read_bin_block(void * buf,int bufsz,FILE * fp)2731 int read_bin_block(void *buf, int bufsz, FILE * fp)
2732 {
2733 int rec_len, rec_len1, rec_8;
2734
2735 if (fread(&rec_len, sizeof(rec_len), 1, fp)!=1) return -1;
2736
2737 if(Fmarker==1)
2738 if (fread(&rec_8, sizeof(rec_8), 1, fp)!=1) return -1;
2739
2740 if (buf)
2741 {
2742 if (rec_len>bufsz) return -1;
2743 if (fread(buf,(unsigned int)rec_len,1,fp)!=1) return -1;
2744 }
2745 else
2746 fseek(fp,rec_len,SEEK_CUR);
2747
2748 if (fread(&rec_len1, sizeof(rec_len1), 1, fp)!=1) return -1;
2749
2750 if(Fmarker==1)
2751 if (fread(&rec_8, sizeof(rec_8), 1, fp)!=1) return -1;
2752
2753 if (rec_len1!=rec_len)
2754 {
2755 printf("ERROR reading BLOCK!\n"); fflush(stdout);
2756 return -1;
2757 }
2758 return rec_len;
2759 }
2760
2761
2762 /***
2763 Parse header, fill a structure with packing params.
2764 ***/
getPackingParams(const char * line,struct packing_params_t * prm)2765 static int getPackingParams(const char* line, struct packing_params_t* prm)
2766 {
2767 int flags;
2768 int i;
2769
2770 if (sscanf(line,"BHeader= %d %lf %lf %lf %lf %d %d %d",
2771 &flags, /* not used yet */
2772 &(prm->xlimit[0]), /* xmin */
2773 &(prm->xlimit[1]), /* xleft */
2774 &(prm->xlimit[2]), /* xright */
2775 &(prm->xlimit[3]), /* xmax */
2776 &(prm->ydelta[0]), /* dy1 */
2777 &(prm->ydelta[1]), /* dy2 */
2778 &(prm->ydelta[2]) /* dy3 */)!=8) return 0;
2779
2780 prm->ileft=0; prm->iright=2;
2781 if (prm->ydelta[0]==0)
2782 {
2783 prm->ileft++;
2784 prm->xdelta[0]=0.;
2785 }
2786 if (prm->ydelta[2]==0)
2787 {
2788 prm->iright--;
2789 prm->xdelta[2]=0.;
2790 }
2791
2792 prm->ylimit[0]=0;
2793 prm->ylimit[1]=prm->ydelta[0];
2794 prm->ylimit[2]=prm->ylimit[1] + prm->ydelta[1];
2795 prm->ylimit[3]=prm->ylimit[2] + prm->ydelta[2];
2796 prm->range=prm->ylimit[3];
2797
2798 prm->nbytes=(prm->range==128)?1:2;
2799
2800 for (i=prm->ileft; i<=prm->iright; i++)
2801 prm->xdelta[i]=prm->xlimit[i+1] - prm->xlimit[i];
2802
2803 fprintf(stderr,"LIMITS:\n"
2804 "range=%d\nxlimit=%f %f %f %f\nxdelta=%f %f %f\n"
2805 "ylimit=%d %d %d %d\nydelta=%d %d %d\n",
2806 prm->range,
2807 prm->xlimit[0],prm->xlimit[1],prm->xlimit[2],prm->xlimit[3],
2808 prm->xdelta[0],prm->xdelta[1],prm->xdelta[2],
2809 prm->ylimit[0],prm->ylimit[1],prm->ylimit[2],prm->ylimit[3],
2810 prm->ydelta[0],prm->ydelta[1],prm->ydelta[2]);
2811 return 1;
2812 }
2813 #endif
2814
abfgets(int isBinary,char * str,int max,FILE * fp)2815 int abfgets (int isBinary, char *str, int max, FILE * fp)
2816 {
2817 int i = 1;
2818 if (isBinary) i = read_bin_stream (str, max, fp);
2819 else if (fgets (str, max, fp) == NULL) i = 0;
2820 return i;
2821 }
2822
2823 /***
2824 Unpack data either from array of bytes or from integer value.
2825 In: dest -- array to unpack into
2826 packed_data -- packed values, represented as 1-byte or MSB-first 2-bytes
2827 packed_val -- alternatively, single integer to unpack
2828 ndata -- length (number of values, not size in bytes!) of packed_data,
2829 or, the same, length (in doubles) of the dest buffer
2830 prm -- packing parameters
2831
2832 ***/
unpackData(double * dest,byte_t * packed_data,unsigned int packed_val,int ndata,const struct packing_params_t * prm)2833 static void unpackData(double* dest,
2834 byte_t* packed_data, unsigned int packed_val,
2835 int ndata,
2836 const struct packing_params_t* prm)
2837 {
2838 if (packed_data==0) { ndata=1; } /* we're called to unpack single value */
2839
2840 while (ndata--)
2841 {
2842 if (packed_data)
2843 { /* we're called to unpack array */
2844 packed_val=*(packed_data++);
2845 if (prm->nbytes==2)
2846 { /* 2 bytes, MSB first!!! */
2847 packed_val<<=8;
2848 packed_val += *(packed_data++);
2849 }
2850 } /* else use supplied single value of packed_val */
2851
2852 /** Unpacking begins. **/
2853 /** Inlining -- for effectiveness. Separate block -- for clarity **/
2854 {
2855 int i;
2856 double x=0.0;
2857 int y=packed_val; /* for brevity */
2858 int minus=0;
2859
2860 y-=prm->range; /* BEWARE: reperesentable int range and sign! */
2861 if (y<0) { y=-y; minus=1; }
2862
2863 if (y<prm->ylimit[0]) { x=prm->xlimit[0]; goto x_done; }
2864 if (y>=prm->ylimit[3]) { x=prm->xlimit[3]; goto x_done; }
2865
2866 for (i=prm->ileft; i<=prm->iright; i++)
2867 {
2868 if (y>=prm->ylimit[i] && y<prm->ylimit[i+1])
2869 {
2870 x=((double)(y - prm->ylimit[i]))/ prm->ydelta[i] * prm->xdelta[i] + prm->xlimit[i];
2871 goto x_done;
2872 }
2873 }
2874 x_done:
2875 *(dest)=(minus?-x:x);
2876 if( *(dest) > prm->xlimit[3]) *(dest)=0;
2877 *(dest++);
2878 }
2879 /** Unpacking ends. **/
2880 } /* cycle through array ends */
2881 }
2882
2883 /** Read string from unformatted Fortran file
2884 Returns:
2885 0 : on failure
2886 1 : on success
2887 ***/
read_bin_stream(char * str,int max,FILE * fp)2888 int read_bin_stream(char *str, int max, FILE * fp)
2889 {
2890 int j;
2891
2892 j = fread(str, (size_t) max-1, 1, fp);
2893
2894 /* j=read_bin_block(str,max-1,fp);*/
2895 if (j==-1) return 0;
2896
2897 str[j]=0;
2898 return 1;
2899 }
2900
current_grid_limits(double * min,double * max,double * guess)2901 void current_grid_limits(double* min, double* max, double* guess)
2902 {
2903 if (m->ngrids)
2904 {
2905 *min = m->grid.minval;
2906 *max = m->grid.maxval;
2907 *guess = m->grid.guess;
2908 }
2909 }
2910
current_grid_data(double ** data)2911 int current_grid_data(double **data)
2912 {
2913 if (m->ngrids)
2914 {
2915 *data = m->grid.values;
2916 return 1;
2917 }
2918 return 0;
2919 }
2920
get_fmarker(void)2921 int get_fmarker(void)
2922 {
2923 return Fmarker;
2924 }
2925
2926