1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <glib.h>
5 #include <math.h>
6
7 typedef struct _AmberParameters AmberParameters;
8 typedef struct _AmberAtomTypes AmberAtomTypes;
9 typedef struct _AmberBondStretchTerms AmberBondStretchTerms;
10 typedef struct _AmberAngleBendTerms AmberAngleBendTerms;
11 typedef struct _AmberDihedralAngleTerms AmberDihedralAngleTerms;
12 typedef struct _AmberImproperTorsionTerms AmberImproperTorsionTerms;
13 typedef struct _AmberNonBondedTerms AmberNonBondedTerms;
14 typedef struct _AmberHydrogenBondedTerms AmberHydrogenBondedTerms;
15
16 /************************************/
17 struct _AmberAtomTypes
18 {
19 gchar* name;
20 gint number;
21 gdouble masse;
22 gdouble polarisability;
23 gchar* c;
24 };
25 /************************************/
26 struct _AmberBondStretchTerms
27 {
28 gint numbers[2];
29 gdouble equilibriumDistance;
30 gdouble forceConstant;
31 gchar* c;
32 };
33 /************************************/
34 struct _AmberAngleBendTerms
35 {
36 gint numbers[3];
37 gdouble equilibriumAngle;
38 gdouble forceConstant;
39 gchar* c;
40 };
41 /************************************/
42 struct _AmberDihedralAngleTerms
43 {
44 gint numbers[4];
45 gint nSomme;
46 gdouble* divisor;
47 gdouble* barrier;
48 gdouble* phase;
49 gdouble* n;
50 gchar* c;
51 };
52 /************************************/
53 struct _AmberImproperTorsionTerms
54 {
55 gint numbers[4];
56 gdouble barrier;
57 gdouble phase;
58 gdouble n;
59 gchar* c;
60 };
61 /************************************/
62 struct _AmberNonBondedTerms
63 {
64 gint number;
65 gdouble r;
66 gdouble epsilon;
67 gchar* c;
68 };
69 /************************************/
70 struct _AmberHydrogenBondedTerms
71 {
72 gint numbers[2];
73 gdouble c;
74 gdouble d;
75 gchar* cc;
76 };
77 /************************************/
78 struct _AmberParameters
79 {
80 gint numberOfTypes;
81 AmberAtomTypes* atomTypes;
82
83 gint numberOfStretchTerms;
84 AmberBondStretchTerms* bondStretchTerms;
85
86 gint numberOfBendTerms;
87 AmberAngleBendTerms* angleBendTerms;
88
89 gint numberOfDihedralTerms;
90 AmberDihedralAngleTerms* dihedralAngleTerms;
91
92 gint numberOfImproperTorsionTerms;
93 AmberImproperTorsionTerms* improperTorsionTerms;
94
95 gint numberOfNonBonded;
96 AmberNonBondedTerms* nonBondedTerms;
97
98 gint numberOfHydrogenBonded;
99 AmberHydrogenBondedTerms* hydrogenBondedTerms;
100 };
101 /************************************/
102
103 static gchar atomTypesTitle[] = "Begin INPUT FOR ATOM TYPES, MASSE AND POLARISABILITIES";
104 static gchar bondStretchTitle[] = "Begin INPUT FOR BOND LENGTH PARAMETERS";
105 static gchar angleBendTitle[] = "Begin INPUT FOR BOND ANGLE PARAMETERS";
106 static gchar hydrogenBondedTitle[] = "Begin INPUT FOR H-BOND 10-12 POTENTIAL PARAMETERS";
107 static gchar improperTorsionTitle[] ="Begin INPUT FOR IMPROPER DIHEDRAL PARAMETERS";
108 static gchar nonBondedTitle[] ="Begin INPUT FOR THE NON-BONDED 6-12 POTENTIAL PARAMETERS";
109 static gchar nonBondedEquivTitle[] ="Begin INPUT FOR EQUIVALENCING ATOM TYPES FOR THE NON-BONDED 6-12 POTENTIAL PARAMETERS";
110 static gchar dihedralAngleTitle[] = "Begin INPUT FOR DIHEDRAL PARAMETERS";
111
112
113 /**********************************************************************/
newAmberParameters()114 AmberParameters newAmberParameters()
115 {
116 AmberParameters amberParameters;
117
118 amberParameters.numberOfTypes = 0;
119 amberParameters.atomTypes = NULL;
120
121 amberParameters.numberOfStretchTerms = 0;
122 amberParameters.bondStretchTerms = NULL;
123
124 amberParameters.numberOfBendTerms = 0;
125 amberParameters.angleBendTerms = NULL;
126
127 amberParameters.numberOfDihedralTerms = 0;
128 amberParameters.dihedralAngleTerms = NULL;
129
130 amberParameters.numberOfImproperTorsionTerms = 0;
131 amberParameters.improperTorsionTerms = NULL;
132
133 amberParameters.numberOfNonBonded = 0;
134 amberParameters.nonBondedTerms = NULL;
135
136 amberParameters.numberOfHydrogenBonded = 0;
137 amberParameters.hydrogenBondedTerms = NULL;
138
139
140 return amberParameters;
141
142 }
143 /**********************************************************************/
freeAmberParameters(AmberParameters * amberParameters)144 void freeAmberParameters(AmberParameters* amberParameters)
145 {
146 gint i;
147
148 for(i=0;i<amberParameters->numberOfTypes;i++)
149 if(amberParameters->atomTypes[i].name)
150 g_free(amberParameters->atomTypes[i].name);
151
152 amberParameters->numberOfTypes = 0;
153 if(amberParameters->atomTypes )
154 g_free(amberParameters->atomTypes );
155 amberParameters->atomTypes = NULL;
156
157 amberParameters->numberOfStretchTerms = 0;
158 if(amberParameters->bondStretchTerms)
159 g_free(amberParameters->bondStretchTerms);
160 amberParameters->bondStretchTerms = NULL;
161
162 amberParameters->numberOfBendTerms = 0;
163 if(amberParameters->angleBendTerms)
164 g_free(amberParameters->angleBendTerms);
165 amberParameters->angleBendTerms = NULL;
166
167 for(i=0;i<amberParameters->numberOfDihedralTerms;i++)
168 {
169 if(amberParameters->dihedralAngleTerms[i].divisor)
170 g_free(amberParameters->dihedralAngleTerms[i].divisor);
171 if(amberParameters->dihedralAngleTerms[i].barrier)
172 g_free(amberParameters->dihedralAngleTerms[i].barrier);
173 if(amberParameters->dihedralAngleTerms[i].phase)
174 g_free(amberParameters->dihedralAngleTerms[i].phase);
175 if(amberParameters->dihedralAngleTerms[i].n)
176 g_free(amberParameters->dihedralAngleTerms[i].n);
177
178 }
179
180 amberParameters->numberOfDihedralTerms = 0;
181 if(amberParameters->dihedralAngleTerms)
182 g_free(amberParameters->dihedralAngleTerms);
183 amberParameters->dihedralAngleTerms = NULL;
184
185 amberParameters->numberOfImproperTorsionTerms = 0;
186 if(amberParameters->improperTorsionTerms)
187 g_free(amberParameters->improperTorsionTerms);
188 amberParameters->improperTorsionTerms = NULL;
189
190 amberParameters->numberOfNonBonded = 0;
191 if(amberParameters->nonBondedTerms)
192 g_free(amberParameters->nonBondedTerms);
193 amberParameters->nonBondedTerms = NULL;
194
195 amberParameters->numberOfHydrogenBonded = 0;
196 if(amberParameters->hydrogenBondedTerms)
197 g_free(amberParameters->hydrogenBondedTerms);
198 amberParameters->hydrogenBondedTerms = NULL;
199 }
200 /**********************************************************************/
getNumberType(AmberParameters * amberParameters,gchar * type)201 gint getNumberType(AmberParameters* amberParameters, gchar* type)
202 {
203 gint i;
204 gint nTypes = amberParameters->numberOfTypes;
205 AmberAtomTypes* types = amberParameters->atomTypes;
206 gint len = strlen(type);
207
208 if(strcmp(type,"X")==0)
209 return -1;
210 for(i=0;i<nTypes;i++)
211 {
212 if(len == (gint)strlen(types[i].name) && strstr(types[i].name,type))
213 return types[i].number;
214
215 }
216 return -2;
217 }
218 /**********************************************************************/
readAmberTypes(AmberParameters * amberParameters,FILE * file)219 gboolean readAmberTypes(AmberParameters* amberParameters, FILE* file)
220 {
221 gchar t[1024];
222
223 gchar dump[1024];
224 gint len = 1024;
225 gboolean Ok = FALSE;
226 gint n = 0;
227 AmberAtomTypes* types = NULL;
228 /* Search Begin INPUT FOR ATOM TYPES */
229
230 while(!feof(file))
231 {
232 if(fgets(t,len,file))
233 {
234 if(strstr(t,atomTypesTitle))
235 {
236 Ok = TRUE;
237 break;
238 }
239 }
240 }
241 if(!Ok)
242 return FALSE;
243
244 types = g_malloc(sizeof(AmberAtomTypes));
245 n = 0;
246 Ok = FALSE;
247 while(!feof(file))
248 {
249 if(fgets(t,len,file))
250 {
251 if(strstr(t,"End"))
252 {
253 Ok = TRUE;
254 break;
255 }
256 }
257 else
258 {
259 Ok = FALSE;
260 break;
261 }
262
263
264
265 sscanf(t,"%s %lf %lf",dump,&types[n].masse,&types[n].polarisability);
266 types[n].name = g_strdup(dump);
267 types[n].c = g_strdup(t+37);
268 types[n].c[strlen(types[n].c)-1]='\0';
269 types[n].number = n;
270
271 n++;
272 types = g_realloc(types,(n+1)*sizeof(AmberAtomTypes));
273 }
274 if(n==0 || !Ok )
275 g_free(types);
276 else
277 {
278 amberParameters->numberOfTypes = n;
279 amberParameters->atomTypes = types;
280 }
281 /* printing for test*/
282 /*
283 printf("umber of types = %d \n",amberParameters->numberOfTypes);
284 for(n=0;n<amberParameters->numberOfTypes;n++)
285 {
286 printf("%s\t %d\t",
287 amberParameters->atomTypes[n].name,
288 amberParameters->atomTypes[n].number
289 );
290 }
291 printf("\n");
292 */
293
294 return TRUE;
295
296
297 }
298 /**********************************************************************/
readAmberBondStretchTerms(AmberParameters * amberParameters,FILE * file)299 gboolean readAmberBondStretchTerms(AmberParameters* amberParameters,FILE* file)
300 {
301 gchar t[1024];
302 gchar dump1[1024];
303 gchar dump2[1024];
304 gint len = 1024;
305 gboolean Ok = FALSE;
306 gint n = 0;
307 AmberBondStretchTerms* terms = NULL;
308
309 /* Search Begin INPUT FOR ATOM TYPES */
310
311 while(!feof(file))
312 {
313 if(fgets(t,len,file))
314 {
315 if(strstr(t,bondStretchTitle))
316 {
317 Ok = TRUE;
318 break;
319 }
320 }
321 }
322 if(!Ok)
323 return FALSE;
324
325 terms = g_malloc(sizeof(AmberBondStretchTerms));
326 n = 0;
327 Ok = FALSE;
328 while(!feof(file))
329 {
330 if(fgets(t,len,file))
331 {
332 if(strstr(t,"End"))
333 {
334 Ok = TRUE;
335 break;
336 }
337 }
338 else
339 {
340 Ok = FALSE;
341 break;
342 }
343
344
345
346 dump1[2]='\0';
347 dump2[2]='\0';
348 sscanf(t,"%c%c-%c%c %lf %lf",&dump1[0],&dump1[1],&dump2[0],&dump2[1],
349 &terms[n].forceConstant,
350 &terms[n].equilibriumDistance);
351 if(strlen(t)>28)
352 terms[n].c = g_strdup(t+28);
353 else
354 terms[n].c = g_strdup(" ");
355
356 terms[n].c[strlen(terms[n].c)-1]='\0';
357 if(dump1[1]==' ') dump1[1]='\0';
358 if(dump2[1]==' ') dump2[1]='\0';
359 terms[n].numbers[0] = getNumberType(amberParameters,dump1);
360 terms[n].numbers[1] = getNumberType(amberParameters,dump2);
361
362 n++;
363 terms = g_realloc(terms,(n+1)*sizeof(AmberBondStretchTerms));
364 }
365 if(n==0 || !Ok )
366 g_free(terms);
367 else
368 {
369 amberParameters->numberOfStretchTerms = n;
370 amberParameters->bondStretchTerms = terms;
371 }
372 /* printing for test*/
373 /*
374 printf("number of bonds = %d \n",amberParameters->numberOfStretchTerms);
375 for(n=0;n<amberParameters->numberOfStretchTerms;n++)
376 {
377 printf("%d %d %f %f\n",
378 amberParameters->bondStretchTerms[n].numbers[0],
379 amberParameters->bondStretchTerms[n].numbers[1],
380 amberParameters->bondStretchTerms[n].forceConstant,
381 amberParameters->bondStretchTerms[n].equilibriumDistance
382 );
383 }
384 printf("\n");
385 */
386
387 return TRUE;
388
389
390 }
391 /**********************************************************************/
readAmberAngleBendTerms(AmberParameters * amberParameters,FILE * file)392 gboolean readAmberAngleBendTerms(AmberParameters* amberParameters,FILE* file)
393 {
394 gchar t[1024];
395 gchar dump1[10];
396 gchar dump2[10];
397 gchar dump3[10];
398 gint len = 1024;
399 gboolean Ok = FALSE;
400 gint n = 0;
401 AmberAngleBendTerms* terms = NULL;
402
403 /* Search Begin INPUT FOR ATOM TYPES */
404
405 while(!feof(file))
406 {
407 if(fgets(t,len,file))
408 {
409 if(strstr(t,angleBendTitle))
410 {
411 Ok = TRUE;
412 break;
413 }
414 }
415 }
416 if(!Ok)
417 return FALSE;
418
419 terms = g_malloc(sizeof(AmberAngleBendTerms));
420 n = 0;
421 Ok = FALSE;
422 while(!feof(file))
423 {
424 if(fgets(t,len,file))
425 {
426 if(strstr(t,"End"))
427 {
428 Ok = TRUE;
429 break;
430 }
431 }
432 else
433 {
434 Ok = FALSE;
435 break;
436 }
437
438
439
440 dump1[2]='\0';
441 dump2[2]='\0';
442 dump3[2]='\0';
443 sscanf(t,"%c%c-%c%c-%c%c %lf %lf",
444 &dump1[0],&dump1[1],&dump2[0],&dump2[1],&dump3[0],&dump3[1],
445 &terms[n].forceConstant,
446 &terms[n].equilibriumAngle);
447
448 if(strlen(t)>32)
449 terms[n].c = g_strdup(t+32);
450 else
451 terms[n].c = g_strdup(" ");
452
453 terms[n].c[strlen(terms[n].c)-1]='\0';
454
455 if(dump1[1]==' ') dump1[1]='\0';
456 if(dump2[1]==' ') dump2[1]='\0';
457 if(dump3[1]==' ') dump3[1]='\0';
458
459 terms[n].numbers[0] = getNumberType(amberParameters,dump1);
460 terms[n].numbers[1] = getNumberType(amberParameters,dump2);
461 terms[n].numbers[2] = getNumberType(amberParameters,dump3);
462
463 n++;
464 terms = g_realloc(terms,(n+1)*sizeof(AmberAngleBendTerms));
465 }
466 if(n==0 || !Ok )
467 g_free(terms);
468 else
469 {
470 amberParameters->numberOfBendTerms = n;
471 amberParameters->angleBendTerms = terms;
472 }
473 /* printing for test*/
474 /*
475 printf("number of bonds = %d \n",amberParameters->numberOfBendTerms);
476 for(n=0;n<amberParameters->numberOfBendTerms;n++)
477 {
478 printf("%d %d %d %f %f\n",
479 amberParameters->angleBendTerms[n].numbers[0],
480 amberParameters->angleBendTerms[n].numbers[1],
481 amberParameters->angleBendTerms[n].numbers[2],
482 amberParameters->angleBendTerms[n].forceConstant,
483 amberParameters->angleBendTerms[n].equilibriumAngle
484 );
485 }
486 printf("\n");
487 */
488
489 return TRUE;
490
491
492 }
493 /**********************************************************************/
readAmberDihedralAngleTerms(AmberParameters * amberParameters,FILE * file)494 gboolean readAmberDihedralAngleTerms(AmberParameters* amberParameters,FILE* file)
495 {
496 gchar t[1024];
497 gchar dump1[10];
498 gchar dump2[10];
499 gchar dump3[10];
500 gchar dump4[10];
501 gint len = 1024;
502 gboolean Ok = FALSE;
503 gint n = 0;
504 AmberDihedralAngleTerms* terms = NULL;
505 gdouble divisor = 1;
506 gdouble barrier = 0;
507 gdouble phase = 0;
508 gdouble nN = 0;
509
510 /* Search Begin INPUT FOR DIHEDRAL PARAMETERS */
511
512 while(!feof(file))
513 {
514 if(fgets(t,len,file))
515 {
516 if(strstr(t,dihedralAngleTitle))
517 {
518 Ok = TRUE;
519 break;
520 }
521 }
522 }
523 if(!Ok)
524 return FALSE;
525
526 terms = g_malloc(sizeof(AmberDihedralAngleTerms));
527 n = 0;
528 Ok = FALSE;
529 while(!feof(file))
530 {
531 if(fgets(t,len,file))
532 {
533 if(strstr(t,"End"))
534 {
535 Ok = TRUE;
536 break;
537 }
538 }
539 else
540 {
541 Ok = FALSE;
542 break;
543 }
544
545
546 terms[n].nSomme = 1;
547 terms[n].divisor = g_malloc(sizeof(gdouble));
548 terms[n].barrier = g_malloc(sizeof(gdouble));
549 terms[n].phase = g_malloc(sizeof(gdouble));
550 terms[n].n = g_malloc(sizeof(gdouble));
551
552 dump1[2]='\0';
553 dump2[2]='\0';
554 dump3[2]='\0';
555 dump4[2]='\0';
556 sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf %lf",
557 &dump1[0],&dump1[1],
558 &dump2[0],&dump2[1],
559 &dump3[0],&dump3[1],
560 &dump4[0],&dump4[1],
561 &divisor,
562 &barrier,
563 &phase,
564 &nN);
565 terms[n].divisor[0] = divisor;
566 terms[n].barrier[0] = barrier;
567 terms[n].phase[0] = phase;
568 terms[n].n[0] = fabs(nN);
569
570 if(strlen(t)>60)
571 terms[n].c = g_strdup(t+60);
572 else
573 terms[n].c = g_strdup(" ");
574
575 terms[n].c[strlen(terms[n].c)-1]='\0';
576
577
578 if(dump1[1]==' ') dump1[1]='\0';
579 if(dump2[1]==' ') dump2[1]='\0';
580 if(dump3[1]==' ') dump3[1]='\0';
581 if(dump4[1]==' ') dump4[1]='\0';
582
583 terms[n].numbers[0] = getNumberType(amberParameters,dump1);
584 terms[n].numbers[1] = getNumberType(amberParameters,dump2);
585 terms[n].numbers[2] = getNumberType(amberParameters,dump3);
586 terms[n].numbers[3] = getNumberType(amberParameters,dump4);
587
588 Ok = TRUE;
589 while(!feof(file) && nN<0)
590 {
591 if(!fgets(t,len,file))
592 {
593 Ok = FALSE;
594 break;
595 }
596
597 terms[n].nSomme++;
598 terms[n].divisor = g_realloc(terms[n].divisor,terms[n].nSomme*sizeof(gdouble));
599 terms[n].barrier = g_realloc(terms[n].barrier,terms[n].nSomme*sizeof(gdouble));
600 terms[n].phase = g_realloc(terms[n].phase,terms[n].nSomme*sizeof(gdouble));
601 terms[n].n = g_realloc(terms[n].n,terms[n].nSomme*sizeof(gdouble));
602
603 sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf %lf",
604 &dump1[0],&dump1[1],
605 &dump2[0],&dump2[1],
606 &dump3[0],&dump3[1],
607 &dump4[0],&dump4[1],
608 &divisor,
609 &barrier,
610 &phase,
611 &nN);
612 terms[n].divisor[terms[n].nSomme-1] = divisor;
613 terms[n].barrier[terms[n].nSomme-1] = barrier;
614 terms[n].phase[terms[n].nSomme-1] = phase;
615 terms[n].n[terms[n].nSomme-1] = fabs(nN);
616 }
617 if(!Ok)
618 break;
619 n++;
620 terms = g_realloc(terms,(n+1)*sizeof(AmberDihedralAngleTerms));
621 }
622 if(n==0 || !Ok )
623 g_free(terms);
624 else
625 {
626 amberParameters->numberOfDihedralTerms = n;
627 amberParameters->dihedralAngleTerms = terms;
628 }
629 /* printing for test*/
630 /*
631 printf("number of dihedral torsion terms = %d \n",
632 amberParameters->numberOfDihedralTerms);
633
634 for(n=0;n<amberParameters->numberOfDihedralTerms;n++)
635 {
636 gint j;
637 printf("%d %d %d %d \t",
638 amberParameters->dihedralAngleTerms[n].numbers[0],
639 amberParameters->dihedralAngleTerms[n].numbers[1],
640 amberParameters->dihedralAngleTerms[n].numbers[2],
641 amberParameters->dihedralAngleTerms[n].numbers[3]
642 );
643 for(j=0;j<amberParameters->dihedralAngleTerms[n].nSomme;j++)
644 {
645 printf("%f %f %f %f\t",
646 amberParameters->dihedralAngleTerms[n].divisor[j],
647 amberParameters->dihedralAngleTerms[n].barrier[j],
648 amberParameters->dihedralAngleTerms[n].phase[j],
649 amberParameters->dihedralAngleTerms[n].n[j]
650 );
651 }
652 printf("\n");
653 }
654 printf("\n");
655 */
656
657 return TRUE;
658
659
660 }
661 /**********************************************************************/
readAmberImproperTorsionTerms(AmberParameters * amberParameters,FILE * file)662 gboolean readAmberImproperTorsionTerms(AmberParameters* amberParameters,FILE* file)
663 {
664 gchar t[1024];
665 gchar dump1[10];
666 gchar dump2[10];
667 gchar dump3[10];
668 gchar dump4[10];
669 gint len = 1024;
670 gboolean Ok = FALSE;
671 gint n = 0;
672 AmberImproperTorsionTerms* terms = NULL;
673
674 /* Search Begin INPUT FOR ATOM TYPES */
675
676 while(!feof(file))
677 {
678 if(fgets(t,len,file))
679 {
680 if(strstr(t,improperTorsionTitle))
681 {
682 Ok = TRUE;
683 break;
684 }
685 }
686 }
687 if(!Ok)
688 return FALSE;
689
690 terms = g_malloc(sizeof(AmberImproperTorsionTerms));
691 n = 0;
692 Ok = FALSE;
693 while(!feof(file))
694 {
695 if(fgets(t,len,file))
696 {
697 if(strstr(t,"End"))
698 {
699 Ok = TRUE;
700 break;
701 }
702 }
703 else
704 {
705 Ok = FALSE;
706 break;
707 }
708
709
710
711 dump1[2]='\0';
712 dump2[2]='\0';
713 dump3[2]='\0';
714 dump4[2]='\0';
715 sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf",
716 &dump1[0],&dump1[1],
717 &dump2[0],&dump2[1],
718 &dump3[0],&dump3[1],
719 &dump4[0],&dump4[1],
720 &terms[n].barrier,
721 &terms[n].phase,
722 &terms[n].n);
723
724 if(strlen(t)>60)
725 terms[n].c = g_strdup(t+60);
726 else
727 terms[n].c = g_strdup(" ");
728
729 terms[n].c[strlen(terms[n].c)-1]='\0';
730
731 if(dump1[1]==' ') dump1[1]='\0';
732 if(dump2[1]==' ') dump2[1]='\0';
733 if(dump3[1]==' ') dump3[1]='\0';
734 if(dump4[1]==' ') dump4[1]='\0';
735
736 terms[n].numbers[0] = getNumberType(amberParameters,dump1);
737 terms[n].numbers[1] = getNumberType(amberParameters,dump2);
738 terms[n].numbers[2] = getNumberType(amberParameters,dump3);
739 terms[n].numbers[3] = getNumberType(amberParameters,dump4);
740
741 n++;
742 terms = g_realloc(terms,(n+1)*sizeof(AmberImproperTorsionTerms));
743 }
744 if(n==0 || !Ok )
745 g_free(terms);
746 else
747 {
748 amberParameters->numberOfImproperTorsionTerms = n;
749 amberParameters->improperTorsionTerms = terms;
750 }
751 /* printing for test*/
752 /*
753 printf("number of improper torsion terms = %d \n",
754 amberParameters->numberOfImproperTorsionTerms);
755
756 for(n=0;n<amberParameters->numberOfImproperTorsionTerms;n++)
757 {
758 printf("%d %d %d %d %f %f %f\n",
759 amberParameters->improperTorsionTerms[n].numbers[0],
760 amberParameters->improperTorsionTerms[n].numbers[1],
761 amberParameters->improperTorsionTerms[n].numbers[2],
762 amberParameters->improperTorsionTerms[n].numbers[3],
763 amberParameters->improperTorsionTerms[n].barrier,
764 amberParameters->improperTorsionTerms[n].phase,
765 amberParameters->improperTorsionTerms[n].n
766 );
767 }
768 printf("\n");
769 */
770
771 return TRUE;
772
773
774 }
775 /**********************************************************************/
readAmberHydrogenBondedTerms(AmberParameters * amberParameters,FILE * file)776 gboolean readAmberHydrogenBondedTerms(AmberParameters* amberParameters,FILE* file)
777 {
778 gchar t[1024];
779 gchar dump1[10];
780 gchar dump2[10];
781 gint len = 1024;
782 gboolean Ok = FALSE;
783 gint n = 0;
784 AmberHydrogenBondedTerms* terms = NULL;
785
786 /* Search Begin INPUT FOR ATOM TYPES */
787
788 while(!feof(file))
789 {
790 if(fgets(t,len,file))
791 {
792 if(strstr(t,hydrogenBondedTitle))
793 {
794 Ok = TRUE;
795 break;
796 }
797 }
798 }
799 if(!Ok)
800 return FALSE;
801
802 terms = g_malloc(sizeof(AmberHydrogenBondedTerms));
803 n = 0;
804 Ok = FALSE;
805 while(!feof(file))
806 {
807 if(fgets(t,len,file))
808 {
809 if(strstr(t,"End"))
810 {
811 Ok = TRUE;
812 break;
813 }
814 }
815 else
816 {
817 Ok = FALSE;
818 break;
819 }
820
821
822
823 dump1[2]='\0';
824 dump2[2]='\0';
825 sscanf(t,"%c%c-%c%c %lf %lf",&dump1[0],&dump1[1],&dump2[0],&dump2[1],
826 &terms[n].c,
827 &terms[n].d);
828
829 if(strlen(t)>31)
830 terms[n].cc = g_strdup(t+31);
831 else
832 terms[n].cc = g_strdup(" ");
833
834 terms[n].cc[strlen(terms[n].cc)-1]='\0';
835
836 if(dump1[1]==' ') dump1[1]='\0';
837 if(dump2[1]==' ') dump2[1]='\0';
838
839 terms[n].numbers[0] = getNumberType(amberParameters,dump1);
840 terms[n].numbers[1] = getNumberType(amberParameters,dump2);
841
842 n++;
843 terms = g_realloc(terms,(n+1)*sizeof(AmberHydrogenBondedTerms));
844 }
845 if(n==0 || !Ok )
846 g_free(terms);
847 else
848 {
849 amberParameters->numberOfHydrogenBonded = n;
850 amberParameters->hydrogenBondedTerms = terms;
851 }
852 /* printing for test*/
853 /*
854 printf("number of hydrogen bonds terms = %d \n",amberParameters->numberOfHydrogenBonded);
855 for(n=0;n<amberParameters->numberOfHydrogenBonded;n++)
856 {
857 printf("%d %d %f %f\n",
858 amberParameters->hydrogenBondedTerms[n].numbers[0],
859 amberParameters->hydrogenBondedTerms[n].numbers[1],
860 amberParameters->hydrogenBondedTerms[n].c,
861 amberParameters->hydrogenBondedTerms[n].d
862 );
863 }
864 printf("\n");
865 */
866
867 return TRUE;
868
869
870 }
871 /**********************************************************************/
readAmberNonBondedTerms(AmberParameters * amberParameters,FILE * file)872 gboolean readAmberNonBondedTerms(AmberParameters* amberParameters,FILE* file)
873 {
874 gchar t[1024];
875 gchar dump1[1024];
876 gchar dump2[1024];
877 gint len = 1024;
878 gboolean Ok = FALSE;
879 gint n = 0;
880 AmberNonBondedTerms* terms = NULL;
881 gint nLigneEquiv = 0;
882 gint* nTypes = NULL;
883 gint **matrixNumbers = NULL;
884 gchar** strLines = NULL;
885 gint i;
886 gint j;
887
888
889 /* Search Equivalence */
890 Ok = FALSE;
891 while(!feof(file))
892 {
893 if(fgets(t,len,file))
894 {
895 if(strstr(t,nonBondedEquivTitle))
896 {
897 Ok = TRUE;
898 break;
899 }
900 }
901 }
902 if(Ok)
903 {
904 /*
905 */
906 strLines = g_malloc(sizeof(gint*));
907 nLigneEquiv = 0;
908 while(!feof(file))
909 {
910 if(fgets(t,len,file))
911 {
912 if(strstr(t,"End"))
913 {
914 Ok = TRUE;
915 break;
916 }
917
918 strLines[nLigneEquiv] = g_strdup(t);
919 nLigneEquiv++;
920 strLines = g_realloc(strLines,(nLigneEquiv+1)*sizeof(gint*));
921 }
922 else
923 {
924 Ok = FALSE;
925 break;
926 }
927 }
928 }
929 if(nLigneEquiv>0)
930 {
931 gchar name[10];
932 gchar* pos;
933
934 strLines = g_realloc(strLines,(nLigneEquiv)*sizeof(gint*));
935 nTypes = g_malloc(nLigneEquiv*sizeof(gint));
936 matrixNumbers = g_malloc(sizeof(gint*));
937 for(i=0;i<nLigneEquiv;i++)
938 {
939 nTypes[i] = (strlen(strLines[i])-2)/4+1;
940 matrixNumbers[i] = g_malloc(nTypes[i]*sizeof(gint));
941 pos = strLines[i];
942 for(j=0;j<nTypes[i];j++)
943 {
944 sscanf(pos,"%s",name);
945 matrixNumbers[i][j] = getNumberType(amberParameters,name);
946 pos += 4;
947
948 }
949
950 }
951 for(i=0;i<nLigneEquiv;i++)
952 if(strLines[i])
953 g_free(strLines[i]);
954
955 if(strLines)
956 g_free(strLines);
957
958 }
959 /* Search Begin INPUT FOR NON-BONDED */
960 Ok = FALSE;
961 while(!feof(file))
962 {
963 if(fgets(t,len,file))
964 {
965 if(strstr(t,nonBondedTitle))
966 {
967 Ok = TRUE;
968 break;
969 }
970 }
971 }
972 if(!Ok)
973 return FALSE;
974
975 terms = g_malloc(sizeof(AmberNonBondedTerms));
976 n = 0;
977 Ok = FALSE;
978 while(!feof(file))
979 {
980 if(fgets(t,len,file))
981 {
982 if(strstr(t,"End"))
983 {
984 Ok = TRUE;
985 break;
986 }
987 }
988 else
989 {
990 Ok = FALSE;
991 break;
992 }
993
994
995
996 dump1[2]='\0';
997 dump2[2]='\0';
998 sscanf(t,"%s %lf %lf",dump1,
999 &terms[n].r,
1000 &terms[n].epsilon);
1001
1002 if(strlen(t)>41)
1003 terms[n].c = g_strdup(t+41);
1004 else
1005 terms[n].c = g_strdup(" ");
1006
1007 terms[n].c[strlen(terms[n].c)-1]='\0';
1008
1009 if(dump1[1]==' ') dump1[1]='\0';
1010 terms[n].number = getNumberType(amberParameters,dump1);
1011 n++;
1012 terms = g_realloc(terms,(n+1)*sizeof(AmberNonBondedTerms));
1013 /* equivalence */
1014 if(nLigneEquiv>0)
1015 {
1016 gint numberOld = terms[n-1].number;
1017 gint numberLine = -1;
1018 for(i=0;i<nLigneEquiv;i++)
1019 {
1020 for(j=0;j<nTypes[i];j++)
1021 {
1022 if(numberOld == matrixNumbers[i][j])
1023 {
1024 numberLine = i;
1025 break;
1026 }
1027 }
1028 }
1029 if(numberLine != -1)
1030 {
1031 for(j=0;j<nTypes[numberLine];j++)
1032 {
1033 if(
1034 numberOld != matrixNumbers[numberLine][j] &&
1035 matrixNumbers[numberLine][j]>-1
1036 )
1037 {
1038 terms[n].number = matrixNumbers[numberLine][j];
1039 terms[n].r = terms[n-1].r;
1040 terms[n].epsilon = terms[n-1].epsilon;
1041 terms[n].c = g_strdup(terms[n-1].c);
1042 n++;
1043 terms =
1044 g_realloc(terms,(n+1)*sizeof(AmberNonBondedTerms));
1045 }
1046
1047 }
1048 }
1049 }
1050 }
1051 if(nLigneEquiv>0)
1052 {
1053 for(i=0;i<nLigneEquiv;i++)
1054 if(matrixNumbers[i])
1055 g_free(matrixNumbers[i]);
1056
1057 }
1058 if(matrixNumbers)
1059 g_free(matrixNumbers);
1060 if(nTypes)
1061 g_free(nTypes);
1062
1063 if(n==0 || !Ok )
1064 g_free(terms);
1065 else
1066 {
1067 amberParameters->numberOfNonBonded = n;
1068 amberParameters->nonBondedTerms = terms;
1069 }
1070 /* printing for test*/
1071 /*
1072 printf("number of non bended terms = %d \n",amberParameters->numberOfNonBonded);
1073 for(n=0;n<amberParameters->numberOfNonBonded;n++)
1074 {
1075 printf("%d %f %f\n",
1076 amberParameters->nonBondedTerms[n].number,
1077 amberParameters->nonBondedTerms[n].r,
1078 amberParameters->nonBondedTerms[n].epsilon
1079 );
1080 }
1081 printf("\n");
1082 */
1083
1084 return TRUE;
1085
1086
1087 }
1088 /**********************************************************************/
readAmberParameters(AmberParameters * amberParameters)1089 void readAmberParameters(AmberParameters* amberParameters)
1090 {
1091 FILE* file;
1092
1093 file = fopen("amber.prm","r");
1094 if(file == NULL)
1095 printf("file amber.prm not found\n");
1096 else
1097 {
1098 readAmberTypes(amberParameters,file);
1099 readAmberBondStretchTerms(amberParameters,file);
1100 readAmberAngleBendTerms(amberParameters,file);
1101 readAmberDihedralAngleTerms(amberParameters,file);
1102 readAmberImproperTorsionTerms(amberParameters,file);
1103 readAmberHydrogenBondedTerms(amberParameters,file);
1104 readAmberNonBondedTerms(amberParameters,file);
1105 fclose(file);
1106 }
1107 }
1108
1109 /**********************************************************************/
getSymbol(gchar * name)1110 gchar* getSymbol(gchar* name)
1111 {
1112 gchar* t;
1113 if(strlen(name)==1)
1114 return name;
1115 if(strcmp(name,"C0")==0)
1116 return "Ca";
1117 if(strcmp(name,"CU")==0)
1118 return "Cu";
1119 if(strcmp(name,"FE")==0)
1120 return "Fe";
1121 if(strcmp(name,"MG")==0)
1122 return "Mg";
1123 if(strcmp(name,"IM")==0)
1124 return "Cl";
1125 if(strcmp(name,"IB")==0)
1126 return "I";
1127
1128 if(strcmp(name,"Li")==0)
1129 return "Li";
1130 if(strcmp(name,"Na")==0)
1131 return "Na";
1132 if(strcmp(name,"IP")==0)
1133 return "Na";
1134 if(strcmp(name,"Cs")==0)
1135 return "Cs";
1136 if(strcmp(name,"Rb")==0)
1137 return "Rb";
1138 if(strcmp(name,"Zn")==0)
1139 return "Zn";
1140 if(strcmp(name,"Cl")==0)
1141 return "Cl";
1142 if(strcmp(name,"Br")==0)
1143 return "Br";
1144
1145 t = g_strdup_printf("%c",name[0]);
1146 return t;
1147 }
1148 /**********************************************************************/
printfMmCFileBond(FILE * file,AmberParameters * amberParameters)1149 void printfMmCFileBond(FILE* file,AmberParameters* amberParameters)
1150 {
1151 gint i;
1152 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",bondStretchTitle);
1153
1154 for(i=0;i<amberParameters->numberOfStretchTerms;i++)
1155 {
1156 fprintf(file,"\tfprintf(fout,\"%d\t%d\t%6.1f\t%7.4f\t\t%s\\n\");\n",
1157 amberParameters->bondStretchTerms[i].numbers[0]+1,
1158 amberParameters->bondStretchTerms[i].numbers[1]+1,
1159 amberParameters->bondStretchTerms[i].forceConstant,
1160 amberParameters->bondStretchTerms[i].equilibriumDistance,
1161 amberParameters->bondStretchTerms[i].c
1162 );
1163 }
1164 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1165 }
1166 /**********************************************************************/
printfMmCFileBend(FILE * file,AmberParameters * amberParameters)1167 void printfMmCFileBend(FILE* file,AmberParameters* amberParameters)
1168 {
1169 gint n;
1170 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",angleBendTitle);
1171
1172 for(n=0;n<amberParameters->numberOfBendTerms;n++)
1173 {
1174 fprintf(file,"\tfprintf(fout,\"%d\t%d\t%d\t%6.3f\t%6.2f\t\t%s\\n\");\n",
1175 amberParameters->angleBendTerms[n].numbers[0]+1,
1176 amberParameters->angleBendTerms[n].numbers[1]+1,
1177 amberParameters->angleBendTerms[n].numbers[2]+1,
1178 amberParameters->angleBendTerms[n].forceConstant,
1179 amberParameters->angleBendTerms[n].equilibriumAngle,
1180 amberParameters->angleBendTerms[n].c
1181 );
1182 }
1183
1184 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1185 }
1186 /**********************************************************************/
printfMmCFileHydrogenBonded(FILE * file,AmberParameters * amberParameters)1187 void printfMmCFileHydrogenBonded(FILE* file,AmberParameters* amberParameters)
1188 {
1189 gint n;
1190 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",hydrogenBondedTitle);
1191 for(n=0;n<amberParameters->numberOfHydrogenBonded;n++)
1192 {
1193 fprintf(file,"\tfprintf(fout,\"%d\t%d\t%8.1f\t%8.1f\t\t%s\\n\");\n",
1194 amberParameters->hydrogenBondedTerms[n].numbers[0]+1,
1195 amberParameters->hydrogenBondedTerms[n].numbers[1]+1,
1196 amberParameters->hydrogenBondedTerms[n].c,
1197 amberParameters->hydrogenBondedTerms[n].d,
1198 amberParameters->hydrogenBondedTerms[n].cc
1199 );
1200 }
1201
1202 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1203 }
1204 /**********************************************************************/
printfMmCFileimproperTorsion(FILE * file,AmberParameters * amberParameters)1205 void printfMmCFileimproperTorsion(FILE* file,AmberParameters* amberParameters)
1206 {
1207 gint n;
1208 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",improperTorsionTitle);
1209 for(n=0;n<amberParameters->numberOfImproperTorsionTerms;n++)
1210 {
1211 fprintf(file,"\tfprintf(fout,\"%d\t%d\t%d\t%d\t%6.3f\t%6.2f\t%4.1f\t\t%s\\n\");\n",
1212 amberParameters->improperTorsionTerms[n].numbers[0]+1,
1213 amberParameters->improperTorsionTerms[n].numbers[1]+1,
1214 amberParameters->improperTorsionTerms[n].numbers[2]+1,
1215 amberParameters->improperTorsionTerms[n].numbers[3]+1,
1216 amberParameters->improperTorsionTerms[n].barrier,
1217 amberParameters->improperTorsionTerms[n].phase,
1218 amberParameters->improperTorsionTerms[n].n,
1219 amberParameters->improperTorsionTerms[n].c
1220 );
1221 }
1222
1223 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1224 }
1225 /**********************************************************************/
printfMmCFileNonBonded(FILE * file,AmberParameters * amberParameters)1226 void printfMmCFileNonBonded(FILE* file,AmberParameters* amberParameters)
1227 {
1228 gint n;
1229 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",nonBondedTitle);
1230
1231 for(n=0;n<amberParameters->numberOfNonBonded;n++)
1232 {
1233 fprintf(file,"\tfprintf(fout,\"%d\t%8.4f\t%8.4f\t\t%s\\n\");\n",
1234 amberParameters->nonBondedTerms[n].number+1,
1235 amberParameters->nonBondedTerms[n].r,
1236 amberParameters->nonBondedTerms[n].epsilon,
1237 amberParameters->nonBondedTerms[n].c
1238 );
1239 }
1240
1241
1242
1243 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1244 }
1245
1246
1247 /**********************************************************************/
printfMmCFileDihedralAngle(FILE * file,AmberParameters * amberParameters)1248 void printfMmCFileDihedralAngle(FILE* file,AmberParameters* amberParameters)
1249 {
1250 gint n;
1251 gdouble nn;
1252
1253 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",dihedralAngleTitle);
1254 for(n=0;n<amberParameters->numberOfDihedralTerms;n++)
1255 {
1256 gint j;
1257 for(j=0;j<amberParameters->dihedralAngleTerms[n].nSomme;j++)
1258 {
1259 if(j==amberParameters->dihedralAngleTerms[n].nSomme-1)
1260 nn = amberParameters->dihedralAngleTerms[n].n[j];
1261 else
1262 nn = -amberParameters->dihedralAngleTerms[n].n[j];
1263
1264 fprintf(file,"\tfprintf(fout,\"%d\t%d\t%d\t%d\t%4.1f\t%6.3f\t%6.2f\t%4.1f\t\t%s\\n\");\n",
1265 amberParameters->dihedralAngleTerms[n].numbers[0]+1,
1266 amberParameters->dihedralAngleTerms[n].numbers[1]+1,
1267 amberParameters->dihedralAngleTerms[n].numbers[2]+1,
1268 amberParameters->dihedralAngleTerms[n].numbers[3]+1,
1269 amberParameters->dihedralAngleTerms[n].divisor[j],
1270 amberParameters->dihedralAngleTerms[n].barrier[j],
1271 amberParameters->dihedralAngleTerms[n].phase[j],
1272 nn,
1273 amberParameters->dihedralAngleTerms[n].c
1274
1275 );
1276 }
1277 }
1278 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1279 }
1280 /**********************************************************************/
printfMmCFileAtomTypes(FILE * file,AmberParameters * amberParameters)1281 void printfMmCFileAtomTypes(FILE* file,AmberParameters* amberParameters)
1282 {
1283 gint i;
1284 fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",atomTypesTitle);
1285
1286 for(i=0;i<amberParameters->numberOfTypes;i++)
1287 {
1288 fprintf(file,"\tfprintf(fout,\"%s\t%s\t%d\t%6.3f\t%6.3f\t\t%s\\n\");\n",
1289 amberParameters->atomTypes[i].name,
1290 getSymbol(amberParameters->atomTypes[i].name),
1291 amberParameters->atomTypes[i].number+1,
1292 amberParameters->atomTypes[i].masse,
1293 amberParameters->atomTypes[i].polarisability,
1294 amberParameters->atomTypes[i].c
1295 );
1296 }
1297
1298 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1299 }
1300
1301 /**********************************************************************/
printfMmCFileTitle(FILE * file)1302 void printfMmCFileTitle(FILE* file)
1303 {
1304 fprintf(file,"\tfprintf(fout,\"Begin Title\\n\");\n");
1305 fprintf(file,"\tfprintf(fout,\"\tAtom Types : Ty(Type) Symbol Numero Masse(C12 UMA) Polarisablities(Ang**3) \\n\");\n");
1306 fprintf(file,"\tfprintf(fout,\"\tBond Length : N1-N2 Force(Kcal/mol/A**2) Re\\n\");\n");
1307 fprintf(file,"\tfprintf(fout,\"\tBond Angle : N1-N2-N3 Force(Kcal/mol/rad**2) Angle(Deg) \\n\");\n");
1308 fprintf(file,"\tfprintf(fout,\"\tDihedral : N1-N2-N3-N4 Idiv Pk Phase(Deg) Pn \\n\");\n");
1309 fprintf(file,"\tfprintf(fout,\"\t E = Pk/Idiv*(1 + cos(P,*Phi - Phase)\\n\");\n");
1310 fprintf(file,"\tfprintf(fout,\"\t Pk = Barrier/2 Kcal/mol\\n\");\n");
1311 fprintf(file,"\tfprintf(fout,\"\t Idiv barrier is divised by Idiv\\n\");\n");
1312 fprintf(file,"\tfprintf(fout,\"\t Pn = periodicity fo the torional barrier\\n\");\n");
1313 fprintf(file,"\tfprintf(fout,\"\t if Pn<0 the tosional potential is \\n\");\n");
1314 fprintf(file,"\tfprintf(fout,\"\t assumed to have more than one term\\n\");\n");
1315 fprintf(file,"\tfprintf(fout,\"\t if Ni=0 => N is a number for any one Type\\n\");\n");
1316 fprintf(file,"\tfprintf(fout,\"\tImproper Dihedral : N1-N2-N3-N4 Idiv Pk Phase(Deg) Pn \\n\");\n");
1317 fprintf(file,"\tfprintf(fout,\"\tH-Bond : N1-N2 A(coef. 1/r**12) B(coef. -B/r**10)\\n\");\n");
1318 fprintf(file,"\tfprintf(fout,\"End\\n\");\n");
1319 }
1320 /**********************************************************************/
printfMmCFile(AmberParameters * amberParameters)1321 void printfMmCFile(AmberParameters* amberParameters)
1322 {
1323 FILE* file;
1324 file = fopen("createMMFile.fc","w");
1325 if(!file)
1326 {
1327 printf("I can not open createMMFile.fc\n");
1328 return;
1329 }
1330 fprintf(file,"gboolean createMMFile(gchar* filename)\n");
1331 fprintf(file,"{\n");
1332 fprintf(file,"\tFILE* fout = fopen(filename,\"w\");\n\n");
1333 fprintf(file,"\tif(fout==NULL)\n");
1334 fprintf(file,"\t{\n");
1335 fprintf(file,"\t\treturn FALSE;\n");
1336 fprintf(file,"\t}\n");
1337 // fprintf(file,"\tfprintf(fout,\"Natoms = %d\\n\");\n",natoms);
1338 printfMmCFileTitle(file);
1339 printfMmCFileAtomTypes(file,amberParameters);
1340 printfMmCFileBond(file,amberParameters);
1341 printfMmCFileBend(file,amberParameters);
1342 printfMmCFileDihedralAngle(file,amberParameters);
1343 printfMmCFileimproperTorsion(file,amberParameters);
1344 printfMmCFileHydrogenBonded(file,amberParameters);
1345 printfMmCFileNonBonded(file,amberParameters);
1346
1347 fprintf(file,"\tfclose(fout);\n");
1348 fprintf(file,"\treturn TRUE;\n");
1349 fprintf(file,"}\n");
1350 fclose(file);
1351 }
main()1352 int main()
1353 {
1354 AmberParameters amberParameters;
1355 readAmberParameters(&amberParameters);
1356 printfMmCFile(&amberParameters);
1357
1358 return 0;
1359
1360 }
1361