ReadPARM.h Source File

00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: ReadPARM.h,v $
00005  *      $Author: johns $        $Locker:  $                $State: Exp $
00006  *      $Revision: 1.15 $      $Date: 2016/11/06 17:54:57 $
00007  *
00008  ***************************************************************************
00009  * DESCRIPTION:
00010  * NOTE:: Significant 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  * Here is what was changed:
00014  *     Functions became Class Methods, data became instance variables
00015  *     The Code to check for compressed files before opening was disabled
00016  *     Methods get_parm_atom, get_parm_bond, get_hydrogen_bond,
00017  *     get_parm_natoms, get_parm_nbonds, get_parm_boxInfo were added in 
00018  *     order to convert from prm.c parlance to VMD conventions.
00019  *     RCS Information headers and footers were added.
00020  ***************************************************************************/
00021 
00022 /*
00023  * COPYRIGHT 1992, REGENTS OF THE UNIVERSITY OF CALIFORNIA
00024  *
00025  *  prm.c - read information from an amber PARM topology file:
00026  *      atom/residue/bond/charge info, plus force field data.
00027  *      This file and the accompanying prm.h may be distributed
00028  *      provided this notice is retained unmodified and provided
00029  *      that any modifications to the rest of the file are noted
00030  *      in comments.
00031  *
00032  *      Bill Ross, UCSF 1994
00033  */
00034 
00035 #ifndef READPARM_H
00036 #define READPARM_H
00037 
00038 // XXX enable the new AMBER reading code, deals with packed 
00039 // fortran 12I6 integer formats.
00040 #define USENEWCODE 1
00041 
00042 #include <stdio.h>
00043 #include <stdlib.h>
00044 #include <ctype.h>
00045 #include <sys/types.h>
00046 #include <sys/stat.h>
00047 #include <errno.h>
00048 #include <string.h>
00049 #include "molfile_plugin.h"  // needed for molfile return codes etc
00050 
00051 #if 0 
00052 #define _REAL           double
00053 #define DBLFMT          "%lf"
00054 #else
00055 #define _REAL           float
00056 #define DBLFMT          "%f"
00057 #endif
00058 
00059 typedef struct parm {
00060         char    ititl[81];
00061         int     IfBox, Nmxrs, IfCap,
00062                  Natom,  Ntypes,  Nbonh,  Mbona,  Ntheth,  Mtheta, 
00063                  Nphih,  Mphia,  Nhparm, Nparm, Nnb, Nres,
00064                  Nbona,  Ntheta,  Nphia,  Numbnd,  Numang,  Nptra,
00065                  Natyp,  Nphb, Nat3, Ntype2d, Nttyp, Nspm, Iptres, Nspsol,
00066                  Ipatm, Natcap;
00067         char    *AtomNames, *ResNames, *AtomSym, *AtomTree;
00068         _REAL   *Charges, *Masses, *Rk, *Req, *Tk, *Teq, *Pk, *Pn, *Phase,
00069                  *Solty, *Cn1, *Cn2, *HB12, *HB6;
00070         _REAL   Box[3], Cutcap, Xcap, Ycap, Zcap;
00071         int     *Iac, *Iblo, *Cno, *Ipres, *ExclAt, *TreeJoin, 
00072                 *AtomRes, *BondHAt1, *BondHAt2, *BondHNum, *BondAt1, *BondAt2, 
00073                 *BondNum, *AngleHAt1, *AngleHAt2, *AngleHAt3, *AngleHNum, 
00074                 *AngleAt1, *AngleAt2, *AngleAt3, *AngleNum, *DihHAt1, 
00075                 *DihHAt2, *DihHAt3, *DihHAt4, *DihHNum, *DihAt1, *DihAt2, 
00076                 *DihAt3, *DihAt4, *DihNum, *Boundary;
00077 } parmstruct;
00078 
00079 // put the class in an anonymous namespace to give it internal linkage
00080 namespace {
00081   class ReadPARM {
00082   public:
00083     ReadPARM() {popn = 0;}
00084     ~ReadPARM(void) {}
00085  
00086     int popn;
00087     parmstruct *prm;
00088     FILE *open_parm_file(const char *name);
00089     void close_parm_file(FILE *fileptr);
00090     char *get(int size);
00091     int preadln(FILE *file, char *string);
00092     int readparm(FILE *file);
00093     void get_parm_atom(int, char *, char *, char *, char *, int *, float *,
00094                        float *);
00095 
00096     void get_parm_bond(int, int fromAtom[], int toAtom[]);
00097     void get_hydrogen_bond(int, int fromAtom[], int toAtom[]);
00098     int get_parm_natoms();
00099     int get_parm_nbonds();
00100     int get_parm_boxInfo();
00101     int read_fortran_12I6(FILE *fp, int *data, int count);  
00102   };
00103 }
00104 
00105 
00106 /*      fortran formats 
00107  *       9118 FORMAT(12I6)
00108  *       9128 FORMAT(5E16.8)
00109 static char     *f9118 = (char *) "%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d\n";
00110  */
00111 
00112 static int readtoeoln(FILE *f) {
00113   int c;
00114 
00115   /* skip to eoln */
00116   while((c = getc(f)) != '\n') {
00117     if (c == EOF) 
00118       return -1;
00119   }
00120 
00121   return 0;
00122 }  
00123 
00124 /***********************************************************************
00125                                 open_parm_file()
00126 ************************************************************************/
00127 
00128 /*
00129  *  open_parm_file() - fopen regular or popen compressed file for reading
00130  */
00131 
00132 FILE *ReadPARM::open_parm_file(const char *name) {
00133         struct stat     buf;
00134         char            cbuf[120];
00135         int             length, compressed;
00136         FILE            *fp;
00137 
00138         length = strlen(name);
00139         compressed = 0;  // Just to start
00140 //      compressed = iscompressed(name);
00141         strcpy(cbuf, name);
00142 
00143         /*
00144          *  if file doesn't exist, maybe it has been compressed/decompressed
00145          */
00146 
00147         if (stat(cbuf, &buf) == -1) {
00148                 switch (errno) {
00149                 case ENOENT:    {
00150                         if (!compressed) {
00151                                 strcat(cbuf, ".Z");
00152                                 if (stat(cbuf, &buf) == -1) {
00153                                         printf("%s, %s: does not exist\n", 
00154                                                 name, cbuf);
00155                                         return(NULL);
00156                                 }
00157                                 compressed++;
00158                                 // Don't modify the filename
00159                                 //strcat(name, ".Z"); /* TODO: add protection */
00160                         } else {
00161                                 cbuf[length-2] = '\0';
00162                                 if (stat(cbuf, &buf) == -1) {
00163                                         printf("%s, %s: does not exist\n", 
00164                                                         name, cbuf);
00165                                         return(NULL);
00166                                 }
00167                                 compressed = 0;
00168                         }
00169                         break;
00170                 }
00171                 default:
00172                         return(NULL);
00173                 }
00174         }
00175 
00176         /*
00177          *  open the file
00178          */
00179 #if defined(_MSC_VER) || defined(__MINGW32__)
00180         if (compressed) {
00181           /* NO "zcat" on Win32 */
00182           printf("Cannot load compressed PARM files on Windows.\n");
00183           return NULL;
00184         }
00185 #else
00186         if (compressed) {
00187                 char pcmd[sizeof(cbuf) + 7];
00188                 popn = 1;
00189 
00190                 sprintf(pcmd, "zcat '%s'", cbuf);
00191                 if ((fp = popen(pcmd, "r")) == NULL) {
00192                         perror(pcmd);
00193                         return NULL;
00194                 }
00195         }
00196 #endif
00197         else {
00198                 if ((fp = fopen(cbuf, "r")) == NULL) {
00199                         perror(cbuf);
00200                         return NULL;
00201                 }
00202         }
00203         return(fp);
00204 }
00205 
00206 /***********************************************************************
00207                             close_parm_file   
00208 ************************************************************************/
00209 
00210 /*
00211  *  close_parm_file() - close fopened or popened file
00212  */
00213 void ReadPARM::close_parm_file(FILE *fileptr) {
00214 #if defined(_MSC_VER) || defined(__MINGW32__)
00215   if (popn) {
00216     printf("pclose() no such function on win32!\n");
00217   } else {
00218     if (fclose(fileptr) == -1)
00219       perror("fclose");
00220   }
00221 #else
00222   if (popn) {
00223     if (pclose(fileptr) == -1)
00224       perror("pclose");
00225   } else {
00226     if (fclose(fileptr) == -1)
00227       perror("fclose");
00228   }
00229 #endif
00230 }
00231 
00232 
00233 /***********************************************************************
00234                                       GET()
00235 ************************************************************************/
00236 char *ReadPARM::get(int size) {
00237         char    *ptr;
00238 
00239 #ifdef DEBUG
00240         printf("malloc %d\n", size);
00241         fflush(stdout);
00242 #endif
00243         if (size ==0)
00244                 return((char *) NULL);
00245 
00246         if ((ptr = (char *) malloc((unsigned)size)) == NULL) {
00247                 printf("malloc %d", size);
00248                 fflush(stdout);
00249                 perror("malloc err:");
00250         fprintf(stderr, "Exiting due to ReadPARM memory allocation error.\n");
00251         }
00252         return(ptr);
00253 }
00254 
00255 /***********************************************************************
00256                                 PREADLN()
00257 ************************************************************************/
00258 int ReadPARM::preadln(FILE *file, char *string) {
00259   int i, j;
00260 
00261   for (i=0; i<81; i++) {
00262     if ((j = getc(file)) == EOF) {
00263       printf("Error: unexpected EOF in Parm file\n");
00264       return MOLFILE_ERROR;
00265     }
00266 
00267     string[i] = (char) j;
00268     if (string[i] == '\n') {
00269       break;
00270     }
00271   }
00272 
00273   if (i == 80  &&  string[i] != '\n') {
00274     printf("Error: line too long in Parm file:\n%.80s", string);
00275     return MOLFILE_ERROR;
00276   }
00277 
00278   return 0; /* success */
00279 }
00280 
00281 /***************************************************************************
00282                                 READPARM()
00283 ****************************************************************************/
00284 /*
00285  * readparm() - instantiate a given parmstruct
00286  */
00287 int ReadPARM::readparm(FILE *file) {
00288         //
00289         // XXX This code leaks memory every time it's run!  prm is allocated
00290         // but its data is never freed, even when ReadPARM itself is deleted.
00291         // I've added exception code so we at least have a chance of 
00292         // recovering gracefully from a memory allocation error.
00293 
00294         _REAL           *H;
00295         int             i, res, ifpert;
00296         int *buffer; // used for reading fortran integer blocks
00297 
00298         prm = (parmstruct *) get(sizeof(parmstruct));
00299         if (prm == NULL) {
00300           return MOLFILE_ERROR; /* allocation failure */
00301         }
00302 
00303         /* READ TITLE */
00304         if (preadln(file, prm->ititl) != 0) {
00305           return MOLFILE_ERROR;
00306         }
00307 
00308         /* READ CONTROL INTEGERS */
00309 #if !defined(USENEWCODE)
00310         fscanf(file, f9118, 
00311                 &prm->Natom,  &prm->Ntypes, &prm->Nbonh, &prm->Mbona, 
00312                 &prm->Ntheth, &prm->Mtheta, &prm->Nphih, &prm->Mphia, 
00313                 &prm->Nhparm, &prm->Nparm,  &prm->Nnb,   &prm->Nres);
00314 
00315         fscanf(file, f9118, 
00316                 &prm->Nbona,  &prm->Ntheta, &prm->Nphia, &prm->Numbnd, 
00317                 &prm->Numang, &prm->Nptra,  &prm->Natyp, &prm->Nphb, 
00318                 &ifpert,      &idum,        &idum,       &idum);
00319 
00320         fscanf(file, " %d %d %d %d %d %d", 
00321                 &idum, &idum,&idum,&prm->IfBox,&prm->Nmxrs,&prm->IfCap);
00322 #else
00323         buffer = new int[30];
00324         if (!read_fortran_12I6(file,buffer,30)) { 
00325           return MOLFILE_ERROR;
00326         }
00327         prm->Natom = buffer[0];
00328         prm->Ntypes = buffer[1];
00329         prm->Nbonh = buffer[2];
00330         prm->Mbona = buffer[3];
00331         prm->Ntheth = buffer[4];
00332         prm->Mtheta = buffer[5];
00333         prm->Nphih = buffer[6];
00334         prm->Mphia = buffer[7];
00335         prm->Nhparm = buffer[8];
00336         prm->Nparm = buffer[9];
00337         prm->Nnb = buffer[10];
00338         prm->Nres = buffer[11];
00339         prm->Nbona = buffer[12];
00340         prm->Ntheta = buffer[13];
00341         prm->Nphia = buffer[14];
00342         prm->Numbnd = buffer[15];
00343         prm->Numang = buffer[16];
00344         prm->Nptra = buffer[17];
00345         prm->Natyp = buffer[18];
00346         prm->Nphb = buffer[19];
00347         ifpert = buffer[20];
00348         // items 21 through 26 are ignored currently.
00349         prm->IfBox = buffer[27];
00350         prm->Nmxrs = buffer[28];
00351         prm->IfCap = buffer[29];
00352         delete [] buffer;
00353 #endif
00354         readtoeoln(file);
00355 
00356 
00357        if (ifpert) {
00358          printf("not equipped to read perturbation prmtop\n");
00359          free(prm);
00360          return MOLFILE_ERROR;
00361        }
00362 
00363 
00364         /* ALLOCATE MEMORY */
00365         prm->Nat3 = 3 * prm->Natom;
00366         prm->Ntype2d = prm->Ntypes * prm->Ntypes;
00367         prm->Nttyp = prm->Ntypes*(prm->Ntypes+1)/2;
00368 
00369         /*
00370          * get most of the indirect stuff; some extra allowed for char arrays
00371          */
00372         prm->AtomNames = (char *) get(4*prm->Natom+81);
00373         prm->Charges = (_REAL *) get(sizeof(_REAL)*prm->Natom);
00374         prm->Masses = (_REAL *) get(sizeof(_REAL)*prm->Natom);
00375         prm->Iac = (int *) get(sizeof(int)*prm->Natom);
00376         prm->Iblo = (int *) get(sizeof(int)*prm->Natom);
00377         prm->Cno = (int *) get(sizeof(int)* prm->Ntype2d);
00378         prm->ResNames = (char *) get(4* prm->Nres+81);
00379         prm->Ipres = (int *) get(sizeof(int)*( prm->Nres+1));
00380         prm->Rk = (_REAL *) get(sizeof(_REAL)* prm->Numbnd);
00381         prm->Req = (_REAL *) get(sizeof(_REAL)* prm->Numbnd);
00382         prm->Tk = (_REAL *) get(sizeof(_REAL)* prm->Numang);
00383         prm->Teq = (_REAL *) get(sizeof(_REAL)* prm->Numang);
00384         prm->Pk = (_REAL *) get(sizeof(_REAL)* prm->Nptra);
00385         prm->Pn = (_REAL *) get(sizeof(_REAL)* prm->Nptra);
00386         prm->Phase = (_REAL *) get(sizeof(_REAL)* prm->Nptra);
00387         prm->Solty = (_REAL *) get(sizeof(_REAL)* prm->Natyp);
00388         prm->Cn1 = (_REAL *) get(sizeof(_REAL)* prm->Nttyp);
00389         prm->Cn2 = (_REAL *) get(sizeof(_REAL)* prm->Nttyp);
00390         prm->BondHAt1 = (int *) get(sizeof(int)* prm->Nbonh);
00391         prm->BondHAt2 = (int *) get(sizeof(int)* prm->Nbonh);
00392         prm->BondHNum = (int *) get(sizeof(int)* prm->Nbonh);
00393         prm->BondAt1 = (int *) get(sizeof(int)* prm->Nbona);
00394         prm->BondAt2 = (int *) get(sizeof(int)* prm->Nbona);
00395         prm->BondNum = (int *) get(sizeof(int)* prm->Nbona);
00396         prm->AngleHAt1 = (int *) get(sizeof(int)* prm->Ntheth);
00397         prm->AngleHAt2 = (int *) get(sizeof(int)* prm->Ntheth);
00398         prm->AngleHAt3 = (int *) get(sizeof(int)* prm->Ntheth);
00399         prm->AngleHNum = (int *) get(sizeof(int)* prm->Ntheth);
00400         prm->AngleAt1 = (int *) get(sizeof(int)* prm->Ntheta);
00401         prm->AngleAt2 = (int *) get(sizeof(int)*prm->Ntheta);
00402         prm->AngleAt3 = (int *) get(sizeof(int)*prm->Ntheta);
00403         prm->AngleNum = (int *) get(sizeof(int)*prm->Ntheta);
00404         prm->DihHAt1 = (int *) get(sizeof(int)*prm->Nphih);
00405         prm->DihHAt2 = (int *) get(sizeof(int)*prm->Nphih);
00406         prm->DihHAt3 = (int *) get(sizeof(int)*prm->Nphih);
00407         prm->DihHAt4 = (int *) get(sizeof(int)*prm->Nphih);
00408         prm->DihHNum = (int *) get(sizeof(int)*prm->Nphih);
00409         prm->DihAt1 = (int *) get(sizeof(int)*prm->Nphia);
00410         prm->DihAt2 = (int *) get(sizeof(int)*prm->Nphia);
00411         prm->DihAt3 = (int *) get(sizeof(int)*prm->Nphia);
00412         prm->DihAt4 = (int *) get(sizeof(int)*prm->Nphia);
00413         prm->DihNum = (int *) get(sizeof(int)*prm->Nphia);
00414         prm->ExclAt = (int *) get(sizeof(int)*prm->Nnb);
00415         prm->HB12 = (_REAL *) get(sizeof(_REAL)*prm->Nphb);
00416         prm->HB6 = (_REAL *) get(sizeof(_REAL)*prm->Nphb);
00417         prm->AtomSym = (char *) get(4*prm->Natom+81);
00418         prm->AtomTree = (char *) get(4*prm->Natom+81);
00419         prm->TreeJoin = (int *) get(sizeof(int)*prm->Natom);
00420         prm->AtomRes = (int *) get(sizeof(int)*prm->Natom);
00421 
00422         /* 
00423          * READ ATOM NAMES -IH(M04)
00424          */
00425         for (i=0; i<(prm->Natom/20 + (prm->Natom%20 ? 1 : 0)); i++) {
00426                 preadln(file, &prm->AtomNames[i*80]);
00427         }
00428 
00429         /* 
00430          * READ ATOM CHARGES -X(L15)
00431          *      (pre-multiplied by an energy factor of 18.2223 == sqrt(332)
00432          *       for faster force field calculations)
00433          */
00434         for (i=0; i<prm->Natom; i++) {
00435           fscanf(file, " " DBLFMT, &prm->Charges[i]);
00436           prm->Charges[i] *= 0.0548778; /* convert back to elementary charge units */
00437         }
00438         readtoeoln(file);
00439 
00440         /* 
00441          * READ ATOM MASSES -X(L20)
00442          */
00443         for (i=0; i<prm->Natom; i++)
00444           fscanf(file, " " DBLFMT, &prm->Masses[i]);
00445         readtoeoln(file);
00446 
00447         /* 
00448          * READ ATOM L-J TYPES -IX(I04)
00449          */
00450 #if !defined(USENEWCODE)
00451         for (i=0; i<prm->Natom; i++) 
00452           fscanf(file, " %d", &prm->Iac[i]);
00453 #else
00454         if (!read_fortran_12I6(file, prm->Iac, prm->Natom)) {
00455           return MOLFILE_ERROR;
00456         }
00457 #endif
00458         readtoeoln(file);
00459 
00460         /* 
00461          * READ ATOM INDEX TO 1st IN EXCLUDED ATOM LIST "NATEX" -IX(I08)
00462          */
00463 #if !defined(USENEWCODE)
00464         for (i=0; i<prm->Natom; i++)
00465           fscanf(file, " %d", &prm->Iblo[i]);
00466 #else
00467         if (!read_fortran_12I6(file, prm->Iblo, prm->Natom)) {
00468           return MOLFILE_ERROR;
00469         }
00470 #endif
00471         readtoeoln(file);
00472 
00473         /* 
00474          * READ TYPE INDEX TO N-B TYPE -IX(I06)
00475          */
00476 #if !defined(USENEWCODE)
00477         for (i=0; i<prm->Ntype2d; i++)
00478           fscanf(file, " %d", &prm->Cno[i]);
00479 #else
00480         if (!read_fortran_12I6(file, prm->Cno, prm->Ntype2d)) {
00481           return MOLFILE_ERROR;
00482         } 
00483 #endif
00484         readtoeoln(file);
00485 
00486         /* 
00487          * READ RES NAMES (4 chars each, 4th blank) -IH(M02)
00488          */
00489         for (i=0; i<(prm->Nres/20 + (prm->Nres%20 ? 1 : 0)); i++)
00490           preadln(file, &prm->ResNames[i*80]);
00491 
00492         /* 
00493          * READ RES POINTERS TO 1st ATOM                -IX(I02)
00494          */
00495 #if !defined(USENEWCODE)
00496         for (i=0; i<prm->Nres; i++) 
00497           fscanf(file, " %d", &prm->Ipres[i]);
00498 #else
00499         if (!read_fortran_12I6(file, prm->Ipres, prm->Nres)) {
00500           return MOLFILE_ERROR;
00501         }
00502 #endif
00503         prm->Ipres[prm->Nres] = prm->Natom + 1;
00504         readtoeoln(file);
00505 
00506         /* 
00507          * READ BOND FORCE CONSTANTS                    -RK()
00508          */
00509         for (i=0; i< prm->Numbnd; i++) 
00510           fscanf(file, " " DBLFMT, &prm->Rk[i]);
00511         readtoeoln(file);
00512 
00513         /* 
00514          * READ BOND LENGTH OF MINIMUM ENERGY           -REQ()
00515          */
00516         for (i=0; i< prm->Numbnd; i++) 
00517           fscanf(file, " " DBLFMT, &prm->Req[i]);
00518         readtoeoln(file);
00519 
00520         /* 
00521          * READ BOND ANGLE FORCE CONSTANTS (following Rk nomen) -TK()
00522          */
00523         for (i=0; i< prm->Numang; i++) 
00524           fscanf(file, " " DBLFMT, &prm->Tk[i]);
00525         readtoeoln(file);
00526 
00527         /* 
00528          * READ BOND ANGLE OF MINIMUM ENERGY (following Req nomen) -TEQ()
00529          */
00530         for (i=0; i< prm->Numang; i++) 
00531           fscanf(file, " " DBLFMT, &prm->Teq[i]);
00532         readtoeoln(file);
00533 
00534         /* 
00535          * READ DIHEDRAL PEAK MAGNITUDE                 -PK()
00536          */
00537         for (i=0; i< prm->Nptra; i++) 
00538           fscanf(file, " " DBLFMT, &prm->Pk[i]);
00539         readtoeoln(file);
00540 
00541         /* 
00542          * READ DIHEDRAL PERIODICITY                    -PN()
00543          */
00544         for (i=0; i< prm->Nptra; i++) 
00545           fscanf(file, " " DBLFMT, &prm->Pn[i]);
00546         readtoeoln(file);
00547 
00548         /* 
00549          * READ DIHEDRAL PHASE                          -PHASE()
00550          */
00551         for (i=0; i< prm->Nptra; i++) 
00552           fscanf(file, " " DBLFMT, &prm->Phase[i]);
00553         readtoeoln(file);
00554 
00555         /* 
00556          * ?? "RESERVED"                                -SOLTY()
00557          */
00558         for (i=0; i< prm->Natyp; i++) 
00559           fscanf(file, " " DBLFMT, &prm->Solty[i]);
00560         readtoeoln(file);
00561 
00562         /* 
00563          * READ L-J R**12 FOR ALL PAIRS OF ATOM TYPES   -CN1()
00564          *      (SHOULD BE 0 WHERE H-BONDS)
00565          */
00566         for (i=0; i< prm->Nttyp; i++) 
00567           fscanf(file, " " DBLFMT, &prm->Cn1[i]);
00568         readtoeoln(file);
00569 
00570         /* 
00571          * READ L-J R**6 FOR ALL PAIRS OF ATOM TYPES    -CN2()
00572          *      (SHOULD BE 0 WHERE H-BONDS)
00573          */
00574         for (i=0; i< prm->Nttyp; i++) 
00575           fscanf(file, " " DBLFMT, &prm->Cn2[i]);
00576         readtoeoln(file);
00577 
00578         /* 
00579          * READ COVALENT BOND W/ HYDROGEN (3*(atnum-1)): 
00580          *      IBH = ATOM1             -IX(I12)
00581          *      JBH = ATOM2             -IX(I14)
00582          *      ICBH = BOND ARRAY PTR   -IX(I16)
00583          */
00584 #if !defined(USENEWCODE)
00585         for (i=0; i<prm->Nbonh; i++) 
00586           fscanf(file, " %d %d %d", 
00587                  &prm->BondHAt1[i], &prm->BondHAt2[i], &prm->BondHNum[i]);
00588 #else
00589         buffer = new int[3*prm->Nbonh];
00590         if (!read_fortran_12I6(file, buffer, 3*prm->Nbonh)) { 
00591           return MOLFILE_ERROR;
00592         }
00593         for (i=0; i<prm->Nbonh; i++) { 
00594           prm->BondHAt1[i] = buffer[3*i];
00595           prm->BondHAt2[i] = buffer[3*i+1];
00596           prm->BondHNum[i] = buffer[3*i+2];
00597         }
00598         delete [] buffer;
00599 #endif
00600         readtoeoln(file);
00601 
00602         /* 
00603          * READ COVALENT BOND W/OUT HYDROGEN (3*(atnum-1)):
00604          *      IB = ATOM1              -IX(I18)
00605          *      JB = ATOM2              -IX(I20)
00606          *      ICB = BOND ARRAY PTR    -IX(I22)
00607          */
00608 #if !defined(USENEWCODE)
00609         for (i=0; i<prm->Nbona; i++)
00610           fscanf(file, " %d %d %d", 
00611                  &prm->BondAt1[i], &prm->BondAt2[i], &prm->BondNum[i]);
00612 #else
00613         buffer = new int[3*prm->Nbona];
00614         if (!read_fortran_12I6(file, buffer, 3*prm->Nbona)) { 
00615           return MOLFILE_ERROR;
00616         }
00617         for (i=0; i<prm->Nbona; i++) { 
00618           prm->BondAt1[i] = buffer[3*i];
00619           prm->BondAt2[i] = buffer[3*i+1];
00620           prm->BondNum[i] = buffer[3*i+2];
00621         }
00622         delete [] buffer;
00623 #endif
00624         readtoeoln(file);
00625 
00626         /* 
00627          * READ ANGLE W/ HYDROGEN: 
00628          *      ITH = ATOM1                     -IX(I24)
00629          *      JTH = ATOM2                     -IX(I26)
00630          *      KTH = ATOM3                     -IX(I28)
00631          *      ICTH = ANGLE ARRAY PTR          -IX(I30)
00632          */
00633 #if !defined(USENEWCODE)
00634         for (i=0; i<prm->Ntheth; i++)
00635           fscanf(file, " %d %d %d %d", 
00636                  &prm->AngleHAt1[i], &prm->AngleHAt2[i], 
00637                  &prm->AngleHAt3[i], &prm->AngleHNum[i]);
00638 #else
00639         buffer = new int[4*prm->Ntheth];
00640         if (!read_fortran_12I6(file, buffer, 4*prm->Ntheth)) { 
00641           return MOLFILE_ERROR;
00642         }
00643         for (i=0; i<prm->Ntheth; i++) { 
00644           prm->AngleHAt1[i] = buffer[4*i];
00645           prm->AngleHAt2[i] = buffer[4*i+1];
00646           prm->AngleHAt3[i] = buffer[4*i+2];
00647           prm->AngleHNum[i] = buffer[4*i+3];
00648         }
00649         delete [] buffer;
00650 #endif
00651         readtoeoln(file);
00652 
00653         /* 
00654          * READ ANGLE W/OUT HYDROGEN: 
00655          *      IT = ATOM1                      -IX(I32)
00656          *      JT = ATOM2                      -IX(I34)
00657          *      KT = ATOM3                      -IX(I36)
00658          *      ICT = ANGLE ARRAY PTR           -IX(I38)
00659          */
00660 #if !defined(USENEWCODE)
00661         for (i=0; i<prm->Ntheta; i++)
00662           fscanf(file, " %d %d %d %d", 
00663                  &prm->AngleAt1[i], &prm->AngleAt2[i], 
00664                  &prm->AngleAt3[i], &prm->AngleNum[i]);
00665 #else
00666         buffer = new int[4*prm->Ntheta];
00667         if (!read_fortran_12I6(file, buffer, 4*prm->Ntheta)) {
00668           return MOLFILE_ERROR;
00669         }
00670         for (i=0; i<prm->Ntheta; i++) { 
00671           prm->AngleAt1[i] = buffer[4*i];
00672           prm->AngleAt2[i] = buffer[4*i+1];
00673           prm->AngleAt3[i] = buffer[4*i+2];
00674           prm->AngleNum[i] = buffer[4*i+3];
00675         }
00676         delete [] buffer;
00677 #endif
00678         readtoeoln(file);
00679 
00680         /* 
00681          * READ DIHEDRAL W/ HYDROGEN: 
00682          *      ITH = ATOM1                     -IX(40)
00683          *      JTH = ATOM2                     -IX(42)
00684          *      KTH = ATOM3                     -IX(44)
00685          *      LTH = ATOM4                     -IX(46)
00686          *      ICTH = DIHEDRAL ARRAY PTR       -IX(48)
00687          */
00688 #if !defined(USENEWCODE)
00689         for (i=0; i<prm->Nphih; i++)
00690           fscanf(file, " %d %d %d %d %d", 
00691                  &prm->DihHAt1[i], &prm->DihHAt2[i], &prm->DihHAt3[i], 
00692                  &prm->DihHAt4[i], &prm->DihHNum[i]);
00693 #else
00694         buffer = new int[5*prm->Nphih];
00695         if (!read_fortran_12I6(file, buffer, 5*prm->Nphih)) {
00696           return MOLFILE_ERROR;
00697         }
00698         for (i=0; i<prm->Nphih; i++) { 
00699           prm->DihHAt1[i] = buffer[5*i];
00700           prm->DihHAt2[i] = buffer[5*i+1];
00701           prm->DihHAt3[i] = buffer[5*i+2];
00702           prm->DihHAt4[i] = buffer[5*i+3];
00703           prm->DihHNum[i] = buffer[5*i+4];
00704         }
00705         delete [] buffer;
00706 #endif
00707         readtoeoln(file);
00708 
00709         /* 
00710          * READ DIHEDRAL W/OUT HYDROGEN: 
00711          *      IT = ATOM1
00712          *      JT = ATOM2
00713          *      KT = ATOM3
00714          *      LT = ATOM4
00715          *      ICT = DIHEDRAL ARRAY PTR
00716          */
00717 #if !defined(USENEWCODE)
00718         for (i=0; i<prm->Nphia; i++) {
00719           fscanf(file, " %d %d %d %d %d", 
00720                  &prm->DihAt1[i], &prm->DihAt2[i], &prm->DihAt3[i], 
00721                  &prm->DihAt4[i], &prm->DihNum[i]);
00722         }
00723 #else
00724         buffer = new int[5*prm->Nphia];
00725         if (!read_fortran_12I6(file, buffer, 5*prm->Nphia)) {
00726           return MOLFILE_ERROR;
00727         }
00728         for (i=0; i<prm->Nphia; i++) {
00729           prm->DihAt1[i] = buffer[5*i];
00730           prm->DihAt2[i] = buffer[5*i+1];
00731           prm->DihAt3[i] = buffer[5*i+2];
00732           prm->DihAt4[i] = buffer[5*i+3];
00733           prm->DihNum[i] = buffer[5*i+4];
00734         }
00735         delete [] buffer;
00736 #endif
00737         readtoeoln(file);
00738 
00739         /*
00740          * READ EXCLUDED ATOM LIST      -IX(I10)
00741          */
00742 #if !defined(USENEWCODE)
00743         for (i=0; i<prm->Nnb; i++)
00744                 fscanf(file, " %d", &prm->ExclAt[i]);
00745 #else
00746         if (!read_fortran_12I6(file, prm->ExclAt, prm->Nnb)) {
00747            return MOLFILE_ERROR;
00748         }
00749 #endif
00750         readtoeoln(file);
00751 
00752         /*
00753          * READ H-BOND R**12 TERM FOR ALL N-B TYPES     -ASOL()
00754          */
00755         for (i=0; i<prm->Nphb; i++) 
00756           fscanf(file, " " DBLFMT, &prm->HB12[i]);
00757         readtoeoln(file);
00758 
00759         /*
00760          * READ H-BOND R**6 TERM FOR ALL N-B TYPES      -BSOL()
00761          */
00762         for (i=0; i<prm->Nphb; i++) 
00763           fscanf(file, " " DBLFMT, &prm->HB6[i]);
00764         readtoeoln(file);
00765 
00766         /*
00767          * READ H-BOND CUTOFF (NOT USED) ??             -HBCUT()
00768          */
00769         H = (_REAL *) get(prm->Nphb * sizeof(_REAL));
00770         for (i=0; i<prm->Nphb; i++) 
00771           fscanf(file, " " DBLFMT, &H[i]);
00772         free((char *)H);
00773         readtoeoln(file);
00774 
00775         /*
00776          * READ ATOM SYMBOLS (FOR ANALYSIS PROGS)       -IH(M06)
00777          */
00778         for (i=0; i<(prm->Natom/20 + (prm->Natom%20 ? 1 : 0)); i++)
00779           preadln(file, &prm->AtomSym[i*80]);
00780 
00781         /*
00782          * READ TREE SYMBOLS (FOR ANALYSIS PROGS)       -IH(M08)
00783          */
00784         for (i=0; i<(prm->Natom/20 + (prm->Natom%20 ? 1 : 0)); i++)
00785           preadln(file, &prm->AtomTree[i*80]);
00786       
00787         /*
00788          * READ TREE JOIN INFO (FOR ANALYSIS PROGS)     -IX(I64)
00789          */
00790 #if !defined(USENEWCODE)
00791         for (i=0; i<prm->Natom; i++)
00792           fscanf(file, " %d", &prm->TreeJoin[i]);
00793 #else
00794         if (!read_fortran_12I6(file, prm->TreeJoin, prm->Natom)) {
00795           return MOLFILE_ERROR;
00796         }
00797 #endif
00798         readtoeoln(file);
00799 
00800         /*
00801          * READ PER-ATOM RES NUMBER                     -IX(I66)
00802          *      NOTE: this appears to be something entirely different
00803          *      NOTE: overwriting this with correct PER-ATOM RES NUMBERs
00804          */
00805 #if !defined(USENEWCODE)
00806         for (i=0; i<prm->Natom; i++)
00807           fscanf(file, " %d", &prm->AtomRes[i]);
00808 #else
00809         if (!read_fortran_12I6(file, prm->AtomRes, prm->Natom)) {
00810           return MOLFILE_ERROR;
00811         }
00812 #endif
00813         res = 0;
00814         for (i=0; i<prm->Natom; i++) {
00815           if (i+1 == prm->Ipres[res+1]) /* atom is 1st of next res */
00816             res++;
00817           prm->AtomRes[i] = res;
00818         }
00819       
00820         /*
00821          * BOUNDARY CONDITION STUFF
00822          */
00823         if (!prm->IfBox) {
00824                 prm->Nspm = 1;
00825                 prm->Boundary = (int *) get(sizeof(int)*prm->Nspm);
00826                 prm->Boundary[0] = prm->Natom;
00827         } else {
00828                 readtoeoln(file);
00829 #if !defined(USENEWCODE)
00830                 fscanf(file, " %d %d %d", &prm->Iptres, &prm->Nspm, 
00831                                                                 &prm->Nspsol);
00832 #else
00833                 buffer = new int[3];
00834                 if (!read_fortran_12I6(file, buffer, 3)) { 
00835                   return MOLFILE_ERROR;
00836                 }
00837                 prm->Iptres = buffer[0];
00838                 prm->Nspm = buffer[1];
00839                 prm->Nspsol = buffer[2];
00840                 delete [] buffer;
00841 #endif
00842                 readtoeoln(file);
00843                 prm->Boundary = (int *) get(sizeof(int)*prm->Nspm);
00844 #if !defined(USENEWCODE)
00845                 for (i=0; i<prm->Nspm; i++)
00846                         fscanf(file, " %d", &prm->Boundary[i]);
00847 #else
00848                 if (!read_fortran_12I6(file, prm->Boundary, prm->Nspm)) {
00849                   return MOLFILE_ERROR;
00850                 }
00851 #endif
00852                 readtoeoln(file);
00853                 fscanf(file, " " DBLFMT " " DBLFMT " " DBLFMT, 
00854                                 &prm->Box[0], &prm->Box[1], &prm->Box[2]);
00855                 readtoeoln(file);
00856                 if (prm->Iptres)
00857                         prm->Ipatm = prm->Ipres[prm->Iptres] - 1; 
00858                 /* IF(IPTRES.GT.0) IPTATM = IX(I02+IPTRES-1+1)-1 */
00859         }
00860 
00861         /*
00862          * ----- LOAD THE CAP INFORMATION IF NEEDED -----
00863          */
00864         if (prm->IfCap) {
00865                 /* if (prm->IfBox) 
00866                         skipeoln(file); */
00867                 fscanf(file, " %d " DBLFMT " " DBLFMT " " DBLFMT " " DBLFMT,
00868                                 &prm->Natcap, &prm->Cutcap, 
00869                                 &prm->Xcap, &prm->Ycap, &prm->Zcap);
00870         }
00871 
00872   return MOLFILE_SUCCESS;
00873 }
00874 
00875 
00876 void ReadPARM::get_parm_atom(int i, char *name, char *atype, char *resname,
00877 char *segname, int *resid, float *q, float *m) {
00878   int nres = prm->Nres;
00879 
00880   int j,k;
00881   int flag = 0;
00882   char *blank = (char *) " ";
00883 
00884   *q = prm->Charges[i];
00885   *m = prm->Masses[i];
00886 
00887   for (k = 0; k < 4; k++) {
00888     if (prm->AtomNames[i*4+k] == *blank)
00889        name[k] = '\0';
00890     else
00891        name[k] = prm->AtomNames[i*4+k];
00892   }
00893   name[k] = '\0';
00894 
00895   for (k = 0; k < 4; k++) {
00896     if ((prm->AtomSym[i*4+k]) == *blank)
00897        atype[k] = '\0';
00898     else
00899        atype[k] = prm->AtomSym[i*4+k];
00900   }
00901   atype[k] = '\0';
00902 
00903   for (j = 0; j < nres-1; j++)
00904     if (((i+1) >= prm->Ipres[j]) && ((i+1) < prm->Ipres[j+1])) {
00905       *resid = j;
00906       resname[0] = prm->ResNames[j*4];
00907       resname[1] = prm->ResNames[j*4+1];
00908       resname[2] = prm->ResNames[j*4+2];
00909       resname[3] = '\0';
00910       flag = 1;
00911     }
00912   if (flag == 0) {
00913      *resid = j;
00914      resname[0] = prm->ResNames[j*4];
00915      resname[1] = prm->ResNames[j*4+1];
00916      resname[2] = prm->ResNames[j*4+2];
00917      resname[3] = '\0';
00918      flag = 1;
00919   }
00920 
00921   segname[0] = '\0'; 
00922 
00923 }
00924 
00925 void ReadPARM::get_parm_bond (int i, int fromAtom[], int toAtom[]) {
00926 
00927   if (i < prm->Nbona) {
00928      fromAtom[i] = (int) ((prm->BondAt1[i])/3 + 1);
00929      toAtom[i] =   (int) ((prm->BondAt2[i])/3 + 1);
00930   }
00931   else get_hydrogen_bond (i, fromAtom, toAtom);
00932 }
00933 
00934 void ReadPARM::get_hydrogen_bond(int i, int fromAtom[], int toAtom[]) {
00935   fromAtom[i] = (int) ((prm->BondHAt1[i-prm->Nbona])/3 + 1);
00936   toAtom[i] = (int) ((prm->BondHAt2[i-prm->Nbona])/3 + 1);
00937 }
00938 
00939 
00940 int ReadPARM::get_parm_natoms() {return prm->Natom;}
00941 
00942 int ReadPARM::get_parm_nbonds() {return (prm->Nbona + prm->Nbonh);}
00943 
00944 int ReadPARM::get_parm_boxInfo() {return(prm->IfBox);}
00945 
00946 
00947 // Read FORTRAN 12I6 format data (no space between adjacent data is assumed)
00948 // One needs to read the whole data block into the buffer here
00949 // fp - file pointer.
00950 // data - buffer to hold the whole block. Should be allocated before the call
00951 // count - number of ints in the block.
00952 int ReadPARM::read_fortran_12I6(FILE *fp, int *data, int count) {
00953   int i, j;
00954   char buf[7];
00955 
00956   for (i=0; i<count; ++i) { 
00957     for (j=0; j<6; ++j) { 
00958       buf[j]=getc(fp);
00959       if (buf[j]=='\n' || buf[j]=='\0' || buf[j]==EOF)
00960         return 0;
00961     }
00962     buf[6] = '\0';
00963 
00964     if (sscanf(buf,"%d",data+i) != 1)
00965       return 0;
00966 
00967     if (i%12==11 && i<count-1)
00968       readtoeoln(fp);
00969   }
00970 
00971   return 1;
00972 }
00973 
00974 
00975 #endif
00976