1 /* GridCP.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19
20 /* See W. Tang et al J. Phys. Condens. Matter 21 (2009) 084204 */
21
22 #include "../../Config.h"
23 #ifdef ENABLE_OMP
24 #include <omp.h>
25 #endif
26 #include "../Display/GlobalOrb.h"
27 #include "../Display/StatusOrb.h"
28 #include "../Display/UtilsOrb.h"
29 #include "../Utils/Utils.h"
30 #include "../Utils/UtilsInterface.h"
31 #include "../Utils/AtomsProp.h"
32 #include "../Utils/Constants.h"
33 #include "../Utils/GabeditTextEdit.h"
34 #include "../Common/Windows.h"
35 #include "../Display/GLArea.h"
36 #include "../Display/AtomicOrbitals.h"
37 #include "../Display/Orbitals.h"
38 #include "../Display/ColorMap.h"
39 #include "../Display/GeomOrbXYZ.h"
40 #include "../Display/BondsOrb.h"
41 #include "../Display/GridCP.h"
42
43 /* Point with a volume number :
44 * i>0 : point of volume # i
45 * i<0 : the critical point of volume # -i
46 * i= : point not yet assigned
47 */
48 #define TOL 1e-12
49
50 /* the g_list_remove use thread */
51 /**************************************************************************/
myg_list_remove(GList * list,gconstpointer data)52 GList* myg_list_remove (GList *list, gconstpointer data)
53 {
54 GList *l;
55
56 l = list;
57 while (l)
58 {
59 if (l->data != data) l = l->next;
60 else
61 {
62 if (l->prev) l->prev->next = l->next;
63 if (l->next) l->next->prev = l->prev;
64 if (list == l) list = list->next;
65 g_free (l);
66 break;
67 }
68 }
69 return list;
70 }
71 /* the g_list_prepend use thread */
72 /**************************************************************************/
myg_list_prepend(GList * list,gpointer data)73 GList* myg_list_prepend (GList *list, gpointer data)
74 {
75 GList *new_list;
76
77 new_list = g_malloc(sizeof(GList));
78 new_list->data = data;
79 new_list->next = list;
80
81 if (list)
82 {
83 new_list->prev = list->prev;
84 if (list->prev) list->prev->next = new_list;
85 list->prev = new_list;
86 }
87 else new_list->prev = NULL;
88
89 return new_list;
90 }
91 /**************************************************************************/
myg_list_free(GList * list)92 void myg_list_free (GList *list)
93 {
94 GList* l = NULL;
95 GList* next = NULL;
96 for(l=list;l!=NULL;l=next)
97 {
98 next = l->next;
99 if(l) g_free(l);
100 }
101 }
102 /**************************************************************************/
destroyListOfPointIndex(GList * listPointIndex)103 void destroyListOfPointIndex(GList* listPointIndex)
104 {
105 GList* list = NULL;
106 for(list=listPointIndex;list!=NULL;list=list->next)
107 {
108 PointIndex* data=(PointIndex*)list->data;
109 if(data) g_free(data);
110 }
111 myg_list_free(listPointIndex);
112 }
113 /**************************************************************************/
newPointIndex(gint i,gint j,gint k)114 PointIndex* newPointIndex(gint i, gint j, gint k)
115 {
116 PointIndex* data=g_malloc(sizeof(PointIndex));
117 data->i = i;
118 data->j = j;
119 data->k = k;
120 return data;
121 }
122 /**************************************************************************/
destroyListOfCriticalPoints(GList * listCriticalPoint)123 void destroyListOfCriticalPoints(GList* listCriticalPoint)
124 {
125 GList* list = NULL;
126
127 for(list=listCriticalPoint;list!=NULL;list=list->next)
128 {
129 CriticalPoint* data=(CriticalPoint*)list->data;
130 if(data) g_free(data);
131 }
132 myg_list_free(listCriticalPoint);
133 }
134 /**************************************************************************/
newCriticalPoint(gint i,gint j,gint k,gint numV)135 CriticalPoint* newCriticalPoint(gint i, gint j, gint k, gint numV)
136 {
137 CriticalPoint* data=g_malloc(sizeof(CriticalPoint));
138 data->index[0] = i;
139 data->index[1] = j;
140 data->index[2] = k;
141 data->rank = 0;
142 data->signature = 0;
143 data->lambda[0] = 0;
144 data->lambda[1] = 0;
145 data->lambda[2] = 0;
146 data->integral = 0;
147 data->volume = 0;
148 data->nuclearCharge = 0;
149 data->numVolume = numV;
150 data->numCenter = 0;
151 return data;
152 }
153 /**************************************************************************/
computeGrad(GridCP * gridCP)154 static void computeGrad(GridCP* gridCP)
155 {
156 gint i,j,k;
157 Grid* grid = NULL;
158 Point5*** points = NULL;
159 gdouble dx,dy,dz;
160 gdouble drx,dry,drz;
161 gint i1,i2, j1,j2, k1,k2;
162 gdouble c;
163
164 if(!gridCP) return;
165 grid = gridCP->grid;
166 if(!grid) return;
167 points = grid->point;
168 progress_orb_txt(0,_("Computing of gradient on each point..., Please wait"),TRUE);
169 for(i=0;i< grid->N[0] ;i++)
170 {
171 i1 = i+1;
172 i2 = i-1;
173 if(i2<0) i2 = i;
174 if(i1>grid->N[0]-1) i1 = i;
175 for(j=0;j< grid->N[1] ;j++)
176 {
177 j1 = j+1;
178 j2 = j-1;
179 if(j2<0) j2 = j;
180 if(j1>grid->N[1]-1) j1 = j;
181 for(k=0;k< grid->N[2] ;k++)
182 {
183 k1 = k+1;
184 k2 = k-1;
185 if(k2<0) k2 = k;
186 if(k1>grid->N[2]-1) k1 = k;
187
188 dx = points[i1][j][k].C[0]-points[i2][j][k].C[0];
189 dy = points[i1][j][k].C[1]-points[i2][j][k].C[1];
190 dz = points[i1][j][k].C[2]-points[i2][j][k].C[2];
191 drx = sqrt(dx*dx+dy*dy+dz*dz);
192 gridCP->grad[0][i][j][k] = (points[i1][j][k].C[3]-points[i2][j][k].C[3])/drx;
193
194
195 dx = points[i][j1][k].C[0]-points[i][j2][k].C[0];
196 dy = points[i][j1][k].C[1]-points[i][j2][k].C[1];
197 dz = points[i][j1][k].C[2]-points[i][j2][k].C[2];
198 dry = sqrt(dx*dx+dy*dy+dz*dz);
199 gridCP->grad[1][i][j][k] = (points[i][j1][k].C[3]-points[i][j2][k].C[3])/dry;
200
201
202 dx = points[i][j][k1].C[0]-points[i][j][k2].C[0];
203 dy = points[i][j][k1].C[1]-points[i][j][k2].C[1];
204 dz = points[i][j][k1].C[2]-points[i][j][k2].C[2];
205 drz = sqrt(dx*dx+dy*dy+dz*dz);
206 gridCP->grad[2][i][j][k] = (points[i][j][k1].C[3]-points[i][j][k2].C[3])/drz;
207
208 if(
209 points[i1][j][k].C[3]<points[i][j][k].C[3] &&
210 points[i2][j][k].C[3]<points[i][j][k].C[3]
211 ) gridCP->grad[0][i][j][k] = 0;
212 if(
213 points[i][j1][k].C[3]<points[i][j][k].C[3] &&
214 points[i][j2][k].C[3]<points[i][j][k].C[3]
215 ) gridCP->grad[1][i][j][k] = 0;
216 if(
217 points[i][j][k1].C[3]<points[i][j][k].C[3] &&
218 points[i][j][k2].C[3]<points[i][j][k].C[3]
219 ) gridCP->grad[2][i][j][k] = 0;
220 /*
221 gridCP->grad[0][i][j][k] /= 2;
222 gridCP->grad[1][i][j][k] /= 2;
223 gridCP->grad[2][i][j][k] /= 2;
224 c = 1;
225 if(fabs(gridCP->grad[0][i][j][k]) >TOL)
226 {
227 dx = fabs((points[i1][j][k].C[0]-points[i2][j][k].C[0])/2/gridCP->grad[0][i][j][k]);
228 if(c>dx) c = dx;
229 }
230 if(fabs(gridCP->grad[1][i][j][k]) >TOL)
231 {
232 dy = fabs((points[i][j1][k].C[1]-points[i][j2][k].C[1])/2/gridCP->grad[1][i][j][k]);
233 if(c>dy) c = dy;
234 }
235 if(fabs(gridCP->grad[2][i][j][k]) >TOL)
236 {
237 dz = fabs((points[i][j][k1].C[2]-points[i][j][k2].C[2])/2/gridCP->grad[2][i][j][k]);
238 if(c>dz) c = dz;
239 }
240 if(c>0)
241 {
242 gridCP->grad[0][i][j][k] *= c;
243 gridCP->grad[1][i][j][k] *= c;
244 gridCP->grad[2][i][j][k] *= c;
245 }
246 */
247 c= fabs(gridCP->grad[0][i][j][k]);
248 if(c < fabs(gridCP->grad[1][i][j][k])) c = fabs(gridCP->grad[1][i][j][k]);
249 if(c < fabs(gridCP->grad[2][i][j][k])) c = fabs(gridCP->grad[2][i][j][k]);
250 if(c>0)
251 {
252 c = 1.0/c;
253 gridCP->grad[0][i][j][k] *= c;
254 gridCP->grad[1][i][j][k] *= c;
255 gridCP->grad[2][i][j][k] *= c;
256 }
257 /*
258 if(points[i][j][k].C[3]>TOL
259 && fabs(gridCP->grad[0][i][j][k])<TOL
260 && fabs(gridCP->grad[1][i][j][k])<TOL
261 && fabs(gridCP->grad[2][i][j][k])<TOL
262 )
263 printf("Grad(%d %d %d) = %f %f %f\n",i,j,k,gridCP->grad[0][i][j][k],gridCP->grad[1][i][j][k],gridCP->grad[2][i][j][k]);
264 */
265 }
266 }
267 }
268 }
269 /**************************************************************************/
resetKnown(GridCP * gridCP)270 static void resetKnown(GridCP* gridCP)
271 {
272 gint i,j,k;
273
274 if(!gridCP) return;
275 for(i=0;i< grid->N[0] ;i++)
276 for(j=0;j< grid->N[1] ;j++)
277 for(k=0;k< grid->N[2] ;k++)
278 gridCP->known[i][j][k] = 0;
279 }
280 /**************************************************************************/
initGridCP(GridCP * gridCP,Grid * grid,Grid * gridAux)281 static void initGridCP(GridCP* gridCP, Grid* grid, Grid* gridAux)
282 {
283 gint i,j,k;
284 gint c;
285
286 if(!gridCP) return;
287 gridCP->grid = grid;
288 gridCP->gridAux = gridAux;
289 gridCP->volumeNumberOfPoints = NULL;
290 gridCP->known = NULL;
291 for(c=0;c<3;c++) gridCP->grad[c] = NULL;
292 gridCP->criticalPoints = NULL;
293 gridCP->integral = 0;
294 gridCP->nuclearCharge = 0;
295 if(!grid) return;
296 if(gridAux)
297 {
298 for(c=0;c<3;c++)
299 if(grid->N[c] != gridAux->N[c])
300 {
301 printf(_("The Cube of the 2 grids should be equals\n"));
302 return;
303 }
304 }
305 for(c=0;c<3;c++) if(grid->N[c]<1) return;
306 gridCP->volumeNumberOfPoints = g_malloc( grid->N[0]*sizeof(gint**));
307 for(i=0;i< grid->N[0] ;i++)
308 {
309 gridCP->volumeNumberOfPoints[i] = g_malloc(grid->N[1]*sizeof(gint*));
310 for(j=0;j< grid->N[1] ;j++)
311 {
312 gridCP->volumeNumberOfPoints[i][j] = g_malloc(grid->N[2]*sizeof(gint));
313 for(k=0;k< grid->N[2] ;k++)
314 gridCP->volumeNumberOfPoints[i][j][k] = 0;
315 }
316 }
317 gridCP->known = g_malloc( grid->N[0]*sizeof(gint**));
318 for(i=0;i< grid->N[0] ;i++)
319 {
320 gridCP->known[i] = g_malloc(grid->N[1]*sizeof(gint*));
321 for(j=0;j< grid->N[1] ;j++)
322 {
323 gridCP->known[i][j] = g_malloc(grid->N[2]*sizeof(gint));
324 for(k=0;k< grid->N[2] ;k++)
325 gridCP->known[i][j][k] = 0;
326 }
327 }
328 for(c=0;c<3;c++)
329 {
330 gridCP->grad[c] = g_malloc( grid->N[0]*sizeof(gdouble**));
331 for(i=0;i< grid->N[0] ;i++)
332 {
333 gridCP->grad[c][i] = g_malloc(grid->N[1]*sizeof(gdouble*));
334 for(j=0;j< grid->N[1] ;j++)
335 {
336 gridCP->grad[c][i][j] = g_malloc(grid->N[2]*sizeof(gdouble));
337 for(k=0;k< grid->N[2] ;k++)
338 gridCP->grad[c][i][j][k] = 0;
339 }
340 }
341 }
342 gridCP->dv = 1;
343 if(grid)
344 {
345 Point5*** points = grid->point;
346 gdouble xx = points[1][0][0].C[0]-points[0][0][0].C[0];
347 gdouble yy = points[0][1][0].C[1]-points[0][0][0].C[1];
348 gdouble zz = points[0][0][1].C[2]-points[0][0][0].C[2];
349 gridCP->dv = fabs(xx*yy*zz);
350 }
351 computeGrad(gridCP);
352
353 }
354 /**************************************************************************/
destroyGridCP(GridCP * gridCP)355 void destroyGridCP(GridCP* gridCP)
356 {
357 gint i,j;
358 gint c;
359
360 if(!gridCP) return;
361
362 if(gridCP->volumeNumberOfPoints)
363 {
364 for(i=0;i< grid->N[0] ;i++)
365 {
366 if(gridCP->volumeNumberOfPoints[i])
367 {
368 for(j=0;j< grid->N[1] ;j++)
369 if(gridCP->volumeNumberOfPoints[i][j]) g_free(gridCP->volumeNumberOfPoints[i][j]);
370 g_free(gridCP->volumeNumberOfPoints[i]);
371 }
372 }
373 g_free(gridCP->volumeNumberOfPoints);
374 }
375 if(gridCP->known)
376 {
377 for(i=0;i< grid->N[0] ;i++)
378 {
379 if(gridCP->known[i])
380 {
381 for(j=0;j< grid->N[1] ;j++)
382 if(gridCP->known[i][j]) g_free(gridCP->known[i][j]);
383 g_free(gridCP->known[i]);
384 }
385 }
386 g_free(gridCP->known);
387 }
388 for(c=0;c<3;c++)
389 if(gridCP->grad[c])
390 {
391 for(i=0;i< grid->N[0] ;i++)
392 {
393 if(gridCP->grad[c][i])
394 {
395 for(j=0;j< grid->N[1] ;j++)
396 if(gridCP->grad[c][i][j]) g_free(gridCP->grad[c][i][j]);
397 g_free(gridCP->grad[c][i]);
398 }
399 }
400 g_free(gridCP->grad[c]);
401 }
402 destroyListOfCriticalPoints(gridCP->criticalPoints);
403 gridCP->criticalPoints = NULL;
404
405 }
406 /**************************************************************************/
setArroundTo(GridCP * gridCP,gint current[3],gboolean kn)407 static gint setArroundTo(GridCP* gridCP, gint current[3], gboolean kn)
408 {
409 gint ic,jc,kc;
410 gint i1,i2;
411 gint j1,j2;
412 gint k1,k2;
413 gint*** vP;
414 gint I[3];
415 gint J[3];
416 gint K[3];
417 gint i, j, k;
418 gint n = 0;
419
420 if(!gridCP) return 0;
421 if(!gridCP->grid) return 0;
422 vP = gridCP->volumeNumberOfPoints;
423
424 i = current[0];
425 j = current[1];
426 k = current[2];
427
428 i1 = i+1;
429 i2 = i-1;
430 if(i2<0) i2 = i;
431 if(i1>gridCP->grid->N[0]-1) i1 = i;
432 j1 = j+1;
433 j2 = j-1;
434 if(j2<0) j2 = j;
435 if(j1>gridCP->grid->N[1]-1) j1 = j;
436 k1 = k+1;
437 k2 = k-1;
438 if(k2<0) k2 = k;
439 if(k1>gridCP->grid->N[2]-1) k1 = k;
440
441 I[0] = i2;
442 I[1] = i;
443 I[2] = i1;
444
445 J[0] = j2;
446 J[1] = j;
447 J[2] = j1;
448
449 K[0] = k2;
450 K[1] = k;
451 K[2] = k1;
452
453 for(ic=0;ic<3;ic++)
454 for(jc=0;jc<3;jc++)
455 for(kc=0;kc<3;kc++)
456 {
457 if(ic==1 && jc==1 && kc ==1) continue;
458 if(gridCP->known[I[ic]][J[jc]][K[kc]] != 1 && vP[I[ic]][J[jc]][K[kc]]==vP[i][j][k])
459 {
460 gridCP->known[I[ic]][J[jc]][K[kc]] = kn;
461 n++;
462 }
463 }
464 return n;
465 }
466 /**************************************************************************/
isVolumeEdge(GridCP * gridCP,gint current[3])467 static gboolean isVolumeEdge(GridCP* gridCP, gint current[3])
468 {
469 gint ic,jc,kc;
470 gint i1,i2;
471 gint j1,j2;
472 gint k1,k2;
473 gint*** vP;
474 gint I[3];
475 gint J[3];
476 gint K[3];
477 gint i, j, k;
478
479 if(!gridCP) return FALSE;
480 if(!gridCP->grid) return FALSE;
481 vP = gridCP->volumeNumberOfPoints;
482
483 i = current[0];
484 j = current[1];
485 k = current[2];
486
487 i1 = i+1;
488 i2 = i-1;
489 if(i2<0) i2 = i;
490 if(i1>gridCP->grid->N[0]-1) i1 = i;
491 j1 = j+1;
492 j2 = j-1;
493 if(j2<0) j2 = j;
494 if(j1>gridCP->grid->N[1]-1) j1 = j;
495 k1 = k+1;
496 k2 = k-1;
497 if(k2<0) k2 = k;
498 if(k1>gridCP->grid->N[2]-1) k1 = k;
499
500 I[0] = i2;
501 I[1] = i;
502 I[2] = i1;
503
504 J[0] = j2;
505 J[1] = j;
506 J[2] = j1;
507
508 K[0] = k2;
509 K[1] = k;
510 K[2] = k1;
511
512
513 for(ic=0;ic<3;ic++)
514 for(jc=0;jc<3;jc++)
515 for(kc=0;kc<3;kc++)
516 if( 2!=gridCP->known[I[ic]][J[jc]][K[kc]] && abs(vP[i][j][k]) != abs(vP[I[ic]][J[jc]][K[kc]]))
517 {
518 return TRUE;
519 }
520
521 return FALSE;
522 }
523 /**************************************************************************/
isMax(GridCP * gridCP,gint current[3])524 static gboolean isMax(GridCP* gridCP, gint current[3])
525 {
526 gint ic,jc,kc;
527 gint i1,i2;
528 gint j1,j2;
529 gint k1,k2;
530 Point5 ***points = NULL;
531 gint I[3];
532 gint J[3];
533 gint K[3];
534 gint i, j, k;
535
536 if(!gridCP) return FALSE;
537 if(!gridCP->grid) return FALSE;
538 points = gridCP->grid->point;
539
540 i = current[0];
541 j = current[1];
542 k = current[2];
543
544 i1 = i+1;
545 i2 = i-1;
546 if(i2<0) i2 = i;
547 if(i1>gridCP->grid->N[0]-1) i1 = i;
548 j1 = j+1;
549 j2 = j-1;
550 if(j2<0) j2 = j;
551 if(j1>gridCP->grid->N[1]-1) j1 = j;
552 k1 = k+1;
553 k2 = k-1;
554 if(k2<0) k2 = k;
555 if(k1>gridCP->grid->N[2]-1) k1 = k;
556
557 I[0] = i2;
558 I[1] = i;
559 I[2] = i1;
560
561 J[0] = j2;
562 J[1] = j;
563 J[2] = j1;
564
565 K[0] = k2;
566 K[1] = k;
567 K[2] = k1;
568
569
570 for(ic=0;ic<3;ic++)
571 for(jc=0;jc<3;jc++)
572 for(kc=0;kc<3;kc++)
573 if(points[I[ic]][J[jc]][K[kc]].C[3]>points[I[1]][J[1]][K[1]].C[3]) return FALSE;
574
575 return TRUE;
576 }
577 /**************************************************************************/
nextPointOnGrid(GridCP * gridCP,gint current[3],gint next[3])578 static void nextPointOnGrid(GridCP* gridCP, gint current[3], gint next[3])
579 {
580 gint i,j,k;
581 gint ic,jc,kc;
582 gint i1,i2;
583 gint j1,j2;
584 gint k1,k2;
585 gdouble dx;
586 gdouble dy;
587 gdouble dz;
588 Point5 ***points = NULL;
589 gint im, jm, km;
590 gdouble rhoCenter;
591 gdouble rhoMax;
592 gdouble rho;
593 gint I[3];
594 gint J[3];
595 gint K[3];
596
597 if(!gridCP) return;
598 if(!gridCP->grid) return;
599 points = gridCP->grid->point;
600
601 i = current[0];
602 j = current[1];
603 k = current[2];
604
605 i1 = i+1;
606 i2 = i-1;
607 if(i2<0) i2 = i;
608 if(i1>gridCP->grid->N[0]-1) i1 = i;
609 j1 = j+1;
610 j2 = j-1;
611 if(j2<0) j2 = j;
612 if(j1>gridCP->grid->N[1]-1) j1 = j;
613 k1 = k+1;
614 k2 = k-1;
615 if(k2<0) k2 = k;
616 if(k1>gridCP->grid->N[2]-1) k1 = k;
617
618 I[0] = i2;
619 I[1] = i;
620 I[2] = i1;
621
622 J[0] = j2;
623 J[1] = j;
624 J[2] = j1;
625
626 K[0] = k2;
627 K[1] = k;
628 K[2] = k1;
629
630 /*
631 printf("I = %d %d %d\n", I[0],I[1],I[2]);
632 printf("J = %d %d %d\n", J[0],J[1],J[2]);
633 printf("K = %d %d %d\n", K[0],K[1],K[2]);
634
635 printf("index = %d %d %d\n", i, j, k);
636 */
637
638 im = 1;
639 jm = 1;
640 km = 1;
641 /*printf("%d %d %d rho = %lf\n",I[im],J[jm],K[km],points[I[im]][J[jm]][K[km]].C[3]);*/
642 rhoCenter = points[I[1]][J[1]][K[1]].C[3];
643 rhoMax = rhoCenter;
644 for(ic=0;ic<3;ic++)
645 for(jc=0;jc<3;jc++)
646 for(kc=0;kc<3;kc++)
647 {
648 /*printf("%d %d %d rho = %lf\n",I[ic],J[jc],K[kc],points[I[ic]][J[jc]][K[kc]].C[3]);*/
649 if(ic==1 && jc==1 && kc==1) continue;
650 if(gridCP->known[I[ic]][J[jc]][K[kc]] >1) continue;
651 rho =points[I[ic]][J[jc]][K[kc]].C[3];
652
653 dx =points[I[ic]][J[jc]][K[kc]].C[0]-points[I[1]][J[1]][K[1]].C[0];
654 dy =points[I[ic]][J[jc]][K[kc]].C[1]-points[I[1]][J[1]][K[1]].C[1];
655 dz =points[I[ic]][J[jc]][K[kc]].C[2]-points[I[1]][J[1]][K[1]].C[2];
656 rho = rhoCenter + (rho-rhoCenter)/sqrt(dx*dx+dy*dy+dz*dz);
657
658 if(rho>rhoMax)
659 {
660 rhoMax = rho;
661 im = ic;
662 jm = jc;
663 km = kc;
664 }
665 }
666 /*printf("indexNEW = %d %d %d\n", IM, JM, KM);*/
667
668 next[0] = I[im];
669 next[1] = J[jm];
670 next[2] = K[km];
671 }
672 /**************************************************************************/
nextPoint(GridCP * gridCP,gdouble deltaR[3],gint current[3],gint next[3])673 static void nextPoint(GridCP* gridCP, gdouble deltaR[3], gint current[3], gint next[3])
674 {
675 gdouble gradrl[3];
676 gint c;
677 gint i,j,k;
678 if(!gridCP) return;
679 if(!gridCP->grid) return;
680
681
682 /*
683 nextPointOnGrid(gridCP, current, next);
684 return;
685 */
686
687 i = current[0];
688 j = current[1];
689 k = current[2];
690
691 for(c=0;c<3;c++) gradrl[c] = gridCP->grad[c][i][j][k];
692 if((gint)rint(gradrl[0]) ==0 && (gint)rint(gradrl[1]) ==0 && (gint)rint(gradrl[2]) ==0)
693 {
694 if(isMax(gridCP,current))
695 {
696 for(c=0;c<3;c++) next[c] = current[c];
697 for(c=0;c<3;c++) deltaR[c] = 0.0;
698 return;
699 }
700 else
701 {
702 nextPointOnGrid(gridCP, current, next);
703 for(c=0;c<3;c++) deltaR[c] = 0.0;
704 }
705 }
706 else
707 {
708 for(c=0;c<3;c++) next[c] = current[c] + (gint)rint(gradrl[c]);
709 for(c=0;c<3;c++) deltaR[c] += gradrl[c]-(gint)rint(gradrl[c]);
710 for(c=0;c<3;c++) next[c] += (gint)rint(deltaR[c]);
711 for(c=0;c<3;c++) deltaR[c] -= (gint)rint(deltaR[c]);
712 for(c=0;c<3;c++) if(next[c]<0 ) next[c] = 0;
713 for(c=0;c<3;c++) if(next[c]>gridCP->grid->N[c]-1) next[c] = gridCP->grid->N[c]-1;
714
715 i = current[0];
716 j = current[1];
717 k = current[2];
718 gridCP->known[i][j][k] = 1;
719 i = next[0];
720 j = next[1];
721 k = next[2];
722 if(gridCP->known[i][j][k]==1)
723 {
724 nextPointOnGrid(gridCP, current, next);
725 for(c=0;c<3;c++) deltaR[c] = 0.0;
726 }
727
728 }
729 }
730 /**************************************************************************/
addSurroundingEqualPoints(GridCP * gridCP,gint current[3],GList * listOfVisitedPoints)731 static GList* addSurroundingEqualPoints(GridCP* gridCP, gint current[3], GList* listOfVisitedPoints)
732 {
733 gint i,j,k;
734 gint ic,jc,kc;
735 gint i1,i2;
736 gint j1,j2;
737 gint k1,k2;
738 Point5 ***points = NULL;
739 gint I[3];
740 gint J[3];
741 gint K[3];
742 gdouble rho0 = 0;
743 gdouble dRho = 0;
744
745 if(!gridCP) return listOfVisitedPoints;
746 if(!gridCP->grid) return listOfVisitedPoints;
747 points = gridCP->grid->point;
748
749 i = current[0];
750 j = current[1];
751 k = current[2];
752
753 i1 = i+1;
754 i2 = i-1;
755 if(i2<0) i2 = i;
756 if(i1>gridCP->grid->N[0]-1) i1 = i;
757 j1 = j+1;
758 j2 = j-1;
759 if(j2<0) j2 = j;
760 if(j1>gridCP->grid->N[1]-1) j1 = j;
761 k1 = k+1;
762 k2 = k-1;
763 if(k2<0) k2 = k;
764 if(k1>gridCP->grid->N[2]-1) k1 = k;
765
766 I[0] = i2;
767 I[1] = i;
768 I[2] = i1;
769
770 J[0] = j2;
771 J[1] = j;
772 J[2] = j1;
773
774 K[0] = k2;
775 K[1] = k;
776 K[2] = k1;
777
778 rho0 = points[I[1]][J[1]][K[1]].C[3];
779 for(ic=0;ic<3;ic++)
780 for(jc=0;jc<3;jc++)
781 for(kc=0;kc<3;kc++)
782 {
783 if(ic==1 && jc==1 && kc==1) continue;
784 dRho =points[I[ic]][J[jc]][K[kc]].C[3]-rho0;
785 if(fabs(dRho)<TOL)
786 {
787 PointIndex* data = newPointIndex( I[ic], J[jc], K[kc]);
788 listOfVisitedPoints = myg_list_prepend(listOfVisitedPoints,data);
789 gridCP->known[I[ic]][J[jc]][K[kc]] = 1;
790 }
791 }
792 return listOfVisitedPoints;
793 }
794 /**************************************************************************/
assentTrajectory(GridCP * gridCP,gint current[3],gboolean ongrid)795 static GList* assentTrajectory(GridCP* gridCP, gint current[3], gboolean ongrid)
796 {
797 GList* listOfVisitedPoints = NULL;
798 /*Point5 ***points = NULL;*/
799 gint next[3];
800 gdouble deltaR[3] = {0,0,0};
801 gint l;
802 gint imax;
803 PointIndex* data;
804 gint c;
805
806 if(!gridCP) return listOfVisitedPoints;
807 if(!gridCP->grid) return listOfVisitedPoints;
808 /* points = gridCP->grid->point;*/
809
810 for(c=0;c<3;c++) if(grid->N[c]<1) return listOfVisitedPoints;
811
812 data = newPointIndex(current[0], current[1], current[2]);
813 listOfVisitedPoints = myg_list_prepend(listOfVisitedPoints,data);
814
815 imax = grid->N[0]* grid->N[1]* grid->N[2];
816 for(l=0;l<imax;l++)
817 {
818
819 if(ongrid) nextPointOnGrid(gridCP, current, next);
820 else nextPoint(gridCP, deltaR, current, next);
821 data = newPointIndex( next[0], next[1], next[2]);
822 listOfVisitedPoints = myg_list_prepend(listOfVisitedPoints,data);
823
824 if(CancelCalcul) break;
825 if(next[0] == current[0] && next[1] == current[1] && next[2] == current[2])
826 {
827 /* new critical point */
828 /*printf("New critical point, index = %d %d %d\n",next[0] , next[1] , next[2] );*/
829 listOfVisitedPoints =addSurroundingEqualPoints(gridCP, current, listOfVisitedPoints);
830 break;
831 }
832 else
833 {
834 /*printf("cur = %d %d %d\n",next[0], next[1], next[2]);*/
835 }
836 /*if(vP[next[0]][next[1]][next[2]] != 0)
837 {
838 found a point from an old detected volume
839
840 printf("found a point from an old detected volume\n");
841 for(c=0;c<3;c++)current[c] = next[c];
842 if(isVolumeEdge(gridCP,current)) continue;
843 break;
844 }
845 */
846 for(c=0;c<3;c++)current[c] = next[c];
847 }
848 return listOfVisitedPoints;
849 }
850
851
852
853
854 /**************************************************************************/
assignGridCP(GridCP * gridCP,gboolean ongrid)855 static void assignGridCP(GridCP* gridCP, gboolean ongrid)
856 {
857 gint i;
858 gint*** vP = NULL;
859 Point5 ***points = NULL;
860 gint numberOfCriticalPoints = 0;
861 gchar* str =_("Assignation of points to volumes... Please wait");
862 gdouble scal;
863
864 if(!gridCP) return;
865 if(!gridCP->grid) return;
866 points = gridCP->grid->point;
867
868 for(i=0;i<3;i++) if(grid->N[i]<1) return;
869 vP = gridCP->volumeNumberOfPoints;
870 if(gridCP->criticalPoints)
871 {
872 destroyListOfCriticalPoints(gridCP->criticalPoints);
873 gridCP->criticalPoints = NULL;
874 }
875 progress_orb_txt(0,str,TRUE);
876
877 resetKnown(gridCP);
878 /*
879 for(i=0;i<grid->N[0];i++)
880 for(j=0;j<grid->N[1];j++)
881 for(k=0;k<grid->N[2];k++)
882 if(points[i][j][k].C[3]<TOL) gridCP->known[i][j][k] = 2;
883 */
884
885
886 scal = 1.1/(grid->N[0]-1);
887 #ifdef ENABLE_OMP
888 #pragma omp parallel for private(i)
889 #endif
890 for(i=0;i<grid->N[0];i++)
891 {
892 gint j,k;
893 gint current[3];
894 #ifdef ENABLE_OMP
895 #ifndef G_OS_WIN32
896 #pragma omp critical
897 progress_orb_txt(scal,str,FALSE);
898 #endif
899 #else
900 progress_orb_txt(scal,str,FALSE);
901 #endif
902 if(!CancelCalcul)
903 for(j=0;j<grid->N[1];j++)
904 {
905 if(!CancelCalcul)
906 for(k=0;k<grid->N[2];k++)
907 {
908 GList* listOfVisitedPoints = NULL;
909 current[0] = i;
910 current[1] = j;
911 current[2] = k;
912 if(vP[i][j][k] != 0) continue;
913 /*if(gridCP->known[i][j][k]!=0) continue;*/
914 if(points[i][j][k].C[3]<TOL) continue;
915
916 if(CancelCalcul) break;
917 listOfVisitedPoints = assentTrajectory(gridCP, current, ongrid);
918
919 #ifdef ENABLE_OMP
920 #pragma omp critical
921 #endif
922 if(vP[current[0]][current[1]][current[2]] != 0)
923 {
924 GList* list=NULL;
925 gint icp = abs(vP[current[0]][current[1]][current[2]]);
926 for(list=listOfVisitedPoints;list!=NULL;list=list->next)
927 {
928 PointIndex* data=(PointIndex*)list->data;
929 vP[data->i] [data->j] [data->k] = icp;
930 }
931 destroyListOfPointIndex(listOfVisitedPoints);
932 listOfVisitedPoints = NULL;
933 }
934 else
935 {
936 numberOfCriticalPoints++;
937 CriticalPoint* data = NULL;
938 GList* list=NULL;
939 gint icp = numberOfCriticalPoints;
940 for(list=listOfVisitedPoints;list!=NULL;list=list->next)
941 {
942 PointIndex* data=(PointIndex*)list->data;
943 vP[data->i] [data->j] [data->k] = icp;
944 }
945 destroyListOfPointIndex(listOfVisitedPoints);
946 listOfVisitedPoints = NULL;
947 vP[current[0]][current[1]][current[2]] = -icp;
948 data = newCriticalPoint(current[0], current[1], current[2],numberOfCriticalPoints);
949 gridCP->criticalPoints = myg_list_prepend(gridCP->criticalPoints,data);
950 }
951
952 }
953 }
954 }
955 progress_orb_txt(0," ",TRUE);
956 resetKnown(gridCP);
957 }
958 /**************************************************************************/
refineEdge(GridCP * gridCP,gboolean ongrid)959 static gint refineEdge(GridCP* gridCP, gboolean ongrid)
960 {
961 gint i,j,k;
962 gint*** vP = NULL;
963 Point5 ***points = NULL;
964 gboolean ***known = NULL;
965 gchar* str ="Refine grid points adjacent to Bader surface... Please wait";
966 gdouble scal;
967 gint current[3];
968 gint ne = 0;
969 GList* list=NULL;
970
971 if(!gridCP) return 0;
972 if(!gridCP->grid) return 0;
973 points = gridCP->grid->point;
974 known = gridCP->known;
975
976 for(i=0;i<3;i++) if(grid->N[i]<1) return 0;
977 vP = gridCP->volumeNumberOfPoints;
978 progress_orb_txt(0,str,TRUE);
979
980 resetKnown(gridCP);
981
982 for(i=0;i<grid->N[0];i++)
983 for(j=0;j<grid->N[1];j++)
984 for(k=0;k<grid->N[2];k++)
985 {
986 current[0] = i;
987 current[1] = j;
988 current[2] = k;
989 if(points[i][j][k].C[3]<TOL) known[i][j][k] = 2;
990 if(vP[i][j][k]<0) known[i][j][k] = 2;
991 }
992 for(i=0;i<grid->N[0];i++)
993 for(j=0;j<grid->N[1];j++)
994 for(k=0;k<grid->N[2];k++)
995 {
996 current[0] = i;
997 current[1] = j;
998 current[2] = k;
999 if(vP[i][j][k]>0 && !isMax(gridCP,current) && isVolumeEdge(gridCP,current))
1000 {
1001 ne++;
1002 setArroundTo(gridCP, current, 0);
1003 known[i][j][k] = 1;
1004 }
1005 else known[i][j][k] = 2;
1006 }
1007
1008 ne = 0;
1009 for(i=0;i<grid->N[0];i++)
1010 for(j=0;j<grid->N[1];j++)
1011 for(k=0;k<grid->N[2];k++)
1012 {
1013 if(known[i][j][k]==1)
1014 {
1015 ne++;
1016 vP[i][j][k] = 0;
1017 known[i][j][k]=0;
1018 }
1019 }
1020
1021 scal = 1.1/(grid->N[0]-1);
1022 #ifdef ENABLE_OMP
1023 #pragma omp parallel for private(i)
1024 #endif
1025 for(i=0;i<grid->N[0];i++)
1026 {
1027 gint current[3];
1028 gint j,k;
1029 GList* listOfVisitedPoints = NULL;
1030 #ifdef ENABLE_OMP
1031 #ifndef G_OS_WIN32
1032 #pragma omp critical
1033 progress_orb_txt(scal,str,FALSE);
1034 #endif
1035 #else
1036 progress_orb_txt(scal,str,FALSE);
1037 #endif
1038 if(!CancelCalcul)
1039 for(j=0;j<grid->N[1];j++)
1040 {
1041 if(CancelCalcul) break;
1042 for(k=0;k<grid->N[2];k++)
1043 {
1044 current[0] = i;
1045 current[1] = j;
1046 current[2] = k;
1047 if(vP[i][j][k]!=0) continue;
1048 if(gridCP->known[i][j][k] !=0 ) continue;
1049
1050 if(CancelCalcul) break;
1051 listOfVisitedPoints = assentTrajectory(gridCP, current, ongrid);
1052 vP[i][j][k] = abs(vP[current[0]][current[1]][current[2]]);
1053 #ifdef ENABLE_OMP
1054 #pragma omp critical
1055 #endif
1056 for(list=listOfVisitedPoints;list!=NULL;list=list->next)
1057 {
1058 PointIndex* data=(PointIndex*)list->data;
1059 known[data->i] [data->j] [data->k] = 0;
1060 }
1061 destroyListOfPointIndex(listOfVisitedPoints);
1062 }
1063 }
1064 }
1065 progress_orb_txt(0," ",TRUE);
1066
1067 resetKnown(gridCP);
1068 return ne;
1069 }
1070 /**************************************************************************/
assignPointsZero(GridCP * gridCP)1071 static void assignPointsZero(GridCP* gridCP)
1072 {
1073 GList* criticalPoint = gridCP->criticalPoints;
1074 gint*** vP = gridCP->volumeNumberOfPoints;
1075 GList* list = NULL;
1076 gint i,j,k;
1077 gdouble scal;
1078 Grid* grid = gridCP->grid;
1079 gchar* str = _("Assignation of points with f = 0..., Please wait");
1080 gdouble dx, dy, dz;
1081 gdouble r;
1082
1083 progress_orb_txt(0,str,TRUE);
1084 scal = 1.1/(grid->N[0]-1);
1085 for(i=0;i<grid->N[0];i++)
1086 {
1087 progress_orb_txt(scal,str,FALSE);
1088 if(CancelCalcul) break;
1089 for(j=0;j<grid->N[1];j++)
1090 for(k=0;k<grid->N[2];k++)
1091 {
1092 gdouble rmin = 0;
1093 CriticalPoint* dataMin = NULL;
1094 if(vP[i][j][k]!=0) continue;
1095 if(CancelCalcul) break;
1096 for(list=criticalPoint;list!=NULL;list=list->next)
1097 {
1098 CriticalPoint* data=(CriticalPoint*)list->data;
1099 gint ii = data->index[0];
1100 gint jj = data->index[1];
1101 gint kk = data->index[2];
1102 dx = gridCP->grid->point[i][j][k].C[0]-gridCP->grid->point[ii][jj][kk].C[0];
1103 dy = gridCP->grid->point[i][j][k].C[1]-gridCP->grid->point[ii][jj][kk].C[1];
1104 dz = gridCP->grid->point[i][j][k].C[2]-gridCP->grid->point[ii][jj][kk].C[2];
1105 r = (dx*dx + dy*dy + dz*dz);
1106 if(dataMin == NULL || r<rmin )
1107 {
1108 dataMin = data;
1109 rmin = r;
1110 }
1111 }
1112 if(dataMin)
1113 {
1114 gint ii = dataMin->index[0];
1115 gint jj = dataMin->index[1];
1116 gint kk = dataMin->index[2];
1117 vP[i][j][k] = abs(vP[ii][jj][kk]);
1118 }
1119 }
1120 }
1121 progress_orb_txt(0," ",TRUE);
1122 }
1123 /**************************************************************************/
addToResult(gchar * result,gchar * tmp)1124 static gchar* addToResult(gchar* result, gchar* tmp)
1125 {
1126 gchar* old = result;
1127 if(!tmp) return result;
1128 if(old)
1129 {
1130 result = g_strdup_printf("%s%s",old,tmp);
1131 g_free(old);
1132 }
1133 else result = g_strdup_printf("%s",tmp);
1134 return result;
1135 }
1136 /**************************************************************************/
removeAttractor0(GridCP * gridCP)1137 static void removeAttractor0(GridCP* gridCP)
1138 {
1139 GList* list=NULL;
1140 GList* next=NULL;
1141 gint*** vP = gridCP->volumeNumberOfPoints;
1142 for(list=gridCP->criticalPoints;list!=NULL;list=next)
1143 {
1144 CriticalPoint* data = (CriticalPoint*) list->data;
1145 next = list->next;
1146 gint i = data->index[0];
1147 gint j = data->index[1];
1148 gint k = data->index[2];
1149 if(gridCP->grid->point[i][j][k].C[3]<TOL)
1150 {
1151 vP[i][j][k] = 0;
1152 gridCP->criticalPoints=myg_list_remove(gridCP->criticalPoints, data);
1153 g_free(data);
1154 }
1155 }
1156 }
1157 /**************************************************************************/
removeNonSignificantAttractor(GridCP * gridCP)1158 static void removeNonSignificantAttractor(GridCP* gridCP)
1159 {
1160 GList* list=NULL;
1161 GList* next=NULL;
1162 gint n = gridCP->grid->N[0]*gridCP->grid->N[1]*gridCP->grid->N[2]/1000;
1163 for(list=gridCP->criticalPoints;list!=NULL;list=next)
1164 {
1165 CriticalPoint* data = (CriticalPoint*) list->data;
1166 next = list->next;
1167 if(data->volume/gridCP->dv<n)
1168 {
1169 /* printf("n = %d\n",(gint)(data->volume/gridCP->dv));*/
1170 gridCP->criticalPoints=myg_list_remove(gridCP->criticalPoints, data);
1171 g_free(data);
1172 }
1173 }
1174 }
1175 /************************************************************************************************************/
setPartialChargeToAIM(GtkWidget * win)1176 static void setPartialChargeToAIM(GtkWidget *win)
1177 {
1178 GridCP* gridCP = NULL;
1179 GList* list = NULL;
1180 if(GTK_IS_WIDGET(win)) gridCP = g_object_get_data(G_OBJECT (win), "GridCP");
1181 if(!gridCP) return;
1182 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1183 {
1184 CriticalPoint* data = (CriticalPoint*) list->data;
1185 gint c= data->numCenter;
1186 GeomOrb[c].partialCharge = data->nuclearCharge-data->integral;
1187 }
1188 glarea_rafresh(GLArea);
1189 }
1190 /************************************************************************************************************/
destroyResultDlg(GtkWidget * win)1191 static void destroyResultDlg(GtkWidget *win)
1192 {
1193 GridCP* gridCP = NULL;
1194 if(GTK_IS_WIDGET(win)) gridCP = g_object_get_data(G_OBJECT (win), "GridCP");
1195 if(gridCP)
1196 {
1197 destroyGridCP(gridCP);
1198 g_free(gridCP);
1199 }
1200 if(GTK_IS_WIDGET(win)) delete_child(win);
1201 if(GTK_IS_WIDGET(win)) gtk_widget_destroy(win);
1202 }
1203 /********************************************************************************/
showResultDlg(gchar * message,gchar * title,GridCP * gridCP)1204 static GtkWidget* showResultDlg(gchar *message,gchar *title,GridCP* gridCP)
1205 {
1206 GtkWidget *dlgWin = NULL;
1207 GtkWidget *frame;
1208 GtkWidget *vboxframe;
1209 GtkWidget *txtWid;
1210 GtkWidget *button;
1211
1212
1213 dlgWin = gtk_dialog_new();
1214 gtk_widget_realize(GTK_WIDGET(dlgWin));
1215
1216 gtk_window_set_title(GTK_WINDOW(dlgWin),title);
1217 gtk_window_set_position(GTK_WINDOW(dlgWin),GTK_WIN_POS_CENTER);
1218 gtk_window_set_modal (GTK_WINDOW (dlgWin), TRUE);
1219 gtk_window_set_transient_for(GTK_WINDOW(dlgWin),GTK_WINDOW(PrincipalWindow));
1220
1221 g_signal_connect(G_OBJECT(dlgWin), "delete_event", (GCallback)destroyResultDlg, NULL);
1222 frame = gtk_frame_new (NULL);
1223 gtk_frame_set_shadow_type( GTK_FRAME(frame),GTK_SHADOW_ETCHED_OUT);
1224
1225 gtk_container_set_border_width (GTK_CONTAINER (frame), 5);
1226 gtk_box_pack_start( GTK_BOX(GTK_DIALOG(dlgWin)->vbox), frame,TRUE,TRUE,0);
1227
1228 gtk_widget_show (frame);
1229
1230 vboxframe = create_vbox(frame);
1231 txtWid = create_text_widget(vboxframe,NULL,&frame);
1232 if(message) gabedit_text_insert (GABEDIT_TEXT(txtWid), NULL, NULL, NULL,message,-1);
1233
1234 gtk_box_set_homogeneous (GTK_BOX( GTK_DIALOG(dlgWin)->action_area), FALSE);
1235
1236 button = create_button(dlgWin,_("Partial charges of molecule <= AIM charges"));
1237 gtk_box_pack_end (GTK_BOX( GTK_DIALOG(dlgWin)->action_area), button, FALSE, TRUE, 5);
1238 GTK_WIDGET_SET_FLAGS(button, GTK_CAN_DEFAULT);
1239 gtk_widget_grab_default(button);
1240 g_signal_connect_swapped(G_OBJECT(button), "clicked", (GCallback)setPartialChargeToAIM, GTK_OBJECT(dlgWin));
1241
1242 button = create_button(dlgWin,"Close");
1243 gtk_box_pack_end (GTK_BOX( GTK_DIALOG(dlgWin)->action_area), button, FALSE, TRUE, 5);
1244 GTK_WIDGET_SET_FLAGS(button, GTK_CAN_DEFAULT);
1245 gtk_widget_grab_default(button);
1246 g_signal_connect_swapped(G_OBJECT(button), "clicked", (GCallback)destroyResultDlg, GTK_OBJECT(dlgWin));
1247
1248 add_button_windows(title,dlgWin);
1249 gtk_window_set_default_size (GTK_WINDOW(dlgWin), (gint)(ScreenHeight*0.6), (gint)(ScreenHeight*0.5));
1250 gtk_widget_show_all(dlgWin);
1251 g_object_set_data(G_OBJECT (dlgWin), "GridCP",gridCP);
1252 return dlgWin;
1253 }
1254 /**************************************************************************/
showGridCP(GridCP * gridCP)1255 static void showGridCP(GridCP* gridCP)
1256 {
1257 GList* list=NULL;
1258 gchar* result =NULL;
1259 gchar* tmp = NULL;
1260 gint nc = 0;
1261 gdouble xx, yy, zz;
1262 gint n0 = 0, n1 = 0, n2 = 0;
1263 gint i,j,k;
1264 gdouble sum = 0;
1265 Point5*** points = gridCP->grid->point;
1266
1267 result = addToResult(result, _("Geometry (Ang)\n"));
1268 result = addToResult(result, "==============\n");
1269 for(j=0; j<(gint)nCenters; j++)
1270 {
1271 tmp = g_strdup_printf("%s[%d] %lf %lf %lf\n",
1272 GeomOrb[j].Prop.symbol,
1273 j+1,
1274 GeomOrb[j].C[0]*BOHR_TO_ANG,
1275 GeomOrb[j].C[1]*BOHR_TO_ANG,
1276 GeomOrb[j].C[2]*BOHR_TO_ANG);
1277 result = addToResult(result, tmp);
1278 if(tmp) g_free(tmp);
1279 }
1280 result = addToResult(result, "---------------------------------------------------------------------\n");
1281
1282 n0 = gridCP->grid->N[0];
1283 n1 = gridCP->grid->N[1];
1284 n2 = gridCP->grid->N[2];
1285
1286 xx = points[n0-1][0][0].C[0]-points[0][0][0].C[0];
1287 yy = points[0][n1-1][0].C[1]-points[0][0][0].C[1];
1288 zz = points[0][0][n2-1].C[2]-points[0][0][0].C[2];
1289 tmp = g_strdup_printf(
1290 _(
1291 "Grid point density (Ang^-1) on the first direction(>10 is recommended) = %lf\n"
1292 "density of the grid(Ang^-1) on the second direction(>10 is recommended) = %lf\n"
1293 "density of the grid(Ang^-1) on the third direction(>10 is recommended) = %lf\n"
1294 ),
1295 n0/(xx*BOHR_TO_ANG),
1296 n1/(yy*BOHR_TO_ANG),
1297 n2/(zz*BOHR_TO_ANG)
1298 );
1299 result = addToResult(result, tmp);
1300 if(tmp) g_free(tmp);
1301 sum = 0;
1302 for(i=0;i<n0 ;i+=n0-1)
1303 for(j=0;j<n1 ;j++)
1304 for(k=0;k<n2 ;k++)
1305 sum += points[i][j][k].C[3];
1306
1307 for(j=0;j<n1 ;j+=n1-1)
1308 for(i=0;i<n0 ;i++)
1309 for(k=0;k<n2 ;k++)
1310 sum += points[i][j][k].C[3];
1311
1312 for(k=0;k<n2 ;k+=n2-1)
1313 for(i=0;i<n0 ;i++)
1314 for(j=0;j<n1 ;j++)
1315 sum += points[i][j][k].C[3];
1316
1317 for(i=0;i<n0 ;i+=n0-1)
1318 for(j=0;j<n1 ;j+=n1-1)
1319 for(k=0;k<n2 ;k+=n2-1)
1320 sum -= 2*points[i][j][k].C[3];
1321
1322 sum *= gridCP->dv;
1323 tmp = g_strdup_printf(
1324 _("sum of values on the 6 faces of the cube(should be near to 0) = %lf \n"),
1325 sum
1326 );
1327 result = addToResult(result, tmp);
1328 if(tmp) g_free(tmp);
1329
1330 result = addToResult(result, "---------------------------------------------------------------------\n");
1331 if(gridCP->criticalPoints)
1332 {
1333 tmp = g_strdup_printf("%14s %14s %28s %10s %10s %s\n", " ",_("Position(Ang)")," ",_("Nearest at."),_(" AIM Charge "), _("Old charge (read from CCP output file)"));
1334 result = addToResult(result, tmp);
1335 if(tmp) g_free(tmp);
1336 }
1337 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1338 {
1339 CriticalPoint* data = (CriticalPoint*) list->data;
1340 gint i = data->index[0];
1341 gint j = data->index[1];
1342 gint k = data->index[2];
1343 gint c= data->numCenter;
1344
1345 tmp = g_strdup_printf("%+14.8f %+14.8f %+14.8f ",
1346 gridCP->grid->point[i][j][k].C[0]*BOHR_TO_ANG,
1347 gridCP->grid->point[i][j][k].C[1]*BOHR_TO_ANG,
1348 gridCP->grid->point[i][j][k].C[2]*BOHR_TO_ANG);
1349 result = addToResult(result, tmp);
1350 if(tmp) g_free(tmp);
1351
1352 tmp = g_strdup_printf("%2s[%d]%5s", GeomOrb[c].Prop.symbol, c+1," ");
1353 result = addToResult(result, tmp);
1354 if(tmp) g_free(tmp);
1355 tmp = g_strdup_printf("%+14.8f %+14.8f \n",
1356 data->nuclearCharge-data->integral,
1357 GeomOrb[c].partialCharge
1358 );
1359 result = addToResult(result, tmp);
1360 if(tmp) g_free(tmp);
1361 }
1362 result = addToResult(result, "---------------------------------------------------------------------\n");
1363
1364 nc = 0;
1365 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1366 {
1367 CriticalPoint* data = (CriticalPoint*) list->data;
1368 gint i = data->index[0];
1369 gint j = data->index[1];
1370 gint k = data->index[2];
1371 gint c= data->numCenter;
1372
1373 if(data->volume/gridCP->dv<=8) continue;
1374 nc++;
1375 tmp = g_strdup_printf(_("Attracteur number %d\n"),nc);
1376 result = addToResult(result, tmp);
1377 if(tmp) g_free(tmp);
1378 result = addToResult(result,"====================\n");
1379
1380 tmp = g_strdup_printf(_("Position(Ang) = %lf %lf %lf\n"),
1381 gridCP->grid->point[i][j][k].C[0]*BOHR_TO_ANG,
1382 gridCP->grid->point[i][j][k].C[1]*BOHR_TO_ANG,
1383 gridCP->grid->point[i][j][k].C[2]*BOHR_TO_ANG);
1384 result = addToResult(result, tmp);
1385 if(tmp) g_free(tmp);
1386
1387 tmp = g_strdup_printf("Index = %d %d %d rho = %14.10e\n",i,j,k,gridCP->grid->point[i][j][k].C[3]);
1388 result = addToResult(result, tmp);
1389 if(tmp) g_free(tmp);
1390
1391 tmp = g_strdup_printf(_("Nearest atom = %s[%d]\n"),
1392 GeomOrb[c].Prop.symbol,
1393 c+1);
1394 result = addToResult(result, tmp);
1395 if(tmp) g_free(tmp);
1396
1397 tmp = g_strdup_printf(_("Number of electrons in the volume of this attractor = %lf\n"), data->integral);
1398 result = addToResult(result, tmp);
1399 if(tmp) g_free(tmp);
1400 tmp = g_strdup_printf(_("Nuclear charge of the nearest atom = %lf\n"),data->nuclearCharge);
1401 result = addToResult(result, tmp);
1402 if(tmp) g_free(tmp);
1403 tmp = g_strdup_printf(_("Charge in the volume of this attractor = %lf\n"),
1404 data->nuclearCharge-data->integral);
1405 result = addToResult(result, tmp);
1406 if(tmp) g_free(tmp);
1407 tmp = g_strdup_printf(_("# of points in this volume = %d\n"), (gint)(data->volume/gridCP->dv));
1408 result = addToResult(result, tmp);
1409 if(tmp) g_free(tmp);
1410 }
1411 result = addToResult(result, "---------------------------------------------------------------------\n");
1412 tmp = g_strdup_printf(_("Total number of electrons = %lf\n"), gridCP->integral);
1413 result = addToResult(result, tmp);
1414 if(tmp) g_free(tmp);
1415 tmp = g_strdup_printf(_("Total nuclear charges = %lf\n"),gridCP->nuclearCharge);
1416 result = addToResult(result, tmp);
1417 if(tmp) g_free(tmp);
1418 tmp = g_strdup_printf(_("Total charge = %lf\n"), gridCP->nuclearCharge-gridCP->integral);
1419 result = addToResult(result, tmp);
1420 if(tmp) g_free(tmp);
1421 result = addToResult(result, "---------------------------------------------------------------------\n");
1422
1423 if(result && !CancelCalcul)
1424 {
1425 showResultDlg(result,_("AIM charges"),gridCP);
1426 }
1427 else if(!result && !CancelCalcul)
1428 {
1429 GtkWidget* message = MessageTxt(_("Oups a problem...."),_("Attractors"));
1430 gtk_window_set_modal (GTK_WINDOW (message), TRUE);
1431 gtk_window_set_transient_for(GTK_WINDOW(message),GTK_WINDOW(PrincipalWindow));
1432 }
1433 if(result) g_free(result);
1434 }
1435 /**************************************************************************/
computeNumCenters(GridCP * gridCP)1436 static void computeNumCenters(GridCP* gridCP)
1437 {
1438 GList* list=NULL;
1439 gint c;
1440
1441 gdouble dx, dy, dz;
1442 gdouble r;
1443 gdouble rold = 0;
1444
1445 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1446 {
1447 CriticalPoint* data = (CriticalPoint*) list->data;
1448 gint i = data->index[0];
1449 gint j = data->index[1];
1450 gint k = data->index[2];
1451 data->numCenter = 0;
1452 for(c=0; c<(gint)nCenters; c++)
1453 {
1454 dx = GeomOrb[c].C[0]-gridCP->grid->point[i][j][k].C[0];
1455 dy = GeomOrb[c].C[1]-gridCP->grid->point[i][j][k].C[1];
1456 dz = GeomOrb[c].C[2]-gridCP->grid->point[i][j][k].C[2];
1457 r = sqrt(dx*dx + dy*dy + dz*dz);
1458 if(c==0 || r<rold )
1459 {
1460 data->numCenter = c;
1461 rold = r;
1462 }
1463 }
1464 }
1465 }
1466 /**************************************************************************/
computeCharges(GridCP * gridCP)1467 void computeCharges(GridCP* gridCP)
1468 {
1469 GList* criticalPoint = gridCP->criticalPoints;
1470 gint*** vP = gridCP->volumeNumberOfPoints;
1471 Point5 ***points = gridCP->grid->point;
1472 GList* list = NULL;
1473 gint i,j,k;
1474 gint nc = 0;
1475 gint nc2 = 0;
1476 gdouble scal=0;
1477 gchar* str = _("Computing of charges..., Please wait");
1478 gdouble* integ = NULL;
1479 gdouble* volume = NULL;
1480
1481 computeNumCenters(gridCP);
1482 nc = 0;
1483 for(i=0;i<grid->N[0];i++)
1484 for(j=0;j<grid->N[1];j++)
1485 for(k=0;k<grid->N[2];k++)
1486 {
1487 gint n = abs(vP[i][j][k]);
1488 if(nc<n) nc = n;
1489 }
1490 for(list=criticalPoint;list!=NULL;list=list->next) nc2++;
1491
1492 if(nc>0) integ = g_malloc(nc*sizeof(gdouble));
1493 for(i=0;i<nc;i++) integ[i] = 0;
1494 if(nc>0) volume = g_malloc(nc*sizeof(gdouble));
1495 for(i=0;i<nc;i++) volume[i] = 0;
1496
1497 scal = 1.1/grid->N[0];
1498 progress_orb_txt(0,str,TRUE);
1499 if(nc>0)
1500 for(i=0;i<grid->N[0];i++)
1501 {
1502 progress_orb_txt(scal,str,FALSE);
1503 for(j=0;j<grid->N[1];j++)
1504 for(k=0;k<grid->N[2];k++)
1505 {
1506 gint n = abs(vP[i][j][k])-1;
1507
1508 if(n>=0) integ[n] += points[i][j][k].C[3];
1509 if(n>=0) volume[n] += 1;
1510 }
1511 }
1512 progress_orb_txt(0," ",TRUE);
1513
1514 for(list=criticalPoint;list!=NULL;list=list->next)
1515 {
1516 CriticalPoint* data=(CriticalPoint*)list->data;
1517 gint numV = data->numVolume;
1518 gint n = abs(numV)-1;
1519 data->integral = 0;
1520 data->volume = 0;
1521
1522 if(n>=0)
1523 {
1524 data->integral = integ[n];
1525 data->volume = volume[n];
1526 }
1527 data->integral *= gridCP->dv;
1528 data->volume *= gridCP->dv;
1529 data->nuclearCharge = GeomOrb[data->numCenter].nuclearCharge;
1530 }
1531 if(integ) g_free(integ);
1532 if(volume) g_free(volume);
1533 for(i=0;i<grid->N[0];i++)
1534 for(j=0;j<grid->N[1];j++)
1535 for(k=0;k<grid->N[2];k++)
1536 gridCP->integral += points[i][j][k].C[3];
1537 gridCP->integral *= gridCP->dv;
1538
1539 for(j=0; j<(gint)nCenters; j++)
1540 gridCP->nuclearCharge += GeomOrb[j].nuclearCharge;
1541 progress_orb_txt(0," ",TRUE);
1542 }
1543 /**************************************************************************/
computeAIMCharges(Grid * grid,gboolean ongrid)1544 void computeAIMCharges(Grid* grid, gboolean ongrid)
1545 {
1546 GridCP* gridCP = g_malloc(sizeof(GridCP));
1547
1548 if(!test_grid_all_positive(grid))
1549 {
1550 Message(_("Sorry\n The current grid is not a grid for electronic density!!!"),_("Error"),TRUE);
1551 return;
1552 }
1553
1554 initGridCP(gridCP, grid,NULL);
1555 assignGridCP(gridCP, ongrid);
1556 if(!CancelCalcul && !ongrid)
1557 {
1558 gint N = gridCP->grid->N[0]*gridCP->grid->N[1]*gridCP->grid->N[2]/1000;
1559 gint nold = refineEdge(gridCP, ongrid);
1560 gint n = refineEdge(gridCP, ongrid);
1561 gint iter = 0;
1562 while(abs(n-nold)>N)
1563 {
1564 iter++;
1565 if(iter>30) break;
1566 nold = n;
1567 n = refineEdge(gridCP, ongrid);
1568 }
1569 }
1570 if(!CancelCalcul) removeAttractor0(gridCP);
1571 if(!CancelCalcul) computeCharges(gridCP);
1572 if(!CancelCalcul) removeNonSignificantAttractor(gridCP);
1573 if(!CancelCalcul) showGridCP(gridCP);
1574 CancelCalcul = FALSE;
1575 progress_orb_txt(0," ",TRUE);
1576 }
1577 /**************************************************************************/
computeELFEletrons(GridCP * gridCP)1578 void computeELFEletrons(GridCP* gridCP)
1579 {
1580 GList* criticalPoint = gridCP->criticalPoints;
1581 gint*** vP = gridCP->volumeNumberOfPoints;
1582 Point5 ***points = gridCP->gridAux->point;
1583 GList* list = NULL;
1584 gint i,j,k;
1585
1586 computeNumCenters(gridCP);
1587
1588 progress_orb_txt(0,_("Computing of the number of electrons at each attractor..., Please wait"),TRUE);
1589 for(list=criticalPoint;list!=NULL;list=list->next)
1590 {
1591 CriticalPoint* data=(CriticalPoint*)list->data;
1592 gint numV = data->numVolume;
1593 data->integral = 0;
1594 data->volume = 0;
1595 for(i=0;i<grid->N[0];i++)
1596 for(j=0;j<grid->N[1];j++)
1597 for(k=0;k<grid->N[2];k++)
1598 {
1599 if(vP[i][j][k]==numV || vP[i][j][k]==-numV)
1600 {
1601 data->integral += points[i][j][k].C[3];
1602 data->volume += 1;
1603 }
1604 }
1605 data->integral *= gridCP->dv;
1606 data->volume *= gridCP->dv;
1607 data->nuclearCharge = GeomOrb[data->numCenter].nuclearCharge;
1608 }
1609 for(i=0;i<grid->N[0];i++)
1610 for(j=0;j<grid->N[1];j++)
1611 for(k=0;k<grid->N[2];k++)
1612 gridCP->integral += points[i][j][k].C[3];
1613 gridCP->integral *= gridCP->dv;
1614
1615 for(j=0; j<(gint)nCenters; j++)
1616 gridCP->nuclearCharge += GeomOrb[j].nuclearCharge;
1617 progress_orb_txt(0," ",TRUE);
1618 }
1619 /**************************************************************************/
showELFGridCP(GridCP * gridCP)1620 static void showELFGridCP(GridCP* gridCP)
1621 {
1622 GList* list=NULL;
1623 gchar* result =NULL;
1624 gchar* tmp = NULL;
1625 gint nc = 0;
1626 gdouble xx, yy, zz;
1627 gint n0 = 0, n1 = 0, n2 = 0;
1628 gint i,j,k;
1629 gdouble sum = 0;
1630 Point5*** points = gridCP->gridAux->point;
1631
1632 result = addToResult(result, _("Geometry (Ang)\n"));
1633 result = addToResult(result, "==============\n");
1634 for(j=0; j<(gint)nCenters; j++)
1635 {
1636 tmp = g_strdup_printf("%s[%d] %lf %lf %lf\n",
1637 GeomOrb[j].Prop.symbol,
1638 j+1,
1639 GeomOrb[j].C[0]*BOHR_TO_ANG,
1640 GeomOrb[j].C[1]*BOHR_TO_ANG,
1641 GeomOrb[j].C[2]*BOHR_TO_ANG);
1642 result = addToResult(result, tmp);
1643 if(tmp) g_free(tmp);
1644 }
1645 result = addToResult(result, "---------------------------------------------------------------------\n");
1646
1647 n0 = gridCP->grid->N[0];
1648 n1 = gridCP->grid->N[1];
1649 n2 = gridCP->grid->N[2];
1650
1651 xx = points[n0-1][0][0].C[0]-points[0][0][0].C[0];
1652 yy = points[0][n1-1][0].C[1]-points[0][0][0].C[1];
1653 zz = points[0][0][n2-1].C[2]-points[0][0][0].C[2];
1654 tmp = g_strdup_printf(
1655 _(
1656 "Grid point density (Ang^-1) on the first direction(>10 is recommended) = %lf\n"
1657 "density of the grid(Ang^-1) on the second direction(>10 is recommended) = %lf\n"
1658 "density of the grid(Ang^-1) on the third direction(>10 is recommended) = %lf\n"
1659 ),
1660 n0/(xx*BOHR_TO_ANG),
1661 n1/(yy*BOHR_TO_ANG),
1662 n2/(zz*BOHR_TO_ANG)
1663 );
1664 result = addToResult(result, tmp);
1665 if(tmp) g_free(tmp);
1666 sum = 0;
1667 for(i=0;i<n0 ;i+=n0-1)
1668 for(j=0;j<n1 ;j++)
1669 for(k=0;k<n2 ;k++)
1670 sum += points[i][j][k].C[3];
1671
1672 for(j=0;j<n1 ;j+=n1-1)
1673 for(i=0;i<n0 ;i++)
1674 for(k=0;k<n2 ;k++)
1675 sum += points[i][j][k].C[3];
1676
1677 for(k=0;k<n2 ;k+=n2-1)
1678 for(i=0;i<n0 ;i++)
1679 for(j=0;j<n1 ;j++)
1680 sum += points[i][j][k].C[3];
1681
1682 for(i=0;i<n0 ;i+=n0-1)
1683 for(j=0;j<n1 ;j+=n1-1)
1684 for(k=0;k<n2 ;k+=n2-1)
1685 sum -= 2*points[i][j][k].C[3];
1686
1687 sum *= gridCP->dv;
1688 tmp = g_strdup_printf(
1689 _("sum of values on the 6 faces of the cube(should be near to 0) = %lf \n"),
1690 sum
1691 );
1692 result = addToResult(result, tmp);
1693 if(tmp) g_free(tmp);
1694
1695 result = addToResult(result, "---------------------------------------------------------------------\n");
1696 if(gridCP->criticalPoints)
1697 {
1698 tmp = g_strdup_printf("%14s %14s %28s %10s %10s\n", " ",_("Position(Ang)")," ",_("Nearest at."),_(" # electrons"));
1699 result = addToResult(result, tmp);
1700 if(tmp) g_free(tmp);
1701 }
1702 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1703 {
1704 CriticalPoint* data = (CriticalPoint*) list->data;
1705 gint i = data->index[0];
1706 gint j = data->index[1];
1707 gint k = data->index[2];
1708 gint c= data->numCenter;
1709
1710 tmp = g_strdup_printf("%+14.8f %+14.8f %+14.8f ",
1711 gridCP->grid->point[i][j][k].C[0]*BOHR_TO_ANG,
1712 gridCP->grid->point[i][j][k].C[1]*BOHR_TO_ANG,
1713 gridCP->grid->point[i][j][k].C[2]*BOHR_TO_ANG);
1714 result = addToResult(result, tmp);
1715 if(tmp) g_free(tmp);
1716
1717 tmp = g_strdup_printf("%2s[%d]%5s", GeomOrb[c].Prop.symbol, c+1," ");
1718 result = addToResult(result, tmp);
1719 if(tmp) g_free(tmp);
1720 tmp = g_strdup_printf("%+14.8f\n", data->integral);
1721 result = addToResult(result, tmp);
1722 if(tmp) g_free(tmp);
1723 }
1724 result = addToResult(result, "---------------------------------------------------------------------\n");
1725
1726 nc = 0;
1727 for(list=gridCP->criticalPoints;list!=NULL;list=list->next)
1728 {
1729 CriticalPoint* data = (CriticalPoint*) list->data;
1730 gint i = data->index[0];
1731 gint j = data->index[1];
1732 gint k = data->index[2];
1733 gint c= data->numCenter;
1734
1735 if(data->volume/gridCP->dv<=8) continue;
1736 nc++;
1737 tmp = g_strdup_printf(_("Attracteur number %d\n"),nc);
1738 result = addToResult(result, tmp);
1739 if(tmp) g_free(tmp);
1740 result = addToResult(result,"====================\n");
1741
1742 tmp = g_strdup_printf(_("Position(Ang) = %lf %lf %lf\n"),
1743 gridCP->grid->point[i][j][k].C[0]*BOHR_TO_ANG,
1744 gridCP->grid->point[i][j][k].C[1]*BOHR_TO_ANG,
1745 gridCP->grid->point[i][j][k].C[2]*BOHR_TO_ANG);
1746 result = addToResult(result, tmp);
1747 if(tmp) g_free(tmp);
1748
1749 tmp = g_strdup_printf("Index = %d %d %d rho = %14.10e\n",i,j,k,gridCP->gridAux->point[i][j][k].C[3]);
1750 result = addToResult(result, tmp);
1751 if(tmp) g_free(tmp);
1752
1753 tmp = g_strdup_printf(_("Nearest atom = %s[%d]\n"),
1754 GeomOrb[c].Prop.symbol,
1755 c+1);
1756 result = addToResult(result, tmp);
1757 if(tmp) g_free(tmp);
1758
1759 tmp = g_strdup_printf(_("Number of electrons in the volume of this attractor = %lf\n"), data->integral);
1760 result = addToResult(result, tmp);
1761 if(tmp) g_free(tmp);
1762 tmp = g_strdup_printf(_("Nuclear charge of the nearest atom = %lf\n"),data->nuclearCharge);
1763 result = addToResult(result, tmp);
1764 if(tmp) g_free(tmp);
1765 tmp = g_strdup_printf(_("# of points in this volume = %d\n"), (gint)(data->volume/gridCP->dv));
1766 result = addToResult(result, tmp);
1767 if(tmp) g_free(tmp);
1768 }
1769 result = addToResult(result, "---------------------------------------------------------------------\n");
1770 tmp = g_strdup_printf(_("Total number of electrons = %lf\n"), gridCP->integral);
1771 result = addToResult(result, tmp);
1772 result = addToResult(result, "---------------------------------------------------------------------\n");
1773
1774 if(result && !CancelCalcul)
1775 {
1776 GtkWidget* message = MessageTxt(result,_("ELF analysis"));
1777 gtk_window_set_modal (GTK_WINDOW (message), TRUE);
1778 gtk_window_set_transient_for(GTK_WINDOW(message),GTK_WINDOW(PrincipalWindow));
1779 }
1780 else if(!result && !CancelCalcul)
1781 {
1782 GtkWidget* message = MessageTxt(_("Oups a problem...."),_("Attractors"));
1783 gtk_window_set_modal (GTK_WINDOW (message), FALSE);
1784 gtk_window_set_transient_for(GTK_WINDOW(message),GTK_WINDOW(PrincipalWindow));
1785 }
1786 if(result) g_free(result);
1787 }
1788 /**************************************************************************/
computeELFAttractors(Grid * gridELF,Grid * gridDens,gboolean ongrid)1789 void computeELFAttractors(Grid* gridELF, Grid* gridDens, gboolean ongrid)
1790 {
1791 GridCP* gridCP = g_malloc(sizeof(GridCP));
1792
1793 if(!test_grid_all_positive(gridELF))
1794 {
1795 Message(_("Sorry\n The current grid is not a grid for ELF!!!"),_("Error"),TRUE);
1796 return;
1797 }
1798 if(!test_grid_all_positive(gridDens))
1799 {
1800 Message(_("Sorry\n The second grid is not a grid for electronic density!!!"),_("Error"),TRUE);
1801 return;
1802 }
1803 if(gridELF && gridDens)
1804 {
1805 gint c;
1806 for(c=0;c<3;c++)
1807 if(gridELF->N[c] != gridDens->N[c])
1808 {
1809 Message(_("Sorry\n The Cubes of the 2 grids should be equals!!!"),_("Error"),TRUE);
1810 return;
1811 }
1812 }
1813
1814 initGridCP(gridCP, gridELF,gridDens);
1815 assignGridCP(gridCP, ongrid);
1816 if(!CancelCalcul) refineEdge(gridCP,ongrid);
1817 if(!CancelCalcul) removeAttractor0(gridCP);
1818 if(!CancelCalcul) assignPointsZero(gridCP);
1819 if(!CancelCalcul) computeELFEletrons(gridCP);
1820 if(!CancelCalcul) removeNonSignificantAttractor(gridCP);
1821 if(!CancelCalcul) showELFGridCP(gridCP);
1822 CancelCalcul = FALSE;
1823 progress_orb_txt(0," ",TRUE);
1824 }
1825 #undef TOL
1826