1 /* MoleculeSE.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 #include "../../Config.h"
20 #include <stdlib.h>
21 #include <math.h>
22 #include <ctype.h>
23
24 #include "../Common/Global.h"
25 #include "../Utils/AtomsProp.h"
26 #include "../Geometry/Fragments.h"
27 #include "../Geometry/DrawGeom.h"
28 #include "../Geometry/Measure.h"
29 #include "../Geometry/GeomGlobal.h"
30 #include "../Geometry/GeomXYZ.h"
31 #include "../Utils/Constants.h"
32 #include "../Utils/Utils.h"
33 #include "AtomSE.h"
34 #include "MoleculeSE.h"
35 void create_GeomXYZ_from_draw_grometry();
36
37 static gboolean** bondedMatrix = NULL;
38
39 #define BOHR_TO_ANG 0.52917726
40
41 /**********************************************************************/
newMoleculeSE()42 MoleculeSE newMoleculeSE()
43 {
44 gint i;
45 MoleculeSE molecule;
46
47 molecule.nAtoms = 0;
48 molecule.totalCharge = 0;
49 molecule.spinMultiplicity = 0;
50 molecule.atoms = NULL;
51 molecule.energy = 0;
52
53 for(i=0;i<3;i++)
54 molecule.gradient[i] = NULL;
55 molecule.numberOf2Connections = 0;
56 for(i=0;i<2;i++)
57 molecule.connected2[i] = NULL;
58 molecule.numberOf3Connections = 0;
59 for(i=0;i<3;i++)
60 molecule.connected3[i] = NULL;
61
62 return molecule;
63
64 }
65 /**********************************************************************/
freeMoleculeSE(MoleculeSE * molecule)66 void freeMoleculeSE(MoleculeSE* molecule)
67 {
68
69 gint i;
70 if(molecule->nAtoms<=0)
71 return;
72
73 if(molecule->atoms != NULL)
74 {
75 for(i=0;i<molecule->nAtoms;i++)
76 {
77 if(molecule->atoms[i].prop.symbol != NULL)
78 g_free(molecule->atoms[i].prop.symbol);
79 if(molecule->atoms[i].mmType !=NULL )
80 g_free(molecule->atoms[i].mmType);
81 if(molecule->atoms[i].pdbType !=NULL )
82 g_free(molecule->atoms[i].pdbType);
83 }
84
85 g_free(molecule->atoms);
86 molecule->atoms = NULL;
87 }
88 molecule->nAtoms = 0;
89 molecule->energy = 0;
90 molecule->totalCharge = 0;
91 molecule->spinMultiplicity = 0;
92
93 for(i=0;i<3;i++)
94 if(molecule->gradient[i] != NULL)
95 {
96 g_free(molecule->gradient[i]);
97 molecule->gradient[i] = NULL;
98 }
99 molecule->numberOf2Connections = 0;
100 for(i=0;i<2;i++)
101 {
102 if(molecule->connected2[i] != NULL)
103 g_free(molecule->connected2[i]);
104 molecule->connected2[i] = NULL;
105 }
106 molecule->numberOf3Connections = 0;
107 for(i=0;i<3;i++)
108 {
109 if(molecule->connected3[i] != NULL)
110 g_free(molecule->connected3[i]);
111 molecule->connected3[i] = NULL;
112 }
113 }
114 /*****************************************************************************/
createBondedMatrix(MoleculeSE * molecule)115 static void createBondedMatrix(MoleculeSE* molecule)
116 {
117 gint nAtoms = molecule->nAtoms;
118 gint i;
119 gint j;
120
121 if(nAtoms<1)
122 return;
123
124 bondedMatrix = g_malloc(nAtoms*sizeof(gboolean*));
125 for(i=0;i<nAtoms;i++)
126 bondedMatrix[i] = g_malloc(nAtoms*sizeof(gboolean));
127
128 for(i=0;i<nAtoms;i++)
129 {
130 for(j=0;j<nAtoms;j++)
131 bondedMatrix[i][j] = FALSE;
132
133 bondedMatrix[i][i] = TRUE;
134 }
135
136 }
137 /*****************************************************************************/
freeBondedMatrix(MoleculeSE * molecule)138 static void freeBondedMatrix(MoleculeSE* molecule)
139 {
140 gint nAtoms = molecule->nAtoms;
141 gint i;
142
143 if(bondedMatrix == NULL)
144 return;
145 for(i=0;i<nAtoms;i++)
146 if(bondedMatrix[i] != NULL)
147 g_free(bondedMatrix[i]);
148
149 g_free(bondedMatrix);
150 bondedMatrix = NULL;
151
152 }
153 /*****************************************************************************/
updatebondedMatrix(gint a1,gint a2)154 static void updatebondedMatrix(gint a1, gint a2)
155 {
156 bondedMatrix[a1][a2] = TRUE;
157 bondedMatrix[a2][a1] = TRUE;
158
159 }
160 /*****************************************************************************/
isConnected2(MoleculeSE * molecule,gint i,gint j)161 static gboolean isConnected2(MoleculeSE* molecule,gint i,gint j)
162 {
163 gdouble distance;
164 gdouble dij;
165 gint k;
166 AtomSE a1 = molecule->atoms[i];
167 AtomSE a2 = molecule->atoms[j];
168
169 if(molecule->atoms[i].typeConnections)
170 {
171 gint nj = molecule->atoms[j].N-1;
172 if(molecule->atoms[i].typeConnections[nj]>0) return TRUE;
173 else return FALSE;
174 }
175 distance = 0;
176 for (k=0;k<3;k++)
177 {
178 dij = a1.coordinates[k]-a2.coordinates[k];
179 distance +=dij*dij;
180 }
181 distance = sqrt(distance)/BOHR_TO_ANG;
182
183 if(distance<(a1.prop.covalentRadii+a2.prop.covalentRadii)) return TRUE;
184 else return FALSE;
185 }
186 /*****************************************************************************/
set2Connections(MoleculeSE * molecule)187 static void set2Connections(MoleculeSE* molecule)
188 {
189 gint i;
190 gint j;
191 gint k=0;
192
193 k = molecule->nAtoms;
194 k = k*(k-1)/2;
195 for(i=0;i<2;i++)
196 molecule->connected2[i] = g_malloc(k*sizeof(gint));
197
198 k=0;
199 for(i=0;i<molecule->nAtoms-1;i++)
200 for(j=i+1;j<molecule->nAtoms;j++)
201 {
202 if(isConnected2(molecule,i,j))
203 {
204 molecule->connected2[0][k]= i;
205 molecule->connected2[1][k]= j;
206
207 updatebondedMatrix(i,j);
208
209 k++;
210
211 }
212 }
213 molecule->numberOf2Connections = k;
214 if(k==0)
215 for(i=0;i<2;i++)
216 {
217 g_free(molecule->connected2[i]);
218 molecule->connected2[i] = NULL;
219 }
220 else
221 for(i=0;i<2;i++)
222 molecule->connected2[i] = g_realloc(molecule->connected2[i],k*sizeof(gint));
223 /* printing for test*/
224 /*
225 printf("%d 2 connections : \n",molecule->numberOf2Connections);
226 for(k=0;k<molecule->numberOf2Connections;k++)
227 {
228
229 i = molecule->connected2[0][k];
230 j = molecule->connected2[1][k];
231 printf("%d-%d ",i,j);
232 }
233 printf("\n");
234 */
235
236
237 }
238 /*****************************************************************************/
permut(gint * a,gint * b)239 static void permut(gint* a,gint *b)
240 {
241 gint c = *a;
242 *a = *b;
243 *b = c;
244 }
245 /*****************************************************************************/
isConnected3(MoleculeSE * molecule,gint n,gint i,gint j,gint k)246 static gboolean isConnected3(MoleculeSE* molecule,gint n,gint i,gint j, gint k)
247 {
248 gint c;
249 gint a1,a2,a3;
250 for(c=0;c<n;c++)
251 {
252 a1 = molecule->connected3[0][c];
253 a2 = molecule->connected3[1][c];
254 a3 = molecule->connected3[2][c];
255 if(a1==i && a2 == j && a3 == k)
256 return TRUE;
257 }
258 return FALSE;
259
260 }
261 /*****************************************************************************/
connect3(MoleculeSE * molecule,gint n,gint i,gint j,gint k)262 static gboolean connect3(MoleculeSE* molecule,gint n,gint i,gint j, gint k)
263 {
264 if(i>k)permut(&i,&k);
265 if(!isConnected3(molecule,n,i,j,k))
266 {
267 molecule->connected3[0][n]= i;
268 molecule->connected3[1][n]= j;
269 molecule->connected3[2][n]= k;
270
271 updatebondedMatrix(i,j);
272 updatebondedMatrix(i,k);
273 updatebondedMatrix(j,k);
274
275 return TRUE;
276 }
277 return FALSE;
278
279 }
280 /*****************************************************************************/
set3Connections(MoleculeSE * molecule)281 static void set3Connections(MoleculeSE* molecule)
282 {
283 gint i;
284 gint j;
285 gint k=0;
286 gint l=0;
287 gint n=0;
288
289 k = molecule->numberOf2Connections*molecule->nAtoms;
290 for(i=0;i<3;i++)
291 molecule->connected3[i] = g_malloc(k*sizeof(gint));
292
293 n=0;
294 for(k=0;k<molecule->numberOf2Connections;k++)
295 {
296 i = molecule->connected2[0][k];
297 j = molecule->connected2[1][k];
298 for(l=0;l<molecule->nAtoms;l++)
299 {
300 if(l!=i && l!=j)
301 {
302 if( isConnected2(molecule,i,l))
303 if( connect3(molecule,n,l,i,j))
304 n++;
305
306 if( isConnected2(molecule,j,l))
307 if( connect3(molecule,n,i,j,l))
308 n++;
309 }
310 }
311
312 }
313 molecule->numberOf3Connections = n;
314 if(n==0)
315 for(i=0;i<3;i++)
316 {
317 g_free(molecule->connected3[i]);
318 molecule->connected3[i] = NULL;
319 }
320 else
321 for(i=0;i<3;i++)
322 molecule->connected3[i] = g_realloc(molecule->connected3[i],n*sizeof(gint));
323 /* printing for test*/
324 /*
325 printf("%d 3 connections : \n",molecule->numberOf3Connections);
326 for(k=0;k<molecule->numberOf3Connections;k++)
327 {
328
329 i = molecule->connected3[0][k];
330 j = molecule->connected3[1][k];
331 l = molecule->connected3[2][k];
332 printf("%d-%d-%d ",i,j,l);
333 }
334 printf("\n");
335 */
336
337
338 }
339 /*****************************************************************************/
setConnectionsMoleculeSE(MoleculeSE * molecule)340 void setConnectionsMoleculeSE(MoleculeSE* molecule)
341 {
342 createBondedMatrix(molecule);
343
344 /* printf("Set Connection\n");*/
345 set_text_to_draw(_("Establishing connectivity : 2 connections..."));
346 set_statubar_operation_str(_("Establishing connectivity : 2 connections..."));
347 drawGeom();
348 while( gtk_events_pending() )
349 gtk_main_iteration();
350 set2Connections(molecule);
351 set_text_to_draw(_("Establishing connectivity : 3 connections..."));
352 set_statubar_operation_str(_("Establishing connectivity : 3 connections..."));
353
354 drawGeom();
355 if(StopCalcul) return;
356 while( gtk_events_pending() ) gtk_main_iteration();
357 set3Connections(molecule);
358
359 freeBondedMatrix(molecule);
360 }
361 /*****************************************************************************/
createMoleculeSE(GeomDef * geom,gint natoms,gint charge,gint spin,gboolean connections)362 MoleculeSE createMoleculeSE(GeomDef* geom,gint natoms, gint charge, gint spin, gboolean connections)
363 {
364
365 gint i;
366 MoleculeSE molecule = newMoleculeSE();
367
368 molecule.nAtoms = natoms;
369 molecule.atoms = g_malloc(molecule.nAtoms*sizeof(AtomSE));
370 for(i=0;i<molecule.nAtoms;i++)
371 {
372 molecule.atoms[i].prop = prop_atom_get(geom[i].Prop.symbol);
373 molecule.atoms[i].coordinates[0] = geom[i].X*BOHR_TO_ANG;
374 molecule.atoms[i].coordinates[1] = geom[i].Y*BOHR_TO_ANG;
375 molecule.atoms[i].coordinates[2] = geom[i].Z*BOHR_TO_ANG;
376 molecule.atoms[i].charge = geom[i].Charge;
377 molecule.atoms[i].mmType = g_strdup(geom[i].mmType);
378 molecule.atoms[i].pdbType = g_strdup(geom[i].pdbType);
379 molecule.atoms[i].residueName = g_strdup(geom[i].Residue);
380 molecule.atoms[i].residueNumber = geom[i].ResidueNumber;
381 molecule.atoms[i].N = geom[i].N;
382 molecule.atoms[i].layer = geom[i].Layer;
383 molecule.atoms[i].show = geom[i].show;
384 molecule.atoms[i].variable = geom[i].Variable;
385 molecule.atoms[i].typeConnections = NULL;
386 if(geom[i].typeConnections)
387 {
388 gint j;
389 molecule.atoms[i].typeConnections = g_malloc(molecule.nAtoms*sizeof(gint));
390 for(j=0;j<molecule.nAtoms;j++)
391 {
392 gint nj = geom[j].N-1;
393 molecule.atoms[i].typeConnections[nj] = geom[i].typeConnections[nj];
394 }
395 }
396 }
397 if(connections) setConnectionsMoleculeSE(&molecule);
398 molecule.totalCharge = charge;
399 molecule.spinMultiplicity = spin;
400
401 for(i=0;i<3;i++) /* x, y and z derivatives */
402 molecule.gradient[i] = g_malloc(molecule.nAtoms*sizeof(gdouble));
403 for(i=0;i<3;i++) molecule.dipole[i] = 0.0;
404
405 return molecule;
406 }
407 /*****************************************************************************/
createFromGeomXYZMoleculeSE(gint charge,gint spin,gboolean connections)408 MoleculeSE createFromGeomXYZMoleculeSE(gint charge, gint spin, gboolean connections)
409 {
410
411 gdouble x,y,z;
412 gint i;
413 MoleculeSE molecule = newMoleculeSE();
414
415 molecule.nAtoms = NcentersXYZ;
416 molecule.atoms = g_malloc(molecule.nAtoms*sizeof(AtomSE));
417
418 for(i=0;i<molecule.nAtoms;i++)
419 {
420 molecule.atoms[i].prop = prop_atom_get(GeomXYZ[i].Symb);
421 if(!test(GeomXYZ[i].X)) x = get_value_variableXYZ(GeomXYZ[i].X);
422 else x = atof(GeomXYZ[i].X);
423 if(!test(GeomXYZ[i].Y)) y = get_value_variableXYZ(GeomXYZ[i].Y);
424 else y = atof(GeomXYZ[i].Y);
425 if(!test(GeomXYZ[i].Z)) z = get_value_variableXYZ(GeomXYZ[i].Z);
426 else z = atof(GeomXYZ[i].Z);
427
428 if(Units==0)
429 {
430 x *= BOHR_TO_ANG;
431 y *= BOHR_TO_ANG;
432 z *= BOHR_TO_ANG;
433 }
434 molecule.atoms[i].coordinates[0] = x;
435 molecule.atoms[i].coordinates[1] = y;
436 molecule.atoms[i].coordinates[2] = z;
437 molecule.atoms[i].charge = atof(GeomXYZ[i].Charge);
438 molecule.atoms[i].mmType = g_strdup(GeomXYZ[i].mmType);
439 molecule.atoms[i].pdbType = g_strdup(GeomXYZ[i].pdbType);
440 molecule.atoms[i].residueName = g_strdup(GeomXYZ[i].Residue);
441 molecule.atoms[i].residueNumber = GeomXYZ[i].ResidueNumber;
442 molecule.atoms[i].N = i+1;
443 if(strstr(GeomXYZ[i].Layer,"Low")) molecule.atoms[i].layer = LOW_LAYER;
444 if(strstr(GeomXYZ[i].Layer,"Medium")) molecule.atoms[i].layer = MEDIUM_LAYER;
445 else molecule.atoms[i].layer = HIGH_LAYER;
446 molecule.atoms[i].show = TRUE;
447 molecule.atoms[i].variable = FALSE;
448 if(
449 test(GeomXYZ[i].X) ||
450 test(GeomXYZ[i].Y) ||
451 test(GeomXYZ[i].Z)
452 )
453 molecule.atoms[i].variable = TRUE;
454 molecule.atoms[i].typeConnections = NULL;
455 if(GeomXYZ[i].typeConnections)
456 {
457 gint j;
458 molecule.atoms[i].typeConnections = g_malloc(molecule.nAtoms*sizeof(gint));
459 for(j=0;j<molecule.nAtoms;j++)
460 {
461 gint nj = j;
462 molecule.atoms[i].typeConnections[nj] = GeomXYZ[i].typeConnections[nj];
463 }
464 }
465 }
466 if(connections) setConnectionsMoleculeSE(&molecule);
467 molecule.totalCharge = charge;
468 molecule.spinMultiplicity = spin;
469
470 for(i=0;i<3;i++) /* x, y and z derivatives */
471 molecule.gradient[i] = g_malloc(molecule.nAtoms*sizeof(gdouble));
472 for(i=0;i<3;i++) molecule.dipole[i] = 0.0;
473
474 return molecule;
475 }
476 /*****************************************************************************/
redrawMoleculeSE(MoleculeSE * molecule,gchar * str)477 void redrawMoleculeSE(MoleculeSE* molecule,gchar* str)
478 {
479 gint i;
480 gdouble C[3] = {0.0,0.0,0.0};
481
482 Free_One_Geom(geometry0,Natoms);
483 Free_One_Geom(geometry ,Natoms);
484 Natoms = 0;
485 geometry0 = NULL;
486 geometry = NULL;
487
488 Natoms = molecule->nAtoms;
489 geometry0 = g_malloc((Natoms)*sizeof(GeomDef));
490 geometry = g_malloc((Natoms)*sizeof(GeomDef));
491
492 for(i=0;i<(gint)Natoms;i++)
493 {
494 geometry0[i].X = molecule->atoms[i].coordinates[0];
495 geometry0[i].Y = molecule->atoms[i].coordinates[1];
496 geometry0[i].Z = molecule->atoms[i].coordinates[2];
497 geometry0[i].Charge = molecule->atoms[i].charge;
498 geometry0[i].Prop = prop_atom_get(molecule->atoms[i].prop.symbol);
499 geometry0[i].pdbType = g_strdup(molecule->atoms[i].pdbType);
500 geometry0[i].mmType = g_strdup(molecule->atoms[i].mmType);
501 geometry0[i].Residue = g_strdup(molecule->atoms[i].residueName);
502 geometry0[i].ResidueNumber = molecule->atoms[i].residueNumber;
503 geometry0[i].show = molecule->atoms[i].show;
504 geometry0[i].Variable = molecule->atoms[i].variable;
505 geometry0[i].Layer = molecule->atoms[i].layer;
506 geometry0[i].N = molecule->atoms[i].N;
507
508 geometry[i].X = molecule->atoms[i].coordinates[0];
509 geometry[i].Y = molecule->atoms[i].coordinates[1];
510 geometry[i].Z = molecule->atoms[i].coordinates[2];
511 geometry[i].Charge = molecule->atoms[i].charge;
512 geometry[i].Prop = prop_atom_get(molecule->atoms[i].prop.symbol);
513 geometry[i].pdbType = g_strdup(molecule->atoms[i].pdbType);
514 geometry[i].mmType = g_strdup(molecule->atoms[i].mmType);
515 geometry[i].Residue = g_strdup(molecule->atoms[i].residueName);
516 geometry[i].ResidueNumber = molecule->atoms[i].residueNumber;
517 geometry[i].show = molecule->atoms[i].show;
518 geometry[i].Variable = molecule->atoms[i].variable;
519 geometry[i].Layer = molecule->atoms[i].layer;
520 geometry[i].N = molecule->atoms[i].N;
521
522 C[0] += geometry0[i].X;
523 C[1] += geometry0[i].Y;
524 C[2] += geometry0[i].Z;
525
526
527 }
528 for(i=0;i<3;i++)
529 C[i] /= Natoms;
530 /* center */
531 for(i=0;i<(gint)Natoms;i++)
532 {
533 geometry0[i].X -= C[0];
534 geometry0[i].Y -= C[1];
535 geometry0[i].Z -= C[2];
536
537 geometry[i].X -= C[0];
538 geometry[i].Y -= C[1];
539 geometry[i].Z -= C[2];
540 }
541
542 for(i=0;i<(gint)Natoms;i++)
543 {
544 geometry0[i].X /=BOHR_TO_ANG;
545 geometry0[i].Y /=BOHR_TO_ANG;
546 geometry0[i].Z /=BOHR_TO_ANG;
547
548 geometry[i].X /=BOHR_TO_ANG;
549 geometry[i].Y /=BOHR_TO_ANG;
550 geometry[i].Z /=BOHR_TO_ANG;
551 }
552 unselect_all_atoms();
553 set_text_to_draw(str);
554 set_statubar_operation_str(str);
555 change_of_center(NULL,NULL);
556 reset_all_connections();
557 create_GeomXYZ_from_draw_grometry();
558 RebuildGeom = TRUE;
559 drawGeom();
560
561 while( gtk_events_pending() )
562 gtk_main_iteration();
563
564 }
565 /********************************************************************************/
copyMoleculeSE(MoleculeSE * m)566 MoleculeSE copyMoleculeSE(MoleculeSE* m)
567 {
568
569 gint i;
570 gint j;
571 MoleculeSE molecule = newMoleculeSE();
572
573 molecule.energy = m->energy;
574 molecule.nAtoms = m->nAtoms;
575 molecule.totalCharge = m->totalCharge;
576 molecule.spinMultiplicity = m->spinMultiplicity;
577 if( molecule.nAtoms>0) molecule.atoms = g_malloc(molecule.nAtoms*sizeof(AtomSE));
578 for(i=0;i<molecule.nAtoms;i++)
579 {
580 molecule.atoms[i].prop = prop_atom_get(m->atoms[i].prop.symbol);
581 for(j=0;j<3;j++) molecule.atoms[i].coordinates[j] = m->atoms[i].coordinates[j];
582 molecule.atoms[i].charge = m->atoms[i].charge;
583 molecule.atoms[i].mmType = g_strdup(m->atoms[i].mmType);
584 molecule.atoms[i].pdbType = g_strdup(m->atoms[i].pdbType);
585 molecule.atoms[i].residueName = g_strdup(m->atoms[i].residueName);
586 molecule.atoms[i].residueNumber = m->atoms[i].residueNumber;
587 molecule.atoms[i].N = m->atoms[i].N;
588 molecule.atoms[i].layer = m->atoms[i].layer;
589 molecule.atoms[i].show = m->atoms[i].show;
590 molecule.atoms[i].variable = m->atoms[i].variable;
591 if(m->atoms[i].typeConnections)
592 {
593 gint j;
594 molecule.atoms[i].typeConnections = g_malloc(molecule.nAtoms*sizeof(gint));
595 for(j=0;j<molecule.nAtoms;j++)
596 {
597 gint nj = m->atoms[j].N-1;
598 molecule.atoms[i].typeConnections[nj] = m->atoms[i].typeConnections[nj];
599 }
600 }
601 }
602
603
604 if(molecule.nAtoms>0)
605 for(j=0;j<3;j++) /* x, y and z derivatives */
606 {
607 molecule.gradient[j] = g_malloc(molecule.nAtoms*sizeof(gdouble));
608 for(i=0;i<molecule.nAtoms;i++)
609 molecule.gradient[j][i] = m->gradient[j][i];
610 }
611 for(i=0;i<3;i++) molecule.dipole[i] = m->dipole[i];
612
613 return molecule;
614 }
615 /******************************************************************************/
save_atom_hin_file(FILE * file,char * name,int atomNumber,char * atomPDBType,char * atomMMType,double x,double y,double z,char * symbol,double charge,int N,int * connection,int * connectionType)616 static void save_atom_hin_file(FILE* file, char* name, int atomNumber, char* atomPDBType, char* atomMMType,
617 double x, double y, double z, char* symbol, double charge,
618 int N, int* connection, int* connectionType)
619 {
620 int i;
621 fprintf(file,"%s %d ",name,atomNumber);
622 fprintf(file,"%s ",atomPDBType);
623 fprintf(file,"%s ",symbol);
624 fprintf(file,"%s - ",atomMMType);
625 fprintf(file,"%0.14f ",charge);
626 fprintf(file,"%0.14f ",x);
627 fprintf(file,"%0.14f ",y);
628 fprintf(file,"%0.14f ",z);
629 if(N>0)
630 {
631 fprintf(file,"%d ",N);
632 for(i=0;i<N;i++)
633 {
634 if(connectionType[i]==3) fprintf(file,"%d t ",connection[i]);
635 else if(connectionType[i]==2) fprintf(file,"%d d ",connection[i]);
636 else fprintf(file,"%d s ",connection[i]);
637 }
638 }
639 fprintf(file,"\n");
640 }
641 /******************************************************************************/
saveMoleculeSEHIN(MoleculeSE * mol,char * fileName)642 gboolean saveMoleculeSEHIN(MoleculeSE* mol, char* fileName)
643 {
644 int i,n,k;
645 FILE* file = fopen(fileName, "w");
646 int nAtoms;
647 AtomSE* atoms;
648 int* connection;
649 int* connectionType;
650 if(!file) return FALSE;
651 nAtoms = mol->nAtoms;
652 atoms = mol->atoms;
653
654 fprintf(file,"forcefield Amber99\n");
655 fprintf(file,"sys 0 0 1\n");
656 fprintf(file,"view 40 0.1272 55 15 0.247224 0.3713666 0.8949677 -0.8641704 0.5022867 0.0302929 -0.4382806 -0.7808937 0.4451014 6.191 0.64575 -54.754\n");
657 fprintf(file,"seed -1108\n");
658 fprintf(file,"mol 1\n");
659
660 connection = malloc(nAtoms*sizeof(int));
661 connectionType = malloc(nAtoms*sizeof(int));
662
663 for(i=0;i<nAtoms;i++)
664 {
665 n = 0;
666 if(atoms[i].typeConnections)
667 for(k=0;k<nAtoms;k++)
668 {
669 if(i==k) continue;
670 if(atoms[i].typeConnections[k]>0)
671 {
672 connection[n] = k+1;
673 connectionType[n] = atoms[i].typeConnections[k];
674 n++;
675 }
676 }
677 save_atom_hin_file(file,"ATOM",i+1,atoms[i].pdbType, atoms[i].mmType,
678 atoms[i].coordinates[0], atoms[i].coordinates[1], atoms[i].coordinates[2],
679 atoms[i].prop.symbol, atoms[i].charge,n,connection, connectionType);
680 }
681 fprintf(file,"endmol 1\n");
682 fclose(file);
683 free(connection);
684 free(connectionType);
685 return TRUE;
686 }
687
688 /*****************************************************************************/
saveMoleculeSEMol2(MoleculeSE * mol,char * fileName)689 gboolean saveMoleculeSEMol2(MoleculeSE* mol, char* fileName)
690 {
691 int i,n,j;
692 FILE* file = fopen(fileName, "w");
693 if(!file) return FALSE;
694 n = 0;
695 for(i=0;i<mol->nAtoms;i++)
696 if(mol->atoms[i].typeConnections)
697 for(j=i+1;j<mol->nAtoms;j++)
698 if(mol->atoms[i].typeConnections[j]) n++;
699
700 fprintf(file,"@<TRIPOS>MOLECULE\n");
701 fprintf(file,"MOL2 : Made in CChemI. mol2 file\n");
702 fprintf(file," %10d %10d %10d\n",mol->nAtoms,n,1);
703 fprintf(file,"SMALL\n");
704 fprintf(file,"GASTEIGER\n");
705 fprintf(file,"\n");
706 fprintf(file,"@<TRIPOS>ATOM\n");
707 for (i=0;i<mol->nAtoms;i++)
708 {
709 fprintf(file,"%7d%1s%-6s%12.4f%10.4f%10.4f%1s%-5s%4d%1s%-8s%10.4f\n",
710 i+1,"",mol->atoms[i].prop.symbol,
711 mol->atoms[i].coordinates[0],
712 mol->atoms[i].coordinates[1],
713 mol->atoms[i].coordinates[2],
714 "",mol->atoms[i].prop.symbol,1," ","LIG111",mol->atoms[i].charge);
715 }
716 fprintf(file,"@<TRIPOS>BOND\n");
717 n = 0;
718 for(i=0;i<mol->nAtoms;i++)
719 if(mol->atoms[i].typeConnections)
720 for(j=i+1;j<mol->nAtoms;j++)
721 if(mol->atoms[i].typeConnections[j])
722 {
723 n++;
724 fprintf(file,"%6d%6d%6d%3s%2d\n",n+1, i+1, j+1, "",mol->atoms[i].typeConnections[j]);
725 }
726
727
728 fclose(file);
729 return TRUE;
730 }
731 /*****************************************************************************/
computeMoleculeSEDipole(MoleculeSE * mol)732 void computeMoleculeSEDipole(MoleculeSE* mol)
733 {
734 int i,k;
735 for(k=0;k<3;k++) mol->dipole[k] = 0;
736 for(i=0;i<mol->nAtoms;i++)
737 for(k=0;k<3;k++)
738 mol->dipole[k] += mol->atoms[i].charge*mol->atoms[i].coordinates[k];
739 for(k=0;k<3;k++) mol->dipole[k] *= ANG_TO_BOHR*AUTODEB;
740 }
741
742 /********************************************************************************/
readGeomMoleculeSEFromOpenBabelOutputFile(MoleculeSE * mol,char * fileName,int numgeometry)743 void readGeomMoleculeSEFromOpenBabelOutputFile(MoleculeSE* mol, char* fileName, int numgeometry)
744 {
745 FILE* file = NULL;
746 char buffer[BSIZE];
747 char* pdest = NULL;
748 //char* energyTag = "TOTAL ENERGY =";
749 char* energyTag = "FINAL ENERGY:";
750 char* geomTag = "Geometry";
751 char* gradTag = "Gradients:";
752 double dum;
753
754 printf("Read geom from %s\n",fileName);
755 file = fopen(fileName, "r");
756 if(!file) return;
757 while(!feof(file))
758 {
759 if(!fgets(buffer,BSIZE,file))break;
760 pdest = strstr( buffer, energyTag);
761 if(pdest &&sscanf(pdest+strlen(energyTag)+1,"%lf",&mol->energy)==1)
762 {
763 if(strstr(pdest,"kJ")) mol->energy /= KCALTOKJ;
764 break;
765 }
766 }
767 while(!feof(file))
768 {
769 if(!fgets(buffer,BSIZE,file))break;
770 if(strstr(buffer, geomTag))
771 {
772 int i;
773 for(i=0;i<mol->nAtoms;i++)
774 {
775 if(!fgets(buffer,BSIZE,file))break;
776 //printf("%s\n",buffer);
777 if(sscanf(buffer,"%lf %lf %lf %lf",
778 &dum,
779 &mol->atoms[i].coordinates[0],
780 &mol->atoms[i].coordinates[1],
781 &mol->atoms[i].coordinates[2]
782 )!=4) break;
783 }
784 break;
785 }
786 }
787 while(!feof(file))
788 {
789 if(!fgets(buffer,BSIZE,file))break;
790 if(strstr(buffer, gradTag))
791 {
792 int i;
793 for(i=0;i<mol->nAtoms;i++)
794 {
795 if(!fgets(buffer,BSIZE,file))break;
796 //printf("%s\n",buffer);
797 if(sscanf(buffer,"%lf %lf %lf",
798 &mol->gradient[0][i],
799 &mol->gradient[1][i],
800 &mol->gradient[2][i]
801 )!=3) break;
802 }
803 break;
804 }
805 }
806 /*
807 {
808 int i;
809 for(i=0;i<mol->nAtoms;i++) mol->atoms[i].typeConnections = NULL;
810 }
811 */
812 fclose(file);
813 }
814 /********************************************************************************/
get_connections_one_atom(gchar * t,gint nAtoms,gint ibeg,gint * connections)815 gint get_connections_one_atom(gchar* t, gint nAtoms, gint ibeg, gint* connections)
816 {
817 gint k;
818 gint nc;
819 gint nj;
820 gchar** ssplit = NULL;
821 gint nA = 0;
822 /* int ibeg = 12;*/
823 for(k=0;k<nAtoms;k++) connections[k] = 0;
824 ssplit = gab_split(t);
825 nA = 0;
826 while(ssplit && ssplit[nA]!=NULL) nA++;
827 if(nA<ibeg)
828 {
829 gab_strfreev(ssplit);
830 return 0;
831 }
832 nc = atoi(ssplit[ibeg-1]);
833 for(k=0;k<2*nc;k+=2)
834 {
835 if(!ssplit[ibeg+k]) break;
836 if(!ssplit[ibeg+k+1]) break;
837 nj = atoi(ssplit[ibeg+k]);
838 connections[nj-1] = atoi(ssplit[ibeg+k+1]);
839 }
840
841 gab_strfreev(ssplit);
842
843 return 1;
844 }
845 /********************************************************************************/
readGeometryFromGenericOutputFile(MoleculeSE * molecule,char * namefile)846 gboolean readGeometryFromGenericOutputFile(MoleculeSE* molecule, char* namefile)
847 {
848 FILE* file = NULL;
849 static gchar *t = NULL;
850 gboolean Ok = FALSE;
851 gint n,ic,is;
852 MoleculeSE mols = newMoleculeSE();
853 MoleculeSE* mol = &mols;
854 #define SZ 50
855 gchar symbol[SZ];
856 gchar mmType[SZ];
857 gchar pdbType[SZ];
858 gchar residueName[SZ];
859 gdouble X,Y,Z;
860 gdouble charge;
861 gint layer;
862 gint i,j,l;
863 gchar* pos;
864
865 if(!namefile)
866 {
867 printf("Sorry I cannot read geometry namefile = NULL\n");
868 exit(1);
869 }
870 if(t==NULL) t = malloc(BSIZE*sizeof(char));
871 file = fopen(namefile,"rb");
872 if(!file)
873 {
874 printf("Sorry I cannot open %s file\n",namefile);
875 exit(1);
876 }
877 rewind(file);
878 while(!feof(file))
879 {
880 if(!fgets(t,BSIZE, file)) break;
881 deleteFirstSpaces(t);
882 if(t[0]=='#') continue;
883 uppercase(t);
884 pos = strstr(t,"GEOMETRY");
885 if(pos)
886 {
887 while(!feof(file))
888 {
889 if(!fgets(t,BSIZE, file)) break;
890 deleteFirstSpaces(t);
891 if(t[0]=='#') continue;
892 break;
893 }
894 if(3==sscanf(t,"%d%d%d",&n,&ic,&is) && n>0 && is>0)
895 {
896 /*
897 printf("Mult = %d\n",is);
898 printf("nAtoms = %d\n",n);
899 */
900 mol->nAtoms = n;
901 mol->spinMultiplicity = is;
902 mol->totalCharge = ic;
903 mol->atoms = g_malloc(mol->nAtoms*sizeof(AtomSE));
904 for(i=0; i<mol->nAtoms; i++) mol->atoms[i].typeConnections = g_malloc(mol->nAtoms*sizeof(gint));
905 Ok = TRUE;
906 }
907 break;
908 }
909 }
910 if(!Ok)
911 {
912 printf("Sorry I cannot read geometry from %s file\n",namefile);
913 exit(1);
914 }
915 for(i=0; i<mol->nAtoms; i++)
916 {
917 int variable = 0;
918 gint ibeg = 12;
919 if(!fgets(t,BSIZE,file))
920 {
921 printf("Sorry I cannot read geometry from %s file.\n",namefile);
922 exit(1);
923 }
924 deleteFirstSpaces(t);
925 if(t[0]=='#') { i--;continue;}
926 sscanf(t,"%s %s %s %s %d %lf %d %d %lf %lf %lf",
927 symbol,mmType,pdbType,residueName,
928 &mol->atoms[i].residueNumber,
929 &charge,&layer,&variable,&X,&Y,&Z);
930 symbol[0]=toupper(symbol[0]);
931 l=strlen(symbol);
932 if (l==2) symbol[1]=tolower(symbol[1]);
933
934 mol->atoms[i].prop = prop_atom_get(symbol);
935 mol->atoms[i].mmType=strdup(mmType);
936 mol->atoms[i].pdbType=strdup(pdbType);
937 mol->atoms[i].residueName=strdup(residueName);
938 mol->atoms[i].N=i+1;
939 mol->atoms[i].layer=layer;
940 mol->atoms[i].variable=variable;
941 mol->atoms[i].show=TRUE;
942 mol->atoms[i].coordinates[0] = X;
943 mol->atoms[i].coordinates[1] = Y;
944 mol->atoms[i].coordinates[2] = Z;
945 mol->atoms[i].charge = charge;
946 if(!get_connections_one_atom(t, mol->nAtoms, ibeg, mol->atoms[i].typeConnections))
947 {
948 printf("Sorry I cannot read the connection for atom # %d from the %s file.\n",i+1,namefile);
949 exit(1);
950 }
951 }
952 fclose(file);
953 if(mol->nAtoms>0)
954 for(j=0;j<3;j++) /* x, y and z derivatives */
955 {
956 mol->gradient[j] = g_malloc(mol->nAtoms*sizeof(gdouble));
957 for(i=0;i<mol->nAtoms;i++) mol->gradient[j][i] = 0.0;
958 }
959 for(i=0;i<3;i++) mol->dipole[i] = 0.0;
960 /*
961 printf("Begin copyMol\n");
962 *molecule = copyMoleculeSE(mol);
963 printf("End copyMol\n");
964 */
965 if(molecule->nAtoms==mol->nAtoms)
966 {
967 for(j=0;j<3;j++)
968 for(i=0;i<mol->nAtoms;i++)
969 molecule->atoms[i].coordinates[j] = mol->atoms[i].coordinates[j];
970
971 }
972 freeMoleculeSE(mol);
973 return TRUE;
974 }
975 /*****************************************************************************/
addGeometryMoleculeSEToGabedit(MoleculeSE * molecule,FILE * file)976 gboolean addGeometryMoleculeSEToGabedit(MoleculeSE* molecule,FILE* file)
977 {
978 int j,k;
979 int nc;
980
981 if(!molecule) return FALSE;
982 fprintf(file,"%d %d %d\n",molecule->nAtoms, molecule->totalCharge, molecule->spinMultiplicity);
983 for(j=0;j<molecule->nAtoms;j++)
984 {
985 nc = 0;
986 for(k=0;k<molecule->nAtoms;k++) if(molecule->atoms[j].typeConnections[k]>0) nc++;
987 fprintf(file," %s %s %s %s %d %0.12lf %d %d %0.12lf %0.12lf %0.12lf %d ",
988 molecule->atoms[j].prop.symbol,
989 molecule->atoms[j].mmType,
990 molecule->atoms[j].pdbType,
991 molecule->atoms[j].residueName,
992 molecule->atoms[j].residueNumber,
993 molecule->atoms[j].charge,
994 molecule->atoms[j].layer,
995 molecule->atoms[j].variable,
996 molecule->atoms[j].coordinates[0],
997 molecule->atoms[j].coordinates[1],
998 molecule->atoms[j].coordinates[2],
999 nc
1000 );
1001 for(k=0;k<molecule->nAtoms;k++)
1002 {
1003 int nk = molecule->atoms[k].N-1;
1004 if(molecule->atoms[j].typeConnections[nk]>0)
1005 fprintf(file," %d %d", nk+1,molecule->atoms[j].typeConnections[nk]);
1006 }
1007 fprintf(file," GRADIENT %0.14f %0.14f %0.14f",molecule->gradient[0][j], molecule->gradient[1][j],molecule->gradient[2][j]);
1008 fprintf(file,"\n");
1009 }
1010 return TRUE;
1011
1012 }
1013 /*****************************************************************************/
addMoleculeSEToFile(MoleculeSE * molecule,FILE * file)1014 gboolean addMoleculeSEToFile(MoleculeSE* molecule,FILE* file)
1015 {
1016 int j,k;
1017 int nc;
1018
1019 if(!molecule) return FALSE;
1020
1021 fprintf(file,"Geometry\n");
1022 return addGeometryMoleculeSEToGabedit(molecule,file);
1023 }
1024 /*****************************************************************************/
saveMoleculeSETypeSave(MoleculeSE * molecule,char * fileName,char * typeSave)1025 gboolean saveMoleculeSETypeSave(MoleculeSE* molecule, char* fileName, char* typeSave)
1026 {
1027 FILE* file = NULL;
1028 int j;
1029 int form = 1;
1030
1031 printf("Save molecule in %s\n",fileName);
1032 if(!molecule) return FALSE;
1033
1034 file = fopen(fileName, typeSave);
1035
1036 if(!file) return FALSE;
1037
1038 fprintf(file,"[Gabedit Format]\n");
1039 fprintf(file,"[GEOCONV]\n");
1040 fprintf(file,"energy\n");
1041 fprintf(file,"%f\n",molecule->energy);
1042 fprintf(file,"max-force\n");
1043 fprintf(file,"%f\n",0.0);
1044 fprintf(file,"rms-force\n");
1045 fprintf(file,"%f\n",0.0);
1046
1047 fprintf(file,"\n");
1048 fprintf(file,"[GEOMETRIES]\n");
1049 {
1050 fprintf(file,"%d\n",molecule->nAtoms);
1051 fprintf(file,"\n");
1052 for(j=0;j<molecule->nAtoms;j++)
1053 fprintf(file," %s %0.8f %0.8f %0.8f\n",
1054 molecule->atoms[j].prop.symbol,
1055 molecule->atoms[j].coordinates[0],
1056 molecule->atoms[j].coordinates[1],
1057 molecule->atoms[j].coordinates[2]
1058 );
1059 }
1060 fprintf(file,"\n");
1061 fprintf(file,"[GEOMS] %d\n",form);
1062 fprintf(file,"%d 3\n",1);
1063 fprintf(file,"energy kcal/mol 1\n");
1064 fprintf(file,"deltaE K 1\n");
1065 fprintf(file,"Dipole Debye 3\n");
1066 //molecule->klass->computeDipole(molecule);
1067 {
1068 fprintf(file,"%0.14f\n",molecule->energy);
1069 fprintf(file,"0\n");
1070 fprintf(file,"%0.14f %0.14f %0.14f\n",molecule->dipole[0],molecule->dipole[1],molecule->dipole[2]);
1071 addGeometryMoleculeSEToGabedit(molecule,file);
1072 }
1073 //addVibrationToFile(molecule, file);
1074 fclose(file);
1075 return TRUE;
1076
1077 }
1078 /*****************************************************************************/
saveMoleculeSE(MoleculeSE * molecule,char * fileName)1079 gboolean saveMoleculeSE(MoleculeSE* molecule, char* fileName)
1080 {
1081 return saveMoleculeSETypeSave(molecule, fileName, "w");
1082 }
1083 /*****************************************************************************/
getGradientNormMoleculeSE(MoleculeSE * molecule)1084 gdouble getGradientNormMoleculeSE(MoleculeSE* molecule)
1085 {
1086 int j,k;
1087 gdouble norm = 0;
1088
1089 if(!molecule) return -1.0;
1090 for(j=0;j<molecule->nAtoms;j++)
1091 for(k=0;k<3;k++)
1092 norm += molecule->gradient[k][j]*molecule->gradient[k][j];
1093 return sqrt(norm);
1094
1095 }
1096