ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
(Generate patch)
# Line 1 | Line 1
1 + /* NOTICE: this source code file has been modified for use with FreeMOL */
2   #define EXTERN extern
3  
4   #include "pcwin.h"
5   #include "pcmod.h"
5 #include "substr.h"
6   #include "atom_k.h"
7   #include "energies.h"
8 + #include "fix.h"
9 +
10 + #define TABULATOR_INCLUDE_IMPLEMENTATION
11 + #include "tabulator.h"
12  
13   EXTERN struct t_files {
14          int nfiles, append, batch, icurrent, ibatno;
# Line 26 | Line 30
30          float etot,htot,stot,gtot,cptot;
31         } vibdata;
32  
33 + // ============================
34 +
35   int make_atom(int, float, float, float,char *);
36   void make_bond(int, int, int);
37   int read_sdf(int,int);
# Line 34 | Line 40
40   void message_alert(char *, char *);
41   void hdel(int);
42   FILE * fopen_path ( char * , char * , char * ) ;
37 void read_schakal(int,int);
38 void write_schakal(void);
39 void mopaco(int,int);
43   int FetchRecord(FILE *, char *);
44   void avgleg(void);
45 <
45 > void get_tag(char *,char *);
46 > // ==================================
47   /*
48   * flags - indicates whether dipole, xlogp or vibrational calculations were done
49   * and if so write them to the output file. Look at the definitions in pcmod.h
# Line 48 | Line 52
52   void write_sdf(int flags)
53   {
54      int i,j,lptest,nbond,junk;
55 <    FILE *wfile;
55 >    FILE *wfile = pcmoutfile;
56      
57      lptest = 0;
58      for( i = 1; i <= natom; i++ )
# Line 74 | Line 78
78         }
79       }
80      /*     now write the concord file */
81 + /*
82      if (files.append)
83          wfile = fopen_path(Savebox.path,Savebox.fname,"a");
84      else
85          wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86 + */
87  
88      j = strlen(Struct_Title);
89      for (i=0; i < j; i++)
# Line 146 | Line 152
152  
153      fprintf(wfile,"\n");
154      fprintf(wfile,"$$$$\n");
155 <    fclose(wfile);
155 > //    fclose(wfile);
156   }
157 +
158   /* ===================================  */
159   // fast read - assume file is open and positioned
160   // read structure and down to end $$$$
161   // and return
162   int rd_sdf(FILE *rfile)
163   {
164 <    int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
165 <    int ncount,junk,junk1,got_title,istereo,iz;
164 >    int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom;
165 >    int ncount,junk,junk1,got_title,istereo,iz,nvalue;
166      int jji, jj1, jj2, jj3, jj4;
167 +    int n_row;
168      int has_Aromatic, xPlus,yPlus,zPlus, planar;
169      int icount,itemp[4];
170      int Ret_Val;
171      float xtmp, ytmp, ztmp, dz;
172 +    float fconst,min,max;
173      char  c1[4],c2[4];
174      char  atomchar[3];
175      char  inputline[150];
176 +    char  tag[30];
177 +    char ***tab;
178  
179       Ret_Val = TRUE;
180       got_title = FALSE;
181       xPlus = yPlus = zPlus = FALSE;
182 <     FetchRecord(rfile,inputline);
183 <     sscanf(inputline,"SDF %s",Struct_Title);
182 >     if ( 0 == FetchRecord(rfile,inputline) )return -1;
183 > //   sscanf(inputline,"SDF %s",Struct_Title);
184 >     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
185       got_title = TRUE;
186       /*     if (strlen(inputline) > 4)
187       {
# Line 225 | Line 237
237            atom[newatom].mmx_type = itype;
238          if (newatom == -1)
239          {
240 <            printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
240 >          fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
241              Ret_Val = FALSE;
242          }
243       }
# Line 251 | Line 263
263          {
264              ibond = 1;
265              has_Aromatic = TRUE;
266 + // TJO
267 +            Ret_Val = FALSE;
268          }
269          make_bond(ia1, ia2, ibond);
270          if (istereo == 1 && planar)
# Line 267 | Line 281
281          }
282       }
283   // read to end of structure
284 + // parse any tags here
285 + // agreed tags
286 + //  <FIXED_ATOMS>
287 + //  <RESTRAINED_ATOMS>
288 + //  <RESTRAINED_DISTANCES>
289 + //  <RESTRAINED_ANGLES>
290 + //  <RESTRAINED_DIHEDRALS>
291 +
292       while (FetchRecord(rfile,inputline))
293       {
294           if (strncasecmp(inputline,"$$$$",4) == 0)
295           {
296              return Ret_Val;
297 +         } else if (inputline[0] == '>')
298 +         {
299 +           get_tag(inputline,tag);
300 +           if (strcasecmp(tag,"fixed_atoms") == 0)
301 +             {
302 +               tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
303 +               if (tab)
304 +                 {
305 +                   fprintf(pcmlogfile,"got table\n");
306 +                   n_row = tabulator_height(tab);
307 +                   for (i=0; i < n_row; i++)
308 +                     {
309 +                       ia1 = atoi(tab[i][0]);
310 +                       if (ia1 > 0 && ia1 <= natom)
311 +                         {
312 +                           fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
313 +                           fx_atom.natom_fix++;
314 +                         }
315 +                     }
316 +                 }
317 +             } else if (strcasecmp(tag,"restrained_atoms") == 0)
318 +             {
319 +               tab = tabulator_new_from_file_using_header(rfile,"ATOM MAX F_CONST X Y Z",0);
320 +               if (tab)
321 +                 {
322 +                   n_row = tabulator_height(tab);
323 +                   for (i=0; i < n_row; i++)
324 +                     {
325 +                       ia1 = atoi(tab[i][0]);
326 +                       max = atof(tab[i][1]);
327 +                       fconst = atof(tab[i][2]);
328 +                       if (ia1 > 0 && ia1 <= natom)
329 +                         {
330 +                           restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
331 +                           restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
332 +                           restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
333 +                           restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
334 +                           restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
335 +                           restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
336 +                           if (strlen(tab[i][3]) > 0)  // x postion
337 +                             {
338 +                               xtmp = atof(tab[i][3]);
339 +                               restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
340 +                             }
341 +                           if (strlen(tab[i][4]) > 0)  // y postion
342 +                             {
343 +                               xtmp = atof(tab[i][4]);
344 +                               restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
345 +                             }
346 +                           if (strlen(tab[i][5]) > 0)  // y postion
347 +                             {
348 +                               xtmp = atof(tab[i][5]);
349 +                               restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
350 +                             }
351 +                           restrain_atom.natom_restrain++;
352 +                         }
353 +                     }
354 +                 }
355 +             } else if (strcasecmp(tag,"restrained_distances") == 0)
356 +             {
357 +               tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
358 +               if (tab)
359 +                 {
360 +                   n_row = tabulator_height(tab);
361 +                   for (i=0; i < n_row; i++)
362 +                     {
363 +                       ia1 = atoi(tab[i][0]);
364 +                       ia2 = atoi(tab[i][1]);
365 +                       min = atof(tab[i][2]);
366 +                       max = atof(tab[i][3]);
367 +                       fconst = atof(tab[i][4]);
368 +                       if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
369 +                         {
370 +                           fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
371 +                           fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
372 +                           fx_dist.fdconst[fx_dist.ndfix] = fconst;
373 +                           fx_dist.min_dist[fx_dist.ndfix] = min;
374 +                           fx_dist.max_dist[fx_dist.ndfix] = max;
375 +                           fx_dist.ndfix++;
376 +                         }
377 +                     }
378 +                 }
379 +             } else if (strcasecmp(tag,"restrained_angles") == 0)
380 +             {
381 +               tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
382 +               if (tab)
383 +                 {
384 +                   n_row = tabulator_height(tab);
385 +                   for (i=0; i < n_row; i++)
386 +                     {
387 +                       ia1 = atoi(tab[i][0]);
388 +                       ia2 = atoi(tab[i][1]);
389 +                       ia3 = atoi(tab[i][2]);
390 +                       min = atof(tab[i][3]);
391 +                       max = atof(tab[i][4]);
392 +                       fconst = atof(tab[i][5]);
393 +                       if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
394 +                         {
395 +                           fx_angle.kafix[fx_angle.nafix][0] = ia1;
396 +                           fx_angle.kafix[fx_angle.nafix][1] = ia2;
397 +                           fx_angle.kafix[fx_angle.nafix][2] = ia3;
398 +                           fx_angle.faconst[fx_angle.nafix] = fconst;
399 +                           fx_angle.min_ang[fx_angle.nafix] = min;
400 +                           fx_angle.max_ang[fx_angle.nafix] = max;
401 +                           fx_angle.nafix++;
402 +                         }
403 +                     }
404 +                 }
405 +             } else if (strcasecmp(tag,"restrained_dihedrals") == 0)
406 +             {
407 +               tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
408 +               if (tab)
409 +                 {
410 +                   n_row = tabulator_height(tab);
411 +                   for (i=0; i < n_row; i++)
412 +                     {
413 +                       ia1 = atoi(tab[i][0]);
414 +                       ia2 = atoi(tab[i][1]);
415 +                       ia3 = atoi(tab[i][2]);
416 +                       ia4 = atoi(tab[i][3]);
417 +                       min = atof(tab[i][4]);
418 +                       max = atof(tab[i][5]);
419 +                       fconst = atof(tab[i][6]);
420 +                       if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
421 +                         {
422 +                           fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
423 +                           fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
424 +                           fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
425 +                           fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
426 +                           fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
427 +                           fx_torsion.min_tor[fx_torsion.ntfix] = min;
428 +                           fx_torsion.max_tor[fx_torsion.ntfix] = max;
429 +                           fx_torsion.ntfix++;
430 +                         }
431 +                     }
432 +                 }
433 +             }
434           }
435       }
436       return Ret_Val;
437   }
438 < /* ===================================  */
439 < int read_sdf(int istruct, int isubred)
438 > // =====================
439 > void get_tag(char *line,char *tag)
440   {
441 <    int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
442 <    int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
443 <    int jji, jj1, jj2, jj3, jj4;
444 <    int has_Aromatic, xPlus,yPlus,zPlus, planar;
445 <    int icount,itemp[4];
446 <    float xtmp, ytmp, ztmp, dz;
447 <    char  c1[4],c2[4];
448 <    char  atomchar[3];
449 <    char  inputline[150];
450 <    FILE   *infile;
451 <    
452 <    infile = fopen_path(Openbox.path,Openbox.fname,"r");
453 <
454 <    if (infile == NULL)
455 <    {
456 <        message_alert("Error opening Concord file","Concord Setup");
457 <        return FALSE;
458 <    }
459 <
460 <    if (isubred == 0)
302 <       ibotptr = 0;
303 <    else
304 <    {
305 <       ibotptr = natom;
306 <       substr.istract[substr.icurstr] = TRUE;
441 >  int i,j,icount;
442 >  icount = 0;
443 >  strcpy(tag,"");
444 >  for (i=0; i < strlen(line); i++)
445 >    {
446 >      if (line[i] == '<')
447 >        {
448 >          for (j=i+1; j < strlen(line); j++)
449 >            {
450 >              if (line[j] == '>')
451 >                {
452 >                  tag[icount] = '\0';
453 >                  return;
454 >                } else
455 >                {
456 >                  tag[icount] = line[j];
457 >                  icount++;
458 >                }
459 >            }
460 >        }
461      }
308
309    if (istruct > 1)
310    {
311        ncount = 0;
312        while (FetchRecord(infile,inputline))
313        {
314            if (strncasecmp(inputline,"$$$$",4) == 0)
315            {
316                ncount++;
317                if (ncount+1 == istruct)
318                    goto L_1;
319            }
320        }
321    }
322 // start file read here
323 L_1:
324     got_title = FALSE;
325     xPlus = yPlus = zPlus = FALSE;
326     FetchRecord(infile,inputline);
327     sscanf(inputline,"SDF %s",Struct_Title);
328     got_title = TRUE;
329     /*     if (strlen(inputline) > 4)
330     {
331         iz = strlen(inputline);
332         if (iz > 60) iz = 59;
333         for (i=4; i < iz; i++)
334           Struct_Title[i-4] = inputline[i];
335         Struct_Title[i] = '\0';
336        got_title = TRUE;
337        } */
338     FetchRecord(infile,inputline); // blank line
339     FetchRecord(infile,inputline); // blank line
340    
341     FetchRecord(infile,inputline); // natom and nbond
342     for (i=0; i <4; i++)
343     {
344        c1[i] = inputline[i];
345        c2[i] = inputline[i+3];
346     }
347     c1[3] = '\0';c2[3] = '\0';
348     niatom = atoi(c1);
349     nibond = atoi(c2);
350     if (niatom == 0)
351       return FALSE;
352
353     for (i=0; i < niatom; i++)
354     {
355        FetchRecord(infile,inputline);
356        sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
357
358        itype = 0;
359        if (xtmp != 0.0) xPlus= TRUE;
360        if (ytmp != 0.0) yPlus= TRUE;
361        if (ztmp != 0.0) zPlus= TRUE;
362
363        iz = strlen(atomchar);
364        if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1)  // nitro fix
365        {
366            if (atomchar[0] == 'N') itype = 8;
367            if (atomchar[0] == 'O') itype = 6;
368            
369            if (junk1 != 0)
370            {
371                if (junk1 == 3 && itype == 8)
372                  itype = 41;  // N+
373                if (junk1 == 5 && itype == 6)
374                  itype = 42;  // O-
375            }
376        }
377
378        newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
379        if (newatom == -1)
380        {
381            printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
382            return FALSE;
383        }
384        if (isubred == 1)
385        {
386            atom[newatom].substr[0] = 0;
387            atom[newatom].substr[0] |= (1L << substr.icurstr);
388        }
389     }
390     has_Aromatic = FALSE;
391    
392     planar = TRUE;
393     if (xPlus && yPlus && zPlus) planar = FALSE;
394    
395     for (i=0; i < nibond; i++)
396     {
397        FetchRecord(infile,inputline);
398        for (j=0; j <4; j++)
399        {
400           c1[j] = inputline[j];
401           c2[j] = inputline[j+3];
402         }  
403        c1[3] = '\0';c2[3] = '\0';
404        ia1 = atoi(c1);
405        ia2 = atoi(c2);
406        istereo = 0;
407        sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
408        if (ibond >= 4)
409        {
410            ibond = 1;
411            has_Aromatic = TRUE;
412        }
413        make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
414        if (istereo == 1 && planar)
415        {
416            atom[ia2+ibotptr].z += 0.7;
417            if (atom[ia2+ibotptr].atomnum == 1)
418              atom[ia2+ibotptr].type = 60;
419        }
420        if (istereo == 6 && planar)
421        {
422            atom[ia2+ibotptr].z -= 0.7;
423            if (atom[ia2+ibotptr].atomnum == 1)
424             atom[ia2+ibotptr].type = 60;
425        }
426     }
427     FetchRecord(infile,inputline);  // M END line
428     FetchRecord(infile,inputline);
429     if (got_title == FALSE)
430        strcpy(Struct_Title,inputline);
431     fclose(infile);
432     ncount = strlen(Struct_Title);
433     for (i=0; i < ncount; i++)
434     {
435         if (Struct_Title[i] == '\n')
436         {
437             Struct_Title[i] = '\0';
438             break;
439         }
440     }
441 //   search for quaternary centers - check if flat
442     dz = 0.0;
443     for (i=1; i <= natom; i++)
444     {
445         dz += atom[i].z;
446         jji = jj1 = jj2 = jj3 = jj4 = 0;
447         if (atom[i].type != 0)
448         {
449             for (j=0; j < 6; j++)
450             {
451                 if (atom[i].iat[j] != 0)
452                   jji++;
453             }
454         }
455         if (jji == 4)
456         {
457             if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
458                   atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0)  // flat center
459             {
460                 for (j=0; j < 4; j++)
461                 {
462                     itemp[j] = 0;
463                     if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
464                     if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
465                     if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
466                     if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
467                 }
468                 icount = 0;
469                 if (jj1 == 1)
470                 {
471                     itemp[icount] = atom[i].iat[0];
472                     icount++;
473                 }
474                 if (jj2 == 1)
475                 {
476                     itemp[icount] = atom[i].iat[1];
477                     icount++;
478                 }
479                 if (jj3 == 1)
480                 {
481                     itemp[icount] = atom[i].iat[2];
482                     icount++;
483                 }
484                 if (jj3 == 1)
485                 {
486                     itemp[icount] = atom[i].iat[3];
487                     icount++;
488                 }
489                 if (icount >= 2)
490                 {
491                     atom[itemp[0]].z += 0.5;
492                     atom[itemp[1]].z -= 0.5;
493                 }
494             }
495         }
496     }
497     for (i=1; i < natom; i++)
498     {
499         for (j=i+1; j <= natom; j++)
500         {
501             xtmp = atom[i].x - atom[j].x;
502             ytmp = atom[i].y - atom[j].y;
503             ztmp = atom[i].z - atom[j].z;
504             if ( (xtmp + ytmp + ztmp) < 0.1)
505             {
506                 atom[i].x += 0.1;
507                 atom[i].y += 0.1;
508                 atom[i].z += 0.1;
509                 atom[j].x -= 0.1;
510                 atom[j].y -= 0.1;
511                 atom[j].z -= 0.1;
512             }
513         }
514     }
515              
516     if (has_Aromatic == TRUE)
517       mopaco(1,natom);
518     return TRUE;
462   }
520 // ==================================================
521 void read_schakal(int isel,int isubred)
522 {
523    int i, itype,  newatom, iz;
524    float xtmp, ytmp, ztmp;
525    char  dummy[30];
526    char  atomchar[6];
527    char  inputline[150];
528    FILE   *infile;
529    
530    infile = fopen_path(Openbox.path,Openbox.fname,"r");
531
532    if (infile == NULL)
533    {
534        message_alert("Error opening Schakal file","Schakal Setup");
535        return;
536    }
537    while ( FetchRecord(infile,inputline))
538    {
539        sscanf(inputline,"%s",dummy);
540        if (strcmp(dummy,"TITLE") == 0)
541        {
542            sscanf(inputline,"%s %s",dummy,Struct_Title);
543        } else if (strcmp(dummy,"ATOM") == 0)
544        {
545            sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
546            // strip off numbers
547            iz = strlen(atomchar);
548            for (i=0; i < iz; i++)
549            {
550                if (isdigit(atomchar[i]))
551                {
552                    atomchar[i] = '\0';
553                    break;
554                }
555            }
556            itype = 0;
557            newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
558        } else if (strcmp(dummy,"END") == 0)
559        {
560            break;
561        }
562    }
563    fclose(infile);
564    mopaco(1,natom);
565 //
566 }
567 // ==================================================
568 void write_schakal()
569 {
570    int i,j,lptest;
571    char atomchar[6];
572    FILE *wfile;
573    
574    lptest = 0;
575    for( i = 1; i <= natom; i++ )
576    {
577        if( atom[i].mmx_type == 20 )
578        {
579            lptest = 1;
580            hdel( lptest );
581            break;
582        }
583    }
584
585    /*     now write the schakal file */
586    if (files.append)
587        wfile = fopen_path(Savebox.path,Savebox.fname,"a");
588    else
589        wfile = fopen_path(Savebox.path,Savebox.fname,"w");
590
591    j = strlen(Struct_Title);
592    for (i=0; i < j; i++)
593    {
594        if (Struct_Title[i] == '\n')
595        {
596            Struct_Title[i] = '\0';
597            break;
598        }
599    }
600    fprintf(wfile,"TITLE  %s\n",Struct_Title);
601    fprintf(wfile,"CELL \n");
602    for (i=1; i <= natom; i++)
603    {
604        sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
605            
606        fprintf(wfile,"ATOM  %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
607           atom[i].y, atom[i].z);
608    }    
609    fprintf(wfile,"END \n");
610    fclose(wfile);
611 }
612
613

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines