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