#include <stdlib.h> #include <stdio.h> #include <string.h> #include <glib.h> #include <math.h> typedef struct _AmberParameters AmberParameters; typedef struct _AmberAtomTypes AmberAtomTypes; typedef struct _AmberBondStretchTerms AmberBondStretchTerms; typedef struct _AmberAngleBendTerms AmberAngleBendTerms; typedef struct _AmberDihedralAngleTerms AmberDihedralAngleTerms; typedef struct _AmberImproperTorsionTerms AmberImproperTorsionTerms; typedef struct _AmberNonBondedTerms AmberNonBondedTerms; typedef struct _AmberHydrogenBondedTerms AmberHydrogenBondedTerms; /************************************/ struct _AmberAtomTypes { gchar* name; gint number; gdouble masse; gdouble polarisability; gchar* c; }; /************************************/ struct _AmberBondStretchTerms { gint numbers[2]; gdouble equilibriumDistance; gdouble forceConstant; gchar* c; }; /************************************/ struct _AmberAngleBendTerms { gint numbers[3]; gdouble equilibriumAngle; gdouble forceConstant; gchar* c; }; /************************************/ struct _AmberDihedralAngleTerms { gint numbers[4]; gint nSomme; gdouble* divisor; gdouble* barrier; gdouble* phase; gdouble* n; gchar* c; }; /************************************/ struct _AmberImproperTorsionTerms { gint numbers[4]; gdouble barrier; gdouble phase; gdouble n; gchar* c; }; /************************************/ struct _AmberNonBondedTerms { gint number; gdouble r; gdouble epsilon; gchar* c; }; /************************************/ struct _AmberHydrogenBondedTerms { gint numbers[2]; gdouble c; gdouble d; gchar* cc; }; /************************************/ struct _AmberParameters { gint numberOfTypes; AmberAtomTypes* atomTypes; gint numberOfStretchTerms; AmberBondStretchTerms* bondStretchTerms; gint numberOfBendTerms; AmberAngleBendTerms* angleBendTerms; gint numberOfDihedralTerms; AmberDihedralAngleTerms* dihedralAngleTerms; gint numberOfImproperTorsionTerms; AmberImproperTorsionTerms* improperTorsionTerms; gint numberOfNonBonded; AmberNonBondedTerms* nonBondedTerms; gint numberOfHydrogenBonded; AmberHydrogenBondedTerms* hydrogenBondedTerms; }; /************************************/ static gchar atomTypesTitle[] = "Begin INPUT FOR ATOM TYPES, MASSE AND POLARISABILITIES"; static gchar bondStretchTitle[] = "Begin INPUT FOR BOND LENGTH PARAMETERS"; static gchar angleBendTitle[] = "Begin INPUT FOR BOND ANGLE PARAMETERS"; static gchar hydrogenBondedTitle[] = "Begin INPUT FOR H-BOND 10-12 POTENTIAL PARAMETERS"; static gchar improperTorsionTitle[] ="Begin INPUT FOR IMPROPER DIHEDRAL PARAMETERS"; static gchar nonBondedTitle[] ="Begin INPUT FOR THE NON-BONDED 6-12 POTENTIAL PARAMETERS"; static gchar nonBondedEquivTitle[] ="Begin INPUT FOR EQUIVALENCING ATOM TYPES FOR THE NON-BONDED 6-12 POTENTIAL PARAMETERS"; static gchar dihedralAngleTitle[] = "Begin INPUT FOR DIHEDRAL PARAMETERS"; /**********************************************************************/ AmberParameters newAmberParameters() { AmberParameters amberParameters; amberParameters.numberOfTypes = 0; amberParameters.atomTypes = NULL; amberParameters.numberOfStretchTerms = 0; amberParameters.bondStretchTerms = NULL; amberParameters.numberOfBendTerms = 0; amberParameters.angleBendTerms = NULL; amberParameters.numberOfDihedralTerms = 0; amberParameters.dihedralAngleTerms = NULL; amberParameters.numberOfImproperTorsionTerms = 0; amberParameters.improperTorsionTerms = NULL; amberParameters.numberOfNonBonded = 0; amberParameters.nonBondedTerms = NULL; amberParameters.numberOfHydrogenBonded = 0; amberParameters.hydrogenBondedTerms = NULL; return amberParameters; } /**********************************************************************/ void freeAmberParameters(AmberParameters* amberParameters) { gint i; for(i=0;i<amberParameters->numberOfTypes;i++) if(amberParameters->atomTypes[i].name) g_free(amberParameters->atomTypes[i].name); amberParameters->numberOfTypes = 0; if(amberParameters->atomTypes ) g_free(amberParameters->atomTypes ); amberParameters->atomTypes = NULL; amberParameters->numberOfStretchTerms = 0; if(amberParameters->bondStretchTerms) g_free(amberParameters->bondStretchTerms); amberParameters->bondStretchTerms = NULL; amberParameters->numberOfBendTerms = 0; if(amberParameters->angleBendTerms) g_free(amberParameters->angleBendTerms); amberParameters->angleBendTerms = NULL; for(i=0;i<amberParameters->numberOfDihedralTerms;i++) { if(amberParameters->dihedralAngleTerms[i].divisor) g_free(amberParameters->dihedralAngleTerms[i].divisor); if(amberParameters->dihedralAngleTerms[i].barrier) g_free(amberParameters->dihedralAngleTerms[i].barrier); if(amberParameters->dihedralAngleTerms[i].phase) g_free(amberParameters->dihedralAngleTerms[i].phase); if(amberParameters->dihedralAngleTerms[i].n) g_free(amberParameters->dihedralAngleTerms[i].n); } amberParameters->numberOfDihedralTerms = 0; if(amberParameters->dihedralAngleTerms) g_free(amberParameters->dihedralAngleTerms); amberParameters->dihedralAngleTerms = NULL; amberParameters->numberOfImproperTorsionTerms = 0; if(amberParameters->improperTorsionTerms) g_free(amberParameters->improperTorsionTerms); amberParameters->improperTorsionTerms = NULL; amberParameters->numberOfNonBonded = 0; if(amberParameters->nonBondedTerms) g_free(amberParameters->nonBondedTerms); amberParameters->nonBondedTerms = NULL; amberParameters->numberOfHydrogenBonded = 0; if(amberParameters->hydrogenBondedTerms) g_free(amberParameters->hydrogenBondedTerms); amberParameters->hydrogenBondedTerms = NULL; } /**********************************************************************/ gint getNumberType(AmberParameters* amberParameters, gchar* type) { gint i; gint nTypes = amberParameters->numberOfTypes; AmberAtomTypes* types = amberParameters->atomTypes; gint len = strlen(type); if(strcmp(type,"X")==0) return -1; for(i=0;i<nTypes;i++) { if(len == (gint)strlen(types[i].name) && strstr(types[i].name,type)) return types[i].number; } return -2; } /**********************************************************************/ gboolean readAmberTypes(AmberParameters* amberParameters, FILE* file) { gchar t[1024]; gchar dump[1024]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberAtomTypes* types = NULL; /* Search Begin INPUT FOR ATOM TYPES */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,atomTypesTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; types = g_malloc(sizeof(AmberAtomTypes)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } sscanf(t,"%s %lf %lf",dump,&types[n].masse,&types[n].polarisability); types[n].name = g_strdup(dump); types[n].c = g_strdup(t+37); types[n].c[strlen(types[n].c)-1]='\0'; types[n].number = n; n++; types = g_realloc(types,(n+1)*sizeof(AmberAtomTypes)); } if(n==0 || !Ok ) g_free(types); else { amberParameters->numberOfTypes = n; amberParameters->atomTypes = types; } /* printing for test*/ /* printf("umber of types = %d \n",amberParameters->numberOfTypes); for(n=0;n<amberParameters->numberOfTypes;n++) { printf("%s\t %d\t", amberParameters->atomTypes[n].name, amberParameters->atomTypes[n].number ); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberBondStretchTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[1024]; gchar dump2[1024]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberBondStretchTerms* terms = NULL; /* Search Begin INPUT FOR ATOM TYPES */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,bondStretchTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberBondStretchTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } dump1[2]='\0'; dump2[2]='\0'; sscanf(t,"%c%c-%c%c %lf %lf",&dump1[0],&dump1[1],&dump2[0],&dump2[1], &terms[n].forceConstant, &terms[n].equilibriumDistance); if(strlen(t)>28) terms[n].c = g_strdup(t+28); else terms[n].c = g_strdup(" "); terms[n].c[strlen(terms[n].c)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; if(dump2[1]==' ') dump2[1]='\0'; terms[n].numbers[0] = getNumberType(amberParameters,dump1); terms[n].numbers[1] = getNumberType(amberParameters,dump2); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberBondStretchTerms)); } if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfStretchTerms = n; amberParameters->bondStretchTerms = terms; } /* printing for test*/ /* printf("number of bonds = %d \n",amberParameters->numberOfStretchTerms); for(n=0;n<amberParameters->numberOfStretchTerms;n++) { printf("%d %d %f %f\n", amberParameters->bondStretchTerms[n].numbers[0], amberParameters->bondStretchTerms[n].numbers[1], amberParameters->bondStretchTerms[n].forceConstant, amberParameters->bondStretchTerms[n].equilibriumDistance ); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberAngleBendTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[10]; gchar dump2[10]; gchar dump3[10]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberAngleBendTerms* terms = NULL; /* Search Begin INPUT FOR ATOM TYPES */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,angleBendTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberAngleBendTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } dump1[2]='\0'; dump2[2]='\0'; dump3[2]='\0'; sscanf(t,"%c%c-%c%c-%c%c %lf %lf", &dump1[0],&dump1[1],&dump2[0],&dump2[1],&dump3[0],&dump3[1], &terms[n].forceConstant, &terms[n].equilibriumAngle); if(strlen(t)>32) terms[n].c = g_strdup(t+32); else terms[n].c = g_strdup(" "); terms[n].c[strlen(terms[n].c)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; if(dump2[1]==' ') dump2[1]='\0'; if(dump3[1]==' ') dump3[1]='\0'; terms[n].numbers[0] = getNumberType(amberParameters,dump1); terms[n].numbers[1] = getNumberType(amberParameters,dump2); terms[n].numbers[2] = getNumberType(amberParameters,dump3); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberAngleBendTerms)); } if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfBendTerms = n; amberParameters->angleBendTerms = terms; } /* printing for test*/ /* printf("number of bonds = %d \n",amberParameters->numberOfBendTerms); for(n=0;n<amberParameters->numberOfBendTerms;n++) { printf("%d %d %d %f %f\n", amberParameters->angleBendTerms[n].numbers[0], amberParameters->angleBendTerms[n].numbers[1], amberParameters->angleBendTerms[n].numbers[2], amberParameters->angleBendTerms[n].forceConstant, amberParameters->angleBendTerms[n].equilibriumAngle ); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberDihedralAngleTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[10]; gchar dump2[10]; gchar dump3[10]; gchar dump4[10]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberDihedralAngleTerms* terms = NULL; gdouble divisor = 1; gdouble barrier = 0; gdouble phase = 0; gdouble nN = 0; /* Search Begin INPUT FOR DIHEDRAL PARAMETERS */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,dihedralAngleTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberDihedralAngleTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } terms[n].nSomme = 1; terms[n].divisor = g_malloc(sizeof(gdouble)); terms[n].barrier = g_malloc(sizeof(gdouble)); terms[n].phase = g_malloc(sizeof(gdouble)); terms[n].n = g_malloc(sizeof(gdouble)); dump1[2]='\0'; dump2[2]='\0'; dump3[2]='\0'; dump4[2]='\0'; sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf %lf", &dump1[0],&dump1[1], &dump2[0],&dump2[1], &dump3[0],&dump3[1], &dump4[0],&dump4[1], &divisor, &barrier, &phase, &nN); terms[n].divisor[0] = divisor; terms[n].barrier[0] = barrier; terms[n].phase[0] = phase; terms[n].n[0] = fabs(nN); if(strlen(t)>60) terms[n].c = g_strdup(t+60); else terms[n].c = g_strdup(" "); terms[n].c[strlen(terms[n].c)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; if(dump2[1]==' ') dump2[1]='\0'; if(dump3[1]==' ') dump3[1]='\0'; if(dump4[1]==' ') dump4[1]='\0'; terms[n].numbers[0] = getNumberType(amberParameters,dump1); terms[n].numbers[1] = getNumberType(amberParameters,dump2); terms[n].numbers[2] = getNumberType(amberParameters,dump3); terms[n].numbers[3] = getNumberType(amberParameters,dump4); Ok = TRUE; while(!feof(file) && nN<0) { if(!fgets(t,len,file)) { Ok = FALSE; break; } terms[n].nSomme++; terms[n].divisor = g_realloc(terms[n].divisor,terms[n].nSomme*sizeof(gdouble)); terms[n].barrier = g_realloc(terms[n].barrier,terms[n].nSomme*sizeof(gdouble)); terms[n].phase = g_realloc(terms[n].phase,terms[n].nSomme*sizeof(gdouble)); terms[n].n = g_realloc(terms[n].n,terms[n].nSomme*sizeof(gdouble)); sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf %lf", &dump1[0],&dump1[1], &dump2[0],&dump2[1], &dump3[0],&dump3[1], &dump4[0],&dump4[1], &divisor, &barrier, &phase, &nN); terms[n].divisor[terms[n].nSomme-1] = divisor; terms[n].barrier[terms[n].nSomme-1] = barrier; terms[n].phase[terms[n].nSomme-1] = phase; terms[n].n[terms[n].nSomme-1] = fabs(nN); } if(!Ok) break; n++; terms = g_realloc(terms,(n+1)*sizeof(AmberDihedralAngleTerms)); } if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfDihedralTerms = n; amberParameters->dihedralAngleTerms = terms; } /* printing for test*/ /* printf("number of dihedral torsion terms = %d \n", amberParameters->numberOfDihedralTerms); for(n=0;n<amberParameters->numberOfDihedralTerms;n++) { gint j; printf("%d %d %d %d \t", amberParameters->dihedralAngleTerms[n].numbers[0], amberParameters->dihedralAngleTerms[n].numbers[1], amberParameters->dihedralAngleTerms[n].numbers[2], amberParameters->dihedralAngleTerms[n].numbers[3] ); for(j=0;j<amberParameters->dihedralAngleTerms[n].nSomme;j++) { printf("%f %f %f %f\t", amberParameters->dihedralAngleTerms[n].divisor[j], amberParameters->dihedralAngleTerms[n].barrier[j], amberParameters->dihedralAngleTerms[n].phase[j], amberParameters->dihedralAngleTerms[n].n[j] ); } printf("\n"); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberImproperTorsionTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[10]; gchar dump2[10]; gchar dump3[10]; gchar dump4[10]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberImproperTorsionTerms* terms = NULL; /* Search Begin INPUT FOR ATOM TYPES */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,improperTorsionTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberImproperTorsionTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } dump1[2]='\0'; dump2[2]='\0'; dump3[2]='\0'; dump4[2]='\0'; sscanf(t,"%c%c-%c%c-%c%c-%c%c %lf %lf %lf", &dump1[0],&dump1[1], &dump2[0],&dump2[1], &dump3[0],&dump3[1], &dump4[0],&dump4[1], &terms[n].barrier, &terms[n].phase, &terms[n].n); if(strlen(t)>60) terms[n].c = g_strdup(t+60); else terms[n].c = g_strdup(" "); terms[n].c[strlen(terms[n].c)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; if(dump2[1]==' ') dump2[1]='\0'; if(dump3[1]==' ') dump3[1]='\0'; if(dump4[1]==' ') dump4[1]='\0'; terms[n].numbers[0] = getNumberType(amberParameters,dump1); terms[n].numbers[1] = getNumberType(amberParameters,dump2); terms[n].numbers[2] = getNumberType(amberParameters,dump3); terms[n].numbers[3] = getNumberType(amberParameters,dump4); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberImproperTorsionTerms)); } if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfImproperTorsionTerms = n; amberParameters->improperTorsionTerms = terms; } /* printing for test*/ /* printf("number of improper torsion terms = %d \n", amberParameters->numberOfImproperTorsionTerms); for(n=0;n<amberParameters->numberOfImproperTorsionTerms;n++) { printf("%d %d %d %d %f %f %f\n", amberParameters->improperTorsionTerms[n].numbers[0], amberParameters->improperTorsionTerms[n].numbers[1], amberParameters->improperTorsionTerms[n].numbers[2], amberParameters->improperTorsionTerms[n].numbers[3], amberParameters->improperTorsionTerms[n].barrier, amberParameters->improperTorsionTerms[n].phase, amberParameters->improperTorsionTerms[n].n ); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberHydrogenBondedTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[10]; gchar dump2[10]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberHydrogenBondedTerms* terms = NULL; /* Search Begin INPUT FOR ATOM TYPES */ while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,hydrogenBondedTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberHydrogenBondedTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } dump1[2]='\0'; dump2[2]='\0'; sscanf(t,"%c%c-%c%c %lf %lf",&dump1[0],&dump1[1],&dump2[0],&dump2[1], &terms[n].c, &terms[n].d); if(strlen(t)>31) terms[n].cc = g_strdup(t+31); else terms[n].cc = g_strdup(" "); terms[n].cc[strlen(terms[n].cc)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; if(dump2[1]==' ') dump2[1]='\0'; terms[n].numbers[0] = getNumberType(amberParameters,dump1); terms[n].numbers[1] = getNumberType(amberParameters,dump2); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberHydrogenBondedTerms)); } if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfHydrogenBonded = n; amberParameters->hydrogenBondedTerms = terms; } /* printing for test*/ /* printf("number of hydrogen bonds terms = %d \n",amberParameters->numberOfHydrogenBonded); for(n=0;n<amberParameters->numberOfHydrogenBonded;n++) { printf("%d %d %f %f\n", amberParameters->hydrogenBondedTerms[n].numbers[0], amberParameters->hydrogenBondedTerms[n].numbers[1], amberParameters->hydrogenBondedTerms[n].c, amberParameters->hydrogenBondedTerms[n].d ); } printf("\n"); */ return TRUE; } /**********************************************************************/ gboolean readAmberNonBondedTerms(AmberParameters* amberParameters,FILE* file) { gchar t[1024]; gchar dump1[1024]; gchar dump2[1024]; gint len = 1024; gboolean Ok = FALSE; gint n = 0; AmberNonBondedTerms* terms = NULL; gint nLigneEquiv = 0; gint* nTypes = NULL; gint **matrixNumbers = NULL; gchar** strLines = NULL; gint i; gint j; /* Search Equivalence */ Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,nonBondedEquivTitle)) { Ok = TRUE; break; } } } if(Ok) { /* */ strLines = g_malloc(sizeof(gint*)); nLigneEquiv = 0; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } strLines[nLigneEquiv] = g_strdup(t); nLigneEquiv++; strLines = g_realloc(strLines,(nLigneEquiv+1)*sizeof(gint*)); } else { Ok = FALSE; break; } } } if(nLigneEquiv>0) { gchar name[10]; gchar* pos; strLines = g_realloc(strLines,(nLigneEquiv)*sizeof(gint*)); nTypes = g_malloc(nLigneEquiv*sizeof(gint)); matrixNumbers = g_malloc(sizeof(gint*)); for(i=0;i<nLigneEquiv;i++) { nTypes[i] = (strlen(strLines[i])-2)/4+1; matrixNumbers[i] = g_malloc(nTypes[i]*sizeof(gint)); pos = strLines[i]; for(j=0;j<nTypes[i];j++) { sscanf(pos,"%s",name); matrixNumbers[i][j] = getNumberType(amberParameters,name); pos += 4; } } for(i=0;i<nLigneEquiv;i++) if(strLines[i]) g_free(strLines[i]); if(strLines) g_free(strLines); } /* Search Begin INPUT FOR NON-BONDED */ Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,nonBondedTitle)) { Ok = TRUE; break; } } } if(!Ok) return FALSE; terms = g_malloc(sizeof(AmberNonBondedTerms)); n = 0; Ok = FALSE; while(!feof(file)) { if(fgets(t,len,file)) { if(strstr(t,"End")) { Ok = TRUE; break; } } else { Ok = FALSE; break; } dump1[2]='\0'; dump2[2]='\0'; sscanf(t,"%s %lf %lf",dump1, &terms[n].r, &terms[n].epsilon); if(strlen(t)>41) terms[n].c = g_strdup(t+41); else terms[n].c = g_strdup(" "); terms[n].c[strlen(terms[n].c)-1]='\0'; if(dump1[1]==' ') dump1[1]='\0'; terms[n].number = getNumberType(amberParameters,dump1); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberNonBondedTerms)); /* equivalence */ if(nLigneEquiv>0) { gint numberOld = terms[n-1].number; gint numberLine = -1; for(i=0;i<nLigneEquiv;i++) { for(j=0;j<nTypes[i];j++) { if(numberOld == matrixNumbers[i][j]) { numberLine = i; break; } } } if(numberLine != -1) { for(j=0;j<nTypes[numberLine];j++) { if( numberOld != matrixNumbers[numberLine][j] && matrixNumbers[numberLine][j]>-1 ) { terms[n].number = matrixNumbers[numberLine][j]; terms[n].r = terms[n-1].r; terms[n].epsilon = terms[n-1].epsilon; terms[n].c = g_strdup(terms[n-1].c); n++; terms = g_realloc(terms,(n+1)*sizeof(AmberNonBondedTerms)); } } } } } if(nLigneEquiv>0) { for(i=0;i<nLigneEquiv;i++) if(matrixNumbers[i]) g_free(matrixNumbers[i]); } if(matrixNumbers) g_free(matrixNumbers); if(nTypes) g_free(nTypes); if(n==0 || !Ok ) g_free(terms); else { amberParameters->numberOfNonBonded = n; amberParameters->nonBondedTerms = terms; } /* printing for test*/ /* printf("number of non bended terms = %d \n",amberParameters->numberOfNonBonded); for(n=0;n<amberParameters->numberOfNonBonded;n++) { printf("%d %f %f\n", amberParameters->nonBondedTerms[n].number, amberParameters->nonBondedTerms[n].r, amberParameters->nonBondedTerms[n].epsilon ); } printf("\n"); */ return TRUE; } /**********************************************************************/ void readAmberParameters(AmberParameters* amberParameters) { FILE* file; file = fopen("amber.prm","r"); if(file == NULL) printf("file amber.prm not found\n"); else { readAmberTypes(amberParameters,file); readAmberBondStretchTerms(amberParameters,file); readAmberAngleBendTerms(amberParameters,file); readAmberDihedralAngleTerms(amberParameters,file); readAmberImproperTorsionTerms(amberParameters,file); readAmberHydrogenBondedTerms(amberParameters,file); readAmberNonBondedTerms(amberParameters,file); fclose(file); } } /**********************************************************************/ gchar* getSymbol(gchar* name) { gchar* t; if(strlen(name)==1) return name; if(strcmp(name,"C0")==0) return "Ca"; if(strcmp(name,"CU")==0) return "Cu"; if(strcmp(name,"FE")==0) return "Fe"; if(strcmp(name,"MG")==0) return "Mg"; if(strcmp(name,"IM")==0) return "Cl"; if(strcmp(name,"IB")==0) return "I"; if(strcmp(name,"Li")==0) return "Li"; if(strcmp(name,"Na")==0) return "Na"; if(strcmp(name,"IP")==0) return "Na"; if(strcmp(name,"Cs")==0) return "Cs"; if(strcmp(name,"Rb")==0) return "Rb"; if(strcmp(name,"Zn")==0) return "Zn"; if(strcmp(name,"Cl")==0) return "Cl"; if(strcmp(name,"Br")==0) return "Br"; t = g_strdup_printf("%c",name[0]); return t; } /**********************************************************************/ void printfMmCFileBond(FILE* file,AmberParameters* amberParameters) { gint i; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",bondStretchTitle); for(i=0;i<amberParameters->numberOfStretchTerms;i++) { fprintf(file,"\tfprintf(fout,\"%d\t%d\t%6.1f\t%7.4f\t\t%s\\n\");\n", amberParameters->bondStretchTerms[i].numbers[0]+1, amberParameters->bondStretchTerms[i].numbers[1]+1, amberParameters->bondStretchTerms[i].forceConstant, amberParameters->bondStretchTerms[i].equilibriumDistance, amberParameters->bondStretchTerms[i].c ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileBend(FILE* file,AmberParameters* amberParameters) { gint n; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",angleBendTitle); for(n=0;n<amberParameters->numberOfBendTerms;n++) { fprintf(file,"\tfprintf(fout,\"%d\t%d\t%d\t%6.3f\t%6.2f\t\t%s\\n\");\n", amberParameters->angleBendTerms[n].numbers[0]+1, amberParameters->angleBendTerms[n].numbers[1]+1, amberParameters->angleBendTerms[n].numbers[2]+1, amberParameters->angleBendTerms[n].forceConstant, amberParameters->angleBendTerms[n].equilibriumAngle, amberParameters->angleBendTerms[n].c ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileHydrogenBonded(FILE* file,AmberParameters* amberParameters) { gint n; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",hydrogenBondedTitle); for(n=0;n<amberParameters->numberOfHydrogenBonded;n++) { fprintf(file,"\tfprintf(fout,\"%d\t%d\t%8.1f\t%8.1f\t\t%s\\n\");\n", amberParameters->hydrogenBondedTerms[n].numbers[0]+1, amberParameters->hydrogenBondedTerms[n].numbers[1]+1, amberParameters->hydrogenBondedTerms[n].c, amberParameters->hydrogenBondedTerms[n].d, amberParameters->hydrogenBondedTerms[n].cc ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileimproperTorsion(FILE* file,AmberParameters* amberParameters) { gint n; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",improperTorsionTitle); for(n=0;n<amberParameters->numberOfImproperTorsionTerms;n++) { 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", amberParameters->improperTorsionTerms[n].numbers[0]+1, amberParameters->improperTorsionTerms[n].numbers[1]+1, amberParameters->improperTorsionTerms[n].numbers[2]+1, amberParameters->improperTorsionTerms[n].numbers[3]+1, amberParameters->improperTorsionTerms[n].barrier, amberParameters->improperTorsionTerms[n].phase, amberParameters->improperTorsionTerms[n].n, amberParameters->improperTorsionTerms[n].c ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileNonBonded(FILE* file,AmberParameters* amberParameters) { gint n; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",nonBondedTitle); for(n=0;n<amberParameters->numberOfNonBonded;n++) { fprintf(file,"\tfprintf(fout,\"%d\t%8.4f\t%8.4f\t\t%s\\n\");\n", amberParameters->nonBondedTerms[n].number+1, amberParameters->nonBondedTerms[n].r, amberParameters->nonBondedTerms[n].epsilon, amberParameters->nonBondedTerms[n].c ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileDihedralAngle(FILE* file,AmberParameters* amberParameters) { gint n; gdouble nn; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",dihedralAngleTitle); for(n=0;n<amberParameters->numberOfDihedralTerms;n++) { gint j; for(j=0;j<amberParameters->dihedralAngleTerms[n].nSomme;j++) { if(j==amberParameters->dihedralAngleTerms[n].nSomme-1) nn = amberParameters->dihedralAngleTerms[n].n[j]; else nn = -amberParameters->dihedralAngleTerms[n].n[j]; 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", amberParameters->dihedralAngleTerms[n].numbers[0]+1, amberParameters->dihedralAngleTerms[n].numbers[1]+1, amberParameters->dihedralAngleTerms[n].numbers[2]+1, amberParameters->dihedralAngleTerms[n].numbers[3]+1, amberParameters->dihedralAngleTerms[n].divisor[j], amberParameters->dihedralAngleTerms[n].barrier[j], amberParameters->dihedralAngleTerms[n].phase[j], nn, amberParameters->dihedralAngleTerms[n].c ); } } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileAtomTypes(FILE* file,AmberParameters* amberParameters) { gint i; fprintf(file,"\tfprintf(fout,\"%s\\n\");\n",atomTypesTitle); for(i=0;i<amberParameters->numberOfTypes;i++) { fprintf(file,"\tfprintf(fout,\"%s\t%s\t%d\t%6.3f\t%6.3f\t\t%s\\n\");\n", amberParameters->atomTypes[i].name, getSymbol(amberParameters->atomTypes[i].name), amberParameters->atomTypes[i].number+1, amberParameters->atomTypes[i].masse, amberParameters->atomTypes[i].polarisability, amberParameters->atomTypes[i].c ); } fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFileTitle(FILE* file) { fprintf(file,"\tfprintf(fout,\"Begin Title\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tAtom Types : Ty(Type) Symbol Numero Masse(C12 UMA) Polarisablities(Ang**3) \\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tBond Length : N1-N2 Force(Kcal/mol/A**2) Re\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tBond Angle : N1-N2-N3 Force(Kcal/mol/rad**2) Angle(Deg) \\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tDihedral : N1-N2-N3-N4 Idiv Pk Phase(Deg) Pn \\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t E = Pk/Idiv*(1 + cos(P,*Phi - Phase)\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t Pk = Barrier/2 Kcal/mol\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t Idiv barrier is divised by Idiv\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t Pn = periodicity fo the torional barrier\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t if Pn<0 the tosional potential is \\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t assumed to have more than one term\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\t if Ni=0 => N is a number for any one Type\\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tImproper Dihedral : N1-N2-N3-N4 Idiv Pk Phase(Deg) Pn \\n\");\n"); fprintf(file,"\tfprintf(fout,\"\tH-Bond : N1-N2 A(coef. 1/r**12) B(coef. -B/r**10)\\n\");\n"); fprintf(file,"\tfprintf(fout,\"End\\n\");\n"); } /**********************************************************************/ void printfMmCFile(AmberParameters* amberParameters) { FILE* file; file = fopen("createMMFile.fc","w"); if(!file) { printf("I can not open createMMFile.fc\n"); return; } fprintf(file,"gboolean createMMFile(gchar* filename)\n"); fprintf(file,"{\n"); fprintf(file,"\tFILE* fout = fopen(filename,\"w\");\n\n"); fprintf(file,"\tif(fout==NULL)\n"); fprintf(file,"\t{\n"); fprintf(file,"\t\treturn FALSE;\n"); fprintf(file,"\t}\n"); // fprintf(file,"\tfprintf(fout,\"Natoms = %d\\n\");\n",natoms); printfMmCFileTitle(file); printfMmCFileAtomTypes(file,amberParameters); printfMmCFileBond(file,amberParameters); printfMmCFileBend(file,amberParameters); printfMmCFileDihedralAngle(file,amberParameters); printfMmCFileimproperTorsion(file,amberParameters); printfMmCFileHydrogenBonded(file,amberParameters); printfMmCFileNonBonded(file,amberParameters); fprintf(file,"\tfclose(fout);\n"); fprintf(file,"\treturn TRUE;\n"); fprintf(file,"}\n"); fclose(file); } int main() { AmberParameters amberParameters; readAmberParameters(&amberParameters); printfMmCFile(&amberParameters); return 0; }