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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines