Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

ReadPARM7.h

Go to the documentation of this file.
00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: ReadPARM7.h,v $
00005  *      $Author: johns $        $Locker:  $                $State: Exp $
00006  *      $Revision: 1.33 $      $Date: 2017/08/30 17:42:33 $
00007  *
00008  ***************************************************************************
00009  * DESCRIPTION:
00010  * NOTE:: Significant modifications were made to the VMD version of 
00011  *        Bill Ross's original code in order to make it easy to hook 
00012  *        into VMD plugin structures.  
00013  *        Further modifications were made to the VMD code to 
00014  *        read amber 7 parm files, courtesy of Brian Bennion
00015  * Here is what has changed:
00016  *     Functions became Class Methods, data became instance variables
00017  *     The Code to check for compressed files before opening was disabled
00018  *     Methods get_parm7_atom, get_parm7_bond, get_hydrogen_bond,
00019  *     get_parm7_natoms, get_parm7_nbonds, get_parm7_boxInfo were added in 
00020  *     order to convert from prm.c parlance to VMD conventions.
00021  ***************************************************************************/
00022 
00023 /*
00024  * COPYRIGHT 1992, REGENTS OF THE UNIVERSITY OF CALIFORNIA
00025  *
00026  *  prm.c - read information from an amber PARM topology file:
00027  *      atom/residue/bond/charge info, plus force field data.
00028  *      This file and the accompanying prm.h may be distributed
00029  *      provided this notice is retained unmodified and provided
00030  *      that any modifications to the rest of the file are noted
00031  *      in comments.
00032  *
00033  *      Bill Ross, UCSF 1994
00034  */
00035 
00036 #ifndef READPARM7_H
00037 #define READPARM7_H
00038 
00039 #include <stdio.h>
00040 #include <stdlib.h>
00041 #include <ctype.h>
00042 #include <sys/types.h>
00043 #include <sys/stat.h>
00044 #include <errno.h>
00045 #include <string.h>
00046 #include "molfile_plugin.h"  // needed for molfile return codes etc
00047 
00048 #if defined(WIN32) || defined(WIN64)
00049 #define strcasecmp  stricmp
00050 #define strncasecmp strnicmp
00051 #endif
00052 
00053 #if 0 
00054 #define _REAL   double
00055 #define DBLFMT  "%lf"
00056 #else
00057 #define _REAL   float
00058 #define DBLFMT  "%f"
00059 #endif
00060 
00061 
00062 typedef struct parm {
00063   char title[85];
00064   char version[85];
00065   int   IfBox, Nmxrs, IfCap,
00066         Natom,  Ntypes,  Nbonds, Nbonh,  Mbona,  Ntheth,  Mtheta, 
00067         Nphih,  Mphia,  Nhparm, Nparm, Nnb, Nres,Mptra,
00068         Nbona,  Ntheta,  Nphia,  Numbnd,  Numang,  Nptra,Jparm,
00069         Natyp,  Nphb, Nat3, Ntype2d, Nttyp, Nspm, Iptres, Nspsol,
00070         Ipatm, Natcap,Ifpert,Nbper,Ngper,Ndper,Mbper,Mgper,Mdper,
00071         Numextra;
00072   _REAL Box[3], Cutcap, Xcap, Ycap, Zcap;
00073 } parmstruct;
00074 
00075 static int read_parm7_flag(FILE *file, const char *flag, const char *format) {
00076   char buf[1024];
00077     
00078   /* read the %FLAG text */
00079   fscanf(file, "%s\n", buf);
00080   if (strcmp("%FLAG", buf)) {
00081     printf("AMBER 7 parm read error, at flag section %s,\n", flag);
00082     printf("        expected %%FLAG but got %s\n", buf);
00083     return 0; /* read of flag data failed */
00084   }
00085 
00086   /* read field name specifier */
00087   fscanf(file, "%s\n", buf);
00088   if (flag != NULL) {
00089     if (strcmp(flag, buf)) {
00090       printf("AMBER 7 parm read error at flag section %s,\n", flag);
00091       printf("      expected flag field %s but got %s\n", flag, buf);
00092       return 0; /* read of flag data failed */
00093     }
00094   }
00095 
00096   /* read format string */
00097   fscanf(file, "%s\n", buf);
00098   if (format != NULL) {
00099     if (strcmp(format, buf)) {
00100       if ( (!strcmp(flag, "TITLE") || !strcmp(flag, "CTITLE") )
00101           && !strcmp(format, "%FORMAT(20a4)")) {
00102         if (!strcmp(buf, "%FORMAT(a80)"))
00103           return 1; /* accept a80 as substitute for 20a4 */
00104       }
00105       printf("AMBER 7 parm read error at flag section %s,\n", flag);
00106       printf("      expected format %s but got %s\n", format, buf);
00107       return 0; /* read of flag data failed */
00108     }
00109   }
00110 
00111   return 1; /* read of flag data succeeded */
00112 }
00113 
00114 /*
00115  *  open_parm7_file() - fopen regular or popen compressed file for reading
00116  *  Return FILE handle on success.
00117  *  set as_pipe to 1 if opened with popen, or 0 if opened with fopen.
00118  */
00119 
00120 static FILE *open_parm7_file(const char *name, int *as_pipe) {
00121   struct stat buf;
00122   char cbuf[8192];
00123   int length;
00124   int &compressed = *as_pipe;
00125   FILE *fp;
00126 
00127   length = strlen(name);
00128   compressed = 0;  // Just to start
00129   strcpy(cbuf, name);
00130 
00131   /*
00132    *  if file doesn't exist, maybe it has been compressed/decompressed
00133    */
00134   if (stat(cbuf, &buf) == -1) {
00135     switch (errno) {
00136       case ENOENT:
00137         if (!compressed) {
00138           strcat(cbuf, ".Z");
00139           if (stat(cbuf, &buf) == -1) {
00140             printf("%s, %s: does not exist\n", name, cbuf);
00141             return(NULL);
00142           }
00143           compressed++;
00144 
00145           // Don't modify the filename
00146           //strcat(name, ".Z"); /* TODO: add protection */
00147         } else {
00148           cbuf[length-2] = '\0';
00149           if (stat(cbuf, &buf) == -1) {
00150             printf("%s, %s: does not exist\n", name, cbuf);
00151             return(NULL);
00152           }
00153           compressed = 0;
00154         }
00155         break;
00156 
00157       default:
00158         return(NULL);
00159     }
00160   }
00161 
00162   /*
00163    *  open the file
00164    */
00165 #if defined(_MSC_VER) || defined(__MINGW32__)
00166   if (compressed) {
00167     /* NO "zcat" on Win32 */
00168     printf("Cannot load compressed PARM files on Windows.\n");
00169     return NULL;
00170   }
00171 #else
00172   if (compressed) {
00173     char pcmd[sizeof(cbuf) + 7];
00174     sprintf(pcmd, "zcat '%s'", cbuf);
00175     if ((fp = popen(pcmd, "r")) == NULL) {
00176       perror(pcmd);
00177       return NULL;
00178     }
00179   }
00180 #endif
00181   else {
00182     if ((fp = fopen(cbuf, "r")) == NULL) {
00183       perror(cbuf);
00184       return NULL;
00185     }
00186   }
00187 
00188   return(fp);
00189 }
00190 
00191 
00192 static int parse_parm7_atoms(const char *fmt, 
00193     int natoms, molfile_atom_t *atoms, FILE *file) {
00194   char buf[85];
00195 
00196   if (strncasecmp(fmt, "%FORMAT(20a4)", 13))
00197     return 0;
00198 
00199   int j=0;
00200   for (int i=0; i<natoms; i++) {
00201     molfile_atom_t *atom = atoms+i;
00202     if (!(i%20)) {
00203       j=0;
00204       fgets(buf, 85, file);
00205     }
00206     strncpy(atom->name, buf+4*j, 4);
00207     atom->name[4]='\0';
00208     j++;
00209   }
00210   return 1;
00211 }
00212 
00213 
00214 static int parse_parm7_charge(const char *fmt, 
00215     int natoms, molfile_atom_t *atoms, FILE *file) {
00216   if (strncasecmp(fmt, "%FORMAT(5E16.8)", 15) &&
00217       strncasecmp(fmt, "%FORMAT(3E24.16)", 16)) 
00218     return 0;
00219 
00220   for (int i=0; i<natoms; i++) {
00221     double q=0;
00222     if (fscanf(file, " %lf", &q) != 1) {
00223       fprintf(stderr, "PARM7: error reading charge at index %d\n", i);
00224       return 0;
00225     }
00226     atoms[i].charge = 0.0548778 * (float)q; /* convert to elementary charge units */
00227   }
00228 
00229   return 1;
00230 }
00231 
00232 
00233 static int parse_parm7_atomicnumber(const char *fmt,
00234     int natoms, molfile_atom_t *atoms, FILE *file) {
00235   if (strncasecmp(fmt, "%FORMAT(10I8)", 13))
00236     return 0;
00237 
00238   for (int i=0; i<natoms; i++) {
00239     long int j=0;
00240     if (fscanf(file, "%li", &j) != 1) {
00241       fprintf(stderr, "PARM7: error reading atomic number at index %d\n", i);
00242       return 0;
00243     }
00244     atoms[i].atomicnumber = (int)j;
00245   }
00246 
00247   return 1;
00248 }
00249 
00250 
00251 static int parse_parm7_mass(const char *fmt,
00252     int natoms, molfile_atom_t *atoms, FILE *file) {
00253   if (strncasecmp(fmt, "%FORMAT(5E16.8)", 15)) return 0;
00254   for (int i=0; i<natoms; i++) {
00255     double m=0;
00256     if (fscanf(file, " %lf", &m) != 1) {
00257       fprintf(stderr, "PARM7: error reading mass at index %d\n", i);
00258       return 0;
00259     }
00260     atoms[i].mass = (float)m;
00261   }
00262   return 1;
00263 }
00264 
00265 
00266 static int parse_parm7_atype(const char *fmt,
00267     int natoms, molfile_atom_t *atoms, FILE *file) {
00268   if (strncasecmp(fmt, "%FORMAT(20a4)", 13)) return 0;
00269   char buf[85];
00270   int j=0;
00271   for (int i=0; i<natoms; i++) {
00272     molfile_atom_t *atom = atoms+i;
00273     if (!(i%20)) {
00274       j=0;
00275       fgets(buf, 85, file);
00276     }
00277     strncpy(atom->type, buf+4*j, 4);
00278     atom->type[4]='\0';
00279     j++;
00280   }
00281   return 1;
00282 }
00283 
00284 
00285 static int parse_parm7_resnames(const char *fmt,
00286     int nres, char *resnames, FILE *file) {
00287   if (strncasecmp(fmt, "%FORMAT(20a4)", 13)) return 0;
00288   char buf[85];
00289   int j=0;
00290   for (int i=0; i<nres; i++) {
00291     if (!(i%20)) {
00292       j=0;
00293       fgets(buf, 85, file);
00294     }
00295     strncpy(resnames, buf+4*j, 4);
00296     resnames += 4;
00297     j++;
00298   }
00299 
00300   return 1;
00301 }
00302 
00303 static int parse_parm7_respointers(const char *fmt, int natoms, 
00304     molfile_atom_t *atoms, int nres, const char *resnames, FILE *file) {
00305   if (strncasecmp(fmt, "%FORMAT(10I8)", 13)) return 0;
00306   int cur, next;
00307   fscanf(file, " %d", &cur);
00308   for (int i=1; i<nres; i++) {
00309     if (fscanf(file, " %d", &next) != 1) {
00310       fprintf(stderr, "PARM7: error reading respointer records at residue %d\n", i);
00311       return 0;
00312     }
00313     while (cur < next) {
00314       if (cur > natoms) {
00315         fprintf(stderr, "invalid atom index: %d\n", cur);
00316         return 0;
00317       }
00318       strncpy(atoms[cur-1].resname, resnames, 4);
00319       atoms[cur-1].resname[4] = '\0';
00320       atoms[cur-1].resid = i;
00321       cur++;
00322     }
00323     resnames += 4;
00324   }
00325   // store the last residue name
00326   while (cur <= natoms) {
00327     strncpy(atoms[cur-1].resname, resnames, 4);
00328     atoms[cur-1].resname[4] = '\0';
00329     atoms[cur-1].resid = nres;
00330     cur++;
00331   }
00332 
00333   return 1;
00334 }
00335 
00336 
00337 static int parse_parm7_bonds(const char *fmt,
00338     int nbonds, int *from, int *to, FILE *file) {
00339   if (strncasecmp(fmt, "%FORMAT(10I8)", 13)) 
00340     return 0;
00341 
00342   int a, b, tmp;
00343   for (int i=0; i<nbonds; i++) {
00344     if (fscanf(file, " %d %d %d", &a, &b, &tmp) != 3) {
00345       fprintf(stderr, "PARM7: error reading bond number %d\n", i);
00346       return 0;
00347     }
00348     from[i] = a/3 + 1;
00349     to[i]   = b/3 + 1;
00350   }
00351 
00352   return 1;
00353 }
00354 
00355 
00356 /*
00357  *  close_parm7_file() - close fopened or popened file
00358  */
00359 static void close_parm7_file(FILE *fileptr, int popn) {
00360 #if defined(_MSC_VER) || defined(__MINGW32__)
00361   if (popn) {
00362     printf("pclose() no such function on win32!\n");
00363   } else {
00364    if (fclose(fileptr) == -1)
00365      perror("fclose");
00366   }
00367 #else
00368   if (popn) {
00369     if (pclose(fileptr) == -1)
00370       perror("pclose");
00371   } else {
00372     if (fclose(fileptr) == -1)
00373       perror("fclose");
00374   }
00375 #endif
00376 }
00377 
00378 static const char *parm7 = "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n";
00379 
00380 static parmstruct *read_parm7_header(FILE *file) {
00381   char sdum[512]; 
00382   parmstruct *prm;
00383   prm = new parmstruct;
00384 
00385   /* READ VERSION */
00386   fgets(sdum, 512, file);
00387 
00388   /* READ TITLE */
00389   fscanf(file, "%s\n", sdum); // "%FLAG"
00390   if (strcmp("%FLAG", sdum)) {
00391     printf("AMBER 7 parm read error, can't find TITLE flag.\n");
00392     printf("        expected %%FLAG, got %s\n", sdum);
00393     delete prm;
00394     return NULL;
00395   }
00396   fscanf(file, "%s\n", sdum); // "TITLE" or "CTITLE"
00397   if (strcmp("TITLE", sdum) && strcmp("CTITLE", sdum)) {
00398     printf("AMBER 7 parm read error, at flag section TITLE,\n");
00399     printf("        expected TITLE or CTITLE but got %s,\n", sdum);
00400     delete prm;
00401     return NULL;
00402   }
00403   fscanf(file, "%s\n", sdum); // "FORMAT (20a4)"
00404   if (strcmp(sdum, "%FORMAT(20a4)") && strcmp(sdum, "%FORMAT(a80)")) {
00405     printf("AMBER 7 parm read error, at flag section TITLE,\n");
00406     printf("        expected %%FLAG but got %s,\n", sdum);
00407     delete prm;
00408     return NULL;
00409   }
00410 
00411   // read the title string itself, and handle empty lines
00412 #if 0
00413   // XXX this code fails with some AMBER 9 test files
00414   fscanf(file, "%s\n", prm->title);
00415 #else
00416   // XXX this hack causes AMBER 9 prmtop files to load
00417   fgets(prm->title, sizeof(prm->title), file);
00418 #endif
00419 
00420   if (strstr(prm->title, "%FLAG") == NULL) {
00421     // Got a title string
00422     if (!read_parm7_flag(file, "POINTERS", "%FORMAT(10I8)")) {
00423       delete prm;
00424       return NULL;
00425     }
00426   } else {
00427     // NO title string, use a special method to pick up next flag
00428 #if 0
00429     fscanf(file,"%s\n", sdum);
00430     if (strcmp("POINTERS", sdum)) {
00431       printf("AMBER 7 parm read error at flag section POINTERS\n");
00432       printf("      expected flag field POINTERS but got %s\n", sdum);
00433 #else
00434     if (strstr(prm->title, "POINTERS") == NULL) {
00435       printf("AMBER 7 parm read error at flag section POINTERS\n");
00436       printf("      expected flag field POINTERS but got %s\n", prm->title);
00437 #endif
00438       delete prm;
00439       return NULL;
00440     }
00441 #if 0
00442     fscanf(file,"%s\n", sdum);
00443     if (strcasecmp("%FORMAT(10I8)", sdum)) {
00444 #else
00445     fgets(sdum, sizeof(sdum), file);
00446     if ((strstr(sdum, "%FORMAT(10I8)") == NULL) &&
00447         (strstr(sdum, "%FORMAT(10i8)") == NULL)) {
00448 #endif
00449       printf("AMBER 7 parm read error at flag section POINTERS,\n");
00450       printf("      expected format %%FORMAT(10I8) but got %s\n", sdum);
00451       delete prm;
00452       return NULL;
00453     }
00454   }
00455 
00456   /* READ POINTERS (CONTROL INTEGERS) */
00457   fscanf(file,parm7,
00458          &prm->Natom,  &prm->Ntypes, &prm->Nbonh, &prm->Nbona,
00459          &prm->Ntheth, &prm->Ntheta, &prm->Nphih, &prm->Nphia,
00460          &prm->Jparm,  &prm->Nparm);
00461   fscanf(file, parm7,  
00462          &prm->Nnb,   &prm->Nres,   &prm->Mbona,  &prm->Mtheta,
00463          &prm->Mphia, &prm->Numbnd, &prm->Numang, &prm->Mptra,
00464          &prm->Natyp, &prm->Nphb);
00465   fscanf(file, parm7,  &prm->Ifpert, &prm->Nbper,  &prm->Ngper,
00466          &prm->Ndper, &prm->Mbper,  &prm->Mgper, &prm->Mdper,
00467          &prm->IfBox, &prm->Nmxrs,  &prm->IfCap);
00468 
00469   fscanf(file,"%8d",&prm->Numextra); //BB
00470   prm->Nptra=prm->Mptra; //BB new to amber 7 files...
00471 
00472   prm->Nat3 = 3 * prm->Natom;
00473   prm->Ntype2d = prm->Ntypes * prm->Ntypes;
00474   prm->Nttyp = prm->Ntypes*(prm->Ntypes+1)/2;
00475 
00476   return prm;
00477 }
00478 
00479 
00480 #endif

Generated on Tue Jan 21 02:56:22 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002