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 80 | Line 84
84      else
85          wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86   */
83        wfile = stdout;
87  
88      j = strlen(Struct_Title);
89      for (i=0; i < j; i++)
# Line 151 | Line 154
154      fprintf(wfile,"$$$$\n");
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;
# Line 229 | 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 273 | 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,"ATOM1 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][3]);
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;
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      }
306
307    if (isubred == 0)
308       ibotptr = 0;
309    else
310    {
311       ibotptr = natom;
312       substr.istract[substr.icurstr] = TRUE;
313    }
314
315    if (istruct > 1)
316    {
317        ncount = 0;
318        while (FetchRecord(infile,inputline))
319        {
320            if (strncasecmp(inputline,"$$$$",4) == 0)
321            {
322                ncount++;
323                if (ncount+1 == istruct)
324                    goto L_1;
325            }
326        }
327    }
328 // start file read here
329 L_1:
330     got_title = FALSE;
331     xPlus = yPlus = zPlus = FALSE;
332     FetchRecord(infile,inputline);
333     sscanf(inputline,"SDF %s",Struct_Title);
334     got_title = TRUE;
335     /*     if (strlen(inputline) > 4)
336     {
337         iz = strlen(inputline);
338         if (iz > 60) iz = 59;
339         for (i=4; i < iz; i++)
340           Struct_Title[i-4] = inputline[i];
341         Struct_Title[i] = '\0';
342        got_title = TRUE;
343        } */
344     FetchRecord(infile,inputline); // blank line
345     FetchRecord(infile,inputline); // blank line
346    
347     FetchRecord(infile,inputline); // natom and nbond
348     for (i=0; i <4; i++)
349     {
350        c1[i] = inputline[i];
351        c2[i] = inputline[i+3];
352     }
353     c1[3] = '\0';c2[3] = '\0';
354     niatom = atoi(c1);
355     nibond = atoi(c2);
356     if (niatom == 0)
357       return FALSE;
358
359     for (i=0; i < niatom; i++)
360     {
361        FetchRecord(infile,inputline);
362        sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
363
364        itype = 0;
365        if (xtmp != 0.0) xPlus= TRUE;
366        if (ytmp != 0.0) yPlus= TRUE;
367        if (ztmp != 0.0) zPlus= TRUE;
368
369        iz = strlen(atomchar);
370        if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1)  // nitro fix
371        {
372            if (atomchar[0] == 'N') itype = 8;
373            if (atomchar[0] == 'O') itype = 6;
374            
375            if (junk1 != 0)
376            {
377                if (junk1 == 3 && itype == 8)
378                  itype = 41;  // N+
379                if (junk1 == 5 && itype == 6)
380                  itype = 42;  // O-
381            }
382        }
383
384        newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
385        if (newatom == -1)
386        {
387            printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
388            return FALSE;
389        }
390        if (isubred == 1)
391        {
392            atom[newatom].substr[0] = 0;
393            atom[newatom].substr[0] |= (1L << substr.icurstr);
394        }
395     }
396     has_Aromatic = FALSE;
397    
398     planar = TRUE;
399     if (xPlus && yPlus && zPlus) planar = FALSE;
400    
401     for (i=0; i < nibond; i++)
402     {
403        FetchRecord(infile,inputline);
404        for (j=0; j <4; j++)
405        {
406           c1[j] = inputline[j];
407           c2[j] = inputline[j+3];
408         }  
409        c1[3] = '\0';c2[3] = '\0';
410        ia1 = atoi(c1);
411        ia2 = atoi(c2);
412        istereo = 0;
413        sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
414        if (ibond >= 4)
415        {
416            ibond = 1;
417            has_Aromatic = TRUE;
418        }
419        make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
420        if (istereo == 1 && planar)
421        {
422            atom[ia2+ibotptr].z += 0.7;
423            if (atom[ia2+ibotptr].atomnum == 1)
424              atom[ia2+ibotptr].type = 60;
425        }
426        if (istereo == 6 && planar)
427        {
428            atom[ia2+ibotptr].z -= 0.7;
429            if (atom[ia2+ibotptr].atomnum == 1)
430             atom[ia2+ibotptr].type = 60;
431        }
432     }
433     FetchRecord(infile,inputline);  // M END line
434     FetchRecord(infile,inputline);
435     if (got_title == FALSE)
436        strcpy(Struct_Title,inputline);
437     fclose(infile);
438     ncount = strlen(Struct_Title);
439     for (i=0; i < ncount; i++)
440     {
441         if (Struct_Title[i] == '\n')
442         {
443             Struct_Title[i] = '\0';
444             break;
445         }
446     }
447 //   search for quaternary centers - check if flat
448     dz = 0.0;
449     for (i=1; i <= natom; i++)
450     {
451         dz += atom[i].z;
452         jji = jj1 = jj2 = jj3 = jj4 = 0;
453         if (atom[i].type != 0)
454         {
455             for (j=0; j < 6; j++)
456             {
457                 if (atom[i].iat[j] != 0)
458                   jji++;
459             }
460         }
461         if (jji == 4)
462         {
463             if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
464                   atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0)  // flat center
465             {
466                 for (j=0; j < 4; j++)
467                 {
468                     itemp[j] = 0;
469                     if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
470                     if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
471                     if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
472                     if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
473                 }
474                 icount = 0;
475                 if (jj1 == 1)
476                 {
477                     itemp[icount] = atom[i].iat[0];
478                     icount++;
479                 }
480                 if (jj2 == 1)
481                 {
482                     itemp[icount] = atom[i].iat[1];
483                     icount++;
484                 }
485                 if (jj3 == 1)
486                 {
487                     itemp[icount] = atom[i].iat[2];
488                     icount++;
489                 }
490                 if (jj3 == 1)
491                 {
492                     itemp[icount] = atom[i].iat[3];
493                     icount++;
494                 }
495                 if (icount >= 2)
496                 {
497                     atom[itemp[0]].z += 0.5;
498                     atom[itemp[1]].z -= 0.5;
499                 }
500             }
501         }
502     }
503     for (i=1; i < natom; i++)
504     {
505         for (j=i+1; j <= natom; j++)
506         {
507             xtmp = atom[i].x - atom[j].x;
508             ytmp = atom[i].y - atom[j].y;
509             ztmp = atom[i].z - atom[j].z;
510             if ( (xtmp + ytmp + ztmp) < 0.1)
511             {
512                 atom[i].x += 0.1;
513                 atom[i].y += 0.1;
514                 atom[i].z += 0.1;
515                 atom[j].x -= 0.1;
516                 atom[j].y -= 0.1;
517                 atom[j].z -= 0.1;
518             }
519         }
520     }
521              
522     if (has_Aromatic == TRUE)
523       mopaco(1,natom);
524     return TRUE;
525 }
526 // ==================================================
527 void read_schakal(int isel,int isubred)
528 {
529    int i, itype,  newatom, iz;
530    float xtmp, ytmp, ztmp;
531    char  dummy[30];
532    char  atomchar[6];
533    char  inputline[150];
534    FILE   *infile;
535    
536    infile = fopen_path(Openbox.path,Openbox.fname,"r");
537
538    if (infile == NULL)
539    {
540        message_alert("Error opening Schakal file","Schakal Setup");
541        return;
542    }
543    while ( FetchRecord(infile,inputline))
544    {
545        sscanf(inputline,"%s",dummy);
546        if (strcmp(dummy,"TITLE") == 0)
547        {
548            sscanf(inputline,"%s %s",dummy,Struct_Title);
549        } else if (strcmp(dummy,"ATOM") == 0)
550        {
551            sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
552            // strip off numbers
553            iz = strlen(atomchar);
554            for (i=0; i < iz; i++)
555            {
556                if (isdigit(atomchar[i]))
557                {
558                    atomchar[i] = '\0';
559                    break;
560                }
561            }
562            itype = 0;
563            newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
564        } else if (strcmp(dummy,"END") == 0)
565        {
566            break;
567        }
568    }
569    fclose(infile);
570    mopaco(1,natom);
571 //
462   }
573 // ==================================================
574 void write_schakal()
575 {
576    int i,j,lptest;
577    char atomchar[6];
578    FILE *wfile;
579    
580    lptest = 0;
581    for( i = 1; i <= natom; i++ )
582    {
583        if( atom[i].mmx_type == 20 )
584        {
585            lptest = 1;
586            hdel( lptest );
587            break;
588        }
589    }
590
591    /*     now write the schakal file */
592    if (files.append)
593        wfile = fopen_path(Savebox.path,Savebox.fname,"a");
594    else
595        wfile = fopen_path(Savebox.path,Savebox.fname,"w");
596
597    j = strlen(Struct_Title);
598    for (i=0; i < j; i++)
599    {
600        if (Struct_Title[i] == '\n')
601        {
602            Struct_Title[i] = '\0';
603            break;
604        }
605    }
606    fprintf(wfile,"TITLE  %s\n",Struct_Title);
607    fprintf(wfile,"CELL \n");
608    for (i=1; i <= natom; i++)
609    {
610        sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
611            
612        fprintf(wfile,"ATOM  %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
613           atom[i].y, atom[i].z);
614    }    
615    fprintf(wfile,"END \n");
616    fclose(wfile);
617 }
618
619

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines