--- trunk/src/mengine/src/read_sdf.c 2008/07/31 16:34:52 48 +++ trunk/src/mengine/src/read_sdf.c 2008/12/01 06:44:59 58 @@ -3,9 +3,12 @@ #include "pcwin.h" #include "pcmod.h" -#include "substr.h" #include "atom_k.h" #include "energies.h" +#include "fix.h" + +#define TABULATOR_INCLUDE_IMPLEMENTATION +#include "tabulator.h" EXTERN struct t_files { int nfiles, append, batch, icurrent, ibatno; @@ -27,6 +30,8 @@ float etot,htot,stot,gtot,cptot; } vibdata; +// ============================ + int make_atom(int, float, float, float,char *); void make_bond(int, int, int); int read_sdf(int,int); @@ -35,12 +40,10 @@ void message_alert(char *, char *); void hdel(int); FILE * fopen_path ( char * , char * , char * ) ; -void read_schakal(int,int); -void write_schakal(void); -void mopaco(int,int); int FetchRecord(FILE *, char *); void avgleg(void); - +void get_tag(char *,char *); +// ================================== /* * flags - indicates whether dipole, xlogp or vibrational calculations were done * and if so write them to the output file. Look at the definitions in pcmod.h @@ -49,7 +52,7 @@ void write_sdf(int flags) { int i,j,lptest,nbond,junk; - FILE *wfile; + FILE *wfile = pcmoutfile; lptest = 0; for( i = 1; i <= natom; i++ ) @@ -81,7 +84,6 @@ else wfile = fopen_path(Savebox.path,Savebox.fname,"w"); */ - wfile = stdout; j = strlen(Struct_Title); for (i=0; i < j; i++) @@ -152,22 +154,27 @@ fprintf(wfile,"$$$$\n"); // fclose(wfile); } + /* =================================== */ // fast read - assume file is open and positioned // read structure and down to end $$$$ // and return int rd_sdf(FILE *rfile) { - int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom; - int ncount,junk,junk1,got_title,istereo,iz; + int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom; + int ncount,junk,junk1,got_title,istereo,iz,nvalue; int jji, jj1, jj2, jj3, jj4; + int n_row; int has_Aromatic, xPlus,yPlus,zPlus, planar; int icount,itemp[4]; int Ret_Val; float xtmp, ytmp, ztmp, dz; + float fconst,min,max; char c1[4],c2[4]; char atomchar[3]; char inputline[150]; + char tag[30]; + char ***tab; Ret_Val = TRUE; got_title = FALSE; @@ -230,7 +237,7 @@ atom[newatom].mmx_type = itype; if (newatom == -1) { - printf("Atom %d %s not recognized structure skipped\n",i,atomchar); + fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar); Ret_Val = FALSE; } } @@ -274,347 +281,182 @@ } } // read to end of structure +// parse any tags here +// agreed tags +// +// +// +// +// + while (FetchRecord(rfile,inputline)) { if (strncasecmp(inputline,"$$$$",4) == 0) { return Ret_Val; + } else if (inputline[0] == '>') + { + get_tag(inputline,tag); + if (strcasecmp(tag,"fixed_atoms") == 0) + { + tab = tabulator_new_from_file_using_header(rfile,"ATOM",0); + if (tab) + { + fprintf(pcmlogfile,"got table\n"); + n_row = tabulator_height(tab); + for (i=0; i < n_row; i++) + { + ia1 = atoi(tab[i][0]); + if (ia1 > 0 && ia1 <= natom) + { + fx_atom.katom_fix[fx_atom.natom_fix] = ia1; + fx_atom.natom_fix++; + } + } + } + } else if (strcasecmp(tag,"restrained_atoms") == 0) + { + tab = tabulator_new_from_file_using_header(rfile,"ATOM1 MAX F_CONST X Y Z",0); + if (tab) + { + n_row = tabulator_height(tab); + for (i=0; i < n_row; i++) + { + ia1 = atoi(tab[i][0]); + max = atof(tab[i][1]); + fconst = atof(tab[i][2]); + if (ia1 > 0 && ia1 <= natom) + { + restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1; + restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst; + restrain_atom.restrain_max[restrain_atom.natom_restrain] = max; + restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions + restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y; + restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z; + if (strlen(tab[i][3]) > 0) // x postion + { + xtmp = atof(tab[i][3]); + restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp; + } + if (strlen(tab[i][4]) > 0) // y postion + { + xtmp = atof(tab[i][4]); + restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp; + } + if (strlen(tab[i][5]) > 0) // y postion + { + xtmp = atof(tab[i][3]); + restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp; + } + restrain_atom.natom_restrain++; + } + } + } + } else if (strcasecmp(tag,"restrained_distances") == 0) + { + tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0); + if (tab) + { + n_row = tabulator_height(tab); + for (i=0; i < n_row; i++) + { + ia1 = atoi(tab[i][0]); + ia2 = atoi(tab[i][1]); + min = atof(tab[i][2]); + max = atof(tab[i][3]); + fconst = atof(tab[i][4]); + if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom)) + { + fx_dist.kdfix[fx_dist.ndfix][0] = ia1; + fx_dist.kdfix[fx_dist.ndfix][1] = ia2; + fx_dist.fdconst[fx_dist.ndfix] = fconst; + fx_dist.min_dist[fx_dist.ndfix] = min; + fx_dist.max_dist[fx_dist.ndfix] = max; + fx_dist.ndfix++; + } + } + } + } else if (strcasecmp(tag,"restrained_angles") == 0) + { + tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0); + if (tab) + { + n_row = tabulator_height(tab); + for (i=0; i < n_row; i++) + { + ia1 = atoi(tab[i][0]); + ia2 = atoi(tab[i][1]); + ia3 = atoi(tab[i][2]); + min = atof(tab[i][3]); + max = atof(tab[i][4]); + fconst = atof(tab[i][5]); + if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom)) + { + fx_angle.kafix[fx_angle.nafix][0] = ia1; + fx_angle.kafix[fx_angle.nafix][1] = ia2; + fx_angle.kafix[fx_angle.nafix][2] = ia3; + fx_angle.faconst[fx_angle.nafix] = fconst; + fx_angle.min_ang[fx_angle.nafix] = min; + fx_angle.max_ang[fx_angle.nafix] = max; + fx_angle.nafix++; + } + } + } + } else if (strcasecmp(tag,"restrained_dihedrals") == 0) + { + tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0); + if (tab) + { + n_row = tabulator_height(tab); + for (i=0; i < n_row; i++) + { + ia1 = atoi(tab[i][0]); + ia2 = atoi(tab[i][1]); + ia3 = atoi(tab[i][2]); + ia4 = atoi(tab[i][3]); + min = atof(tab[i][4]); + max = atof(tab[i][5]); + fconst = atof(tab[i][6]); + if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom)) + { + fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1; + fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2; + fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3; + fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4; + fx_torsion.ftconst[fx_torsion.ntfix] = fconst; + fx_torsion.min_tor[fx_torsion.ntfix] = min; + fx_torsion.max_tor[fx_torsion.ntfix] = max; + fx_torsion.ntfix++; + } + } + } + } } } return Ret_Val; } -/* =================================== */ -int read_sdf(int istruct, int isubred) +// ===================== +void get_tag(char *line,char *tag) { - int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom; - int ncount, ibotptr,junk,junk1,got_title,istereo,iz; - int jji, jj1, jj2, jj3, jj4; - int has_Aromatic, xPlus,yPlus,zPlus, planar; - int icount,itemp[4]; - float xtmp, ytmp, ztmp, dz; - char c1[4],c2[4]; - char atomchar[3]; - char inputline[150]; - FILE *infile; - - infile = fopen_path(Openbox.path,Openbox.fname,"r"); - - if (infile == NULL) - { - message_alert("Error opening Concord file","Concord Setup"); - return FALSE; + int i,j,icount; + icount = 0; + strcpy(tag,""); + for (i=0; i < strlen(line); i++) + { + if (line[i] == '<') + { + for (j=i+1; j < strlen(line); j++) + { + if (line[j] == '>') + { + tag[icount] = '\0'; + return; + } else + { + tag[icount] = line[j]; + icount++; + } + } + } } - - if (isubred == 0) - ibotptr = 0; - else - { - ibotptr = natom; - substr.istract[substr.icurstr] = TRUE; - } - - if (istruct > 1) - { - ncount = 0; - while (FetchRecord(infile,inputline)) - { - if (strncasecmp(inputline,"$$$$",4) == 0) - { - ncount++; - if (ncount+1 == istruct) - goto L_1; - } - } - } -// start file read here -L_1: - got_title = FALSE; - xPlus = yPlus = zPlus = FALSE; - FetchRecord(infile,inputline); - sscanf(inputline,"SDF %s",Struct_Title); - got_title = TRUE; - /* if (strlen(inputline) > 4) - { - iz = strlen(inputline); - if (iz > 60) iz = 59; - for (i=4; i < iz; i++) - Struct_Title[i-4] = inputline[i]; - Struct_Title[i] = '\0'; - got_title = TRUE; - } */ - FetchRecord(infile,inputline); // blank line - FetchRecord(infile,inputline); // blank line - - FetchRecord(infile,inputline); // natom and nbond - for (i=0; i <4; i++) - { - c1[i] = inputline[i]; - c2[i] = inputline[i+3]; - } - c1[3] = '\0';c2[3] = '\0'; - niatom = atoi(c1); - nibond = atoi(c2); - if (niatom == 0) - return FALSE; - - for (i=0; i < niatom; i++) - { - FetchRecord(infile,inputline); - sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1); - - itype = 0; - if (xtmp != 0.0) xPlus= TRUE; - if (ytmp != 0.0) yPlus= TRUE; - if (ztmp != 0.0) zPlus= TRUE; - - iz = strlen(atomchar); - if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix - { - if (atomchar[0] == 'N') itype = 8; - if (atomchar[0] == 'O') itype = 6; - - if (junk1 != 0) - { - if (junk1 == 3 && itype == 8) - itype = 41; // N+ - if (junk1 == 5 && itype == 6) - itype = 42; // O- - } - } - - newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar); - if (newatom == -1) - { - printf("Atom %d %s not recognized structure skipped\n",i,atomchar); - return FALSE; - } - if (isubred == 1) - { - atom[newatom].substr[0] = 0; - atom[newatom].substr[0] |= (1L << substr.icurstr); - } - } - has_Aromatic = FALSE; - - planar = TRUE; - if (xPlus && yPlus && zPlus) planar = FALSE; - - for (i=0; i < nibond; i++) - { - FetchRecord(infile,inputline); - for (j=0; j <4; j++) - { - c1[j] = inputline[j]; - c2[j] = inputline[j+3]; - } - c1[3] = '\0';c2[3] = '\0'; - ia1 = atoi(c1); - ia2 = atoi(c2); - istereo = 0; - sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo); - if (ibond >= 4) - { - ibond = 1; - has_Aromatic = TRUE; - } - make_bond(ia1+ibotptr, ia2+ibotptr, ibond); - if (istereo == 1 && planar) - { - atom[ia2+ibotptr].z += 0.7; - if (atom[ia2+ibotptr].atomnum == 1) - atom[ia2+ibotptr].type = 60; - } - if (istereo == 6 && planar) - { - atom[ia2+ibotptr].z -= 0.7; - if (atom[ia2+ibotptr].atomnum == 1) - atom[ia2+ibotptr].type = 60; - } - } - FetchRecord(infile,inputline); // M END line - FetchRecord(infile,inputline); - if (got_title == FALSE) - strcpy(Struct_Title,inputline); - fclose(infile); - ncount = strlen(Struct_Title); - for (i=0; i < ncount; i++) - { - if (Struct_Title[i] == '\n') - { - Struct_Title[i] = '\0'; - break; - } - } -// search for quaternary centers - check if flat - dz = 0.0; - for (i=1; i <= natom; i++) - { - dz += atom[i].z; - jji = jj1 = jj2 = jj3 = jj4 = 0; - if (atom[i].type != 0) - { - for (j=0; j < 6; j++) - { - if (atom[i].iat[j] != 0) - jji++; - } - } - if (jji == 4) - { - if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 && - atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0) // flat center - { - for (j=0; j < 4; j++) - { - itemp[j] = 0; - if (atom[atom[i].iat[0]].iat[j] != 0) jj1++; - if (atom[atom[i].iat[1]].iat[j] != 0) jj2++; - if (atom[atom[i].iat[2]].iat[j] != 0) jj3++; - if (atom[atom[i].iat[3]].iat[j] != 0) jj4++; - } - icount = 0; - if (jj1 == 1) - { - itemp[icount] = atom[i].iat[0]; - icount++; - } - if (jj2 == 1) - { - itemp[icount] = atom[i].iat[1]; - icount++; - } - if (jj3 == 1) - { - itemp[icount] = atom[i].iat[2]; - icount++; - } - if (jj3 == 1) - { - itemp[icount] = atom[i].iat[3]; - icount++; - } - if (icount >= 2) - { - atom[itemp[0]].z += 0.5; - atom[itemp[1]].z -= 0.5; - } - } - } - } - for (i=1; i < natom; i++) - { - for (j=i+1; j <= natom; j++) - { - xtmp = atom[i].x - atom[j].x; - ytmp = atom[i].y - atom[j].y; - ztmp = atom[i].z - atom[j].z; - if ( (xtmp + ytmp + ztmp) < 0.1) - { - atom[i].x += 0.1; - atom[i].y += 0.1; - atom[i].z += 0.1; - atom[j].x -= 0.1; - atom[j].y -= 0.1; - atom[j].z -= 0.1; - } - } - } - - if (has_Aromatic == TRUE) - mopaco(1,natom); - return TRUE; -} -// ================================================== -void read_schakal(int isel,int isubred) -{ - int i, itype, newatom, iz; - float xtmp, ytmp, ztmp; - char dummy[30]; - char atomchar[6]; - char inputline[150]; - FILE *infile; - - infile = fopen_path(Openbox.path,Openbox.fname,"r"); - - if (infile == NULL) - { - message_alert("Error opening Schakal file","Schakal Setup"); - return; - } - while ( FetchRecord(infile,inputline)) - { - sscanf(inputline,"%s",dummy); - if (strcmp(dummy,"TITLE") == 0) - { - sscanf(inputline,"%s %s",dummy,Struct_Title); - } else if (strcmp(dummy,"ATOM") == 0) - { - sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp); - // strip off numbers - iz = strlen(atomchar); - for (i=0; i < iz; i++) - { - if (isdigit(atomchar[i])) - { - atomchar[i] = '\0'; - break; - } - } - itype = 0; - newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar); - } else if (strcmp(dummy,"END") == 0) - { - break; - } - } - fclose(infile); - mopaco(1,natom); -// } -// ================================================== -void write_schakal() -{ - int i,j,lptest; - char atomchar[6]; - FILE *wfile; - - lptest = 0; - for( i = 1; i <= natom; i++ ) - { - if( atom[i].mmx_type == 20 ) - { - lptest = 1; - hdel( lptest ); - break; - } - } - - /* now write the schakal file */ - if (files.append) - wfile = fopen_path(Savebox.path,Savebox.fname,"a"); - else - wfile = fopen_path(Savebox.path,Savebox.fname,"w"); - - j = strlen(Struct_Title); - for (i=0; i < j; i++) - { - if (Struct_Title[i] == '\n') - { - Struct_Title[i] = '\0'; - break; - } - } - fprintf(wfile,"TITLE %s\n",Struct_Title); - fprintf(wfile,"CELL \n"); - for (i=1; i <= natom; i++) - { - sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i); - - fprintf(wfile,"ATOM %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x, - atom[i].y, atom[i].z); - } - fprintf(wfile,"END \n"); - fclose(wfile); -} - -