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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines