parm.C

Go to the documentation of this file.
00001 
00002 /*
00003  * COPYRIGHT 1992, REGENTS OF THE UNIVERSITY OF CALIFORNIA
00004  *
00005  *  prm.c - read information from an amber PARM topology file: 
00006  *      atom/residue/bond/charge info, plus force field data. 
00007  *      This file and the accompanying prm.h may be distributed 
00008  *      provided this notice is retained unmodified and provided 
00009  *      that any modifications to the rest of the file are noted 
00010  *      in comments.
00011  *
00012  *      Bill Ross, UCSF 1994
00013  */
00014 
00015 // The original implementation requires space between adjacent
00016 // data. This is not true for AMBER 6 format files if it's a
00017 // large system. So a large part of the code is modified. All
00018 // reads of integer in 12I6 format is now implemented using the
00019 // new method, which doesn't assume space between entries.
00020 // Also, code has been added to read AMBER 7 format files.
00021 // Apr. 9, 2002
00022 
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <ctype.h>
00026 #include <sys/types.h>
00027 #include <sys/stat.h>
00028 #include <errno.h>
00029 #include <string.h>
00030 #include "strlib.h"             //  For strcasecmp and strncasecmp
00031 
00032 #include "common.h"
00033 #include "InfoStream.h"
00034 #include "parm.h"
00035 
00036 static int      debug = 0;      /* set it if you want */
00037 
00038 /*      fortran formats 
00039  *       9118 FORMAT(12I6)
00040  *       9128 FORMAT(5E16.8)
00041  */
00042 
00043 // This trick is misleading and still needs to have space between
00044 // data thus won't work for 6 digit integers. Don't use it.
00045 //const char    *f9118 = "%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d\n";
00046 
00047 
00048 // This function was borrowed from VMD code in "ReadPARM.C". Here it
00049 // is used to replace the "skipeoln()" in the original code, which
00050 // seems to be undefined.
00051 static int readtoeoln(FILE *f) {
00052   char c;
00053 
00054   /* skip to eoln */
00055   while((c = getc(f)) != '\n') {
00056     if (c == EOF) 
00057       return -1;
00058   }
00059 
00060   return 0;
00061 }  
00062 
00063 /***********************************************************************
00064                                                         GENOPEN()
00065 ************************************************************************/
00066 
00067 /*
00068  *  genopen() - fopen regular or popen compressed file for reading
00069  */
00070 //  replaced with NAMD implementation
00071 
00072 FILE *Ambertoppar::genopen(const char *name)
00073 {
00074         return(Fopen(name,"r"));
00075 }
00076 
00077 /***********************************************************************
00078                                                         GENCLOSE()
00079 ************************************************************************/
00080 
00081 /*
00082  *  genclose() - close fopened or popened file
00083  */
00084 //  replaced with NAMD implementation
00085 
00086 void Ambertoppar::genclose(FILE *fileptr)
00087 {
00088         Fclose(fileptr);
00089 }
00090 
00091 
00092 /***********************************************************************
00093                                                         GET()
00094 ************************************************************************/
00095 
00096 char *Ambertoppar::get(int size)
00097 {
00098         char    *ptr;
00099 
00100 #ifdef DEBUG
00101         printf("malloc %d\n", size);
00102 #ifndef NAMD_NO_STDOUT_FLUSH
00103         fflush(stdout);
00104 #endif
00105 #endif
00106         if (size ==0)
00107                 return((char *) NULL);
00108 
00109         if ((ptr = (char *) malloc((unsigned)size)) == NULL) {
00110                 printf("malloc %d", size);
00111 #ifndef NAMD_NO_STDOUT_FLUSH
00112                 fflush(stdout);
00113 #endif
00114                 NAMD_die("Memory allocation error in Ambertoppar::get()");
00115         }
00116         return(ptr);
00117 }
00118 
00119 /***********************************************************************
00120                                                         PREADLN()
00121 ************************************************************************/
00122 
00123 void Ambertoppar::preadln(FILE *file, const char *name, char *string)
00124 {
00125         int     i, j;
00126 
00127         for (i=0; i<81; i++) {
00128                 if ((j = getc(file)) == EOF) {
00129                         NAMD_die("Unexpected EOF in Amber parm file");
00130                         printf("Error: unexpected EOF in %s\n", name);
00131                 }
00132                 string[i] = (char) j;
00133                 if (string[i] == '\n') {
00134                         break;
00135                 }
00136         }
00137         if (i == 80  &&  string[i] != '\n') {
00138                 NAMD_die("Line too long in Amber parm file");
00139                 printf("Error: line too long in %s:\n%.80s", name, string);
00140         }
00141 }
00142 
00143 /***************************************************************************
00144                                                                 READPARM()
00145 ****************************************************************************/
00146 
00147 /*
00148  * readparm() - instantiate a given Ambertoppar
00149  */
00150 
00151 int Ambertoppar::readparm(char *name)
00152 {
00153         _REAL           *H;
00154         int             i, idum, res, ifpert;
00155         int             *buffer, amber7_format;
00156         FILE            *file;
00157 
00158         if (data_read)
00159         { iout << "Duplicate parm data in one object!\n" << endi;
00160           return(0);
00161         }
00162 
00163 //      printf("Reading parm file (%s)\n", name);
00164         iout << "Reading parm file (" << name << ") ...\n" << endi;
00165 
00166 //      if ((file = genopen(name, "parm")) == NULL) 
00167         if ((file = genopen(name)) == NULL)
00168                 return(0);
00169 
00170         /* READ TITLE */
00171 
00172         preadln(file, name, ititl);
00173 // "ititle" doesn't guarantee to have '\0' (as the end of a string),
00174 // so the following is disabled in order to avoid strange output
00175 //      printf("%s title:\n%s", name, ititl);
00176 
00177         // Check whether it's in AMBER 7 format or in old format
00178         
00179         if (strncmp(ititl,"%VERSION",8))
00180                 amber7_format = 0;  // old format
00181         else
00182         {       amber7_format = 1;  // AMBER 7 format
00183                 iout << "PARM file in AMBER 7 format\n" << endi;
00184                 if (!moveto(file,"TITLE"))
00185                 { genclose(file);
00186                   return 0;
00187                 }
00188                 preadln(file, name, ititl);
00189         }
00190 
00191         /* READ CONTROL INTEGERS */
00192 
00193         if (amber7_format)
00194         {       if (!moveto(file,"POINTERS"))
00195                 { genclose(file);
00196                   return 0;
00197                 }
00198 //      fscanf(file, f9118, 
00199                 fscanf(file, "%d%d%d%d%d%d%d%d%d%d%d%d",
00200                 &Natom,  &Ntypes, &Nbonh, &Mbona, 
00201                 &Ntheth, &Mtheta, &Nphih, &Mphia, 
00202                 &Nhparm, &Nparm,  &Nnb,   &Nres);
00203 
00204 //      fscanf(file, f9118, 
00205                 fscanf(file, "%d%d%d%d%d%d%d%d%d%d%d%d",
00206                 &Nbona,  &Ntheta, &Nphia, &Numbnd, 
00207                 &Numang, &Nptra,  &Natyp, &Nphb, 
00208                 &ifpert,      &idum,        &idum,       &idum);
00209 
00210                 fscanf(file, " %d %d %d %d %d %d", 
00211                 &idum, &idum,&idum,&IfBox,&Nmxrs,&IfCap);
00212         }
00213         else
00214         {       buffer = new int[30];
00215                 if (!read_fortran_12I6(file,buffer,30))
00216                 { genclose(file);
00217                   return 0;
00218                 }
00219                 Natom = buffer[0];
00220                 Ntypes = buffer[1];
00221                 Nbonh = buffer[2];
00222                 Mbona = buffer[3];
00223                 Ntheth = buffer[4];
00224                 Mtheta = buffer[5];
00225                 Nphih = buffer[6];
00226                 Mphia = buffer[7];
00227                 Nhparm = buffer[8];
00228                 Nparm = buffer[9];
00229                 Nnb = buffer[10];
00230                 Nres = buffer[11];
00231                 Nbona = buffer[12];
00232                 Ntheta = buffer[13];
00233                 Nphia = buffer[14];
00234                 Numbnd = buffer[15];
00235                 Numang = buffer[16];
00236                 Nptra = buffer[17];
00237                 Natyp = buffer[18];
00238                 Nphb = buffer[19];
00239                 ifpert = buffer[20];
00240                 IfBox = buffer[27];
00241                 Nmxrs = buffer[28];
00242                 IfCap = buffer[29];
00243                 delete [] buffer;
00244 //      skipeoln(file);
00245                 readtoeoln(file);
00246         }
00247 
00248         if (ifpert) {
00249                 printf("not equipped to read perturbation prmtop\n");
00250                 return(0);
00251         }
00252 
00253         /* ALLOCATE MEMORY */
00254 
00255         Nat3 = 3 * Natom;
00256         Ntype2d = Ntypes * Ntypes;
00257         Nttyp = Ntypes*(Ntypes+1)/2;
00258 
00259         /*
00260          * get most of the indirect stuff; some extra allowed for char arrays
00261          */
00262 
00263         AtomNames = (char *) get(4*Natom+81);
00264         Charges = (_REAL *) get(sizeof(_REAL)*Natom);
00265         Masses = (_REAL *) get(sizeof(_REAL)*Natom);
00266         Iac = (int *) get(sizeof(int)*Natom);
00267         Iblo = (int *) get(sizeof(int)*Natom);
00268         Cno = (int *) get(sizeof(int)* Ntype2d);
00269         ResNames = (char *) get(4* Nres+81);
00270         Ipres = (int *) get(sizeof(int)*( Nres+1));
00271         Rk = (_REAL *) get(sizeof(_REAL)* Numbnd);
00272         Req = (_REAL *) get(sizeof(_REAL)* Numbnd);
00273         Tk = (_REAL *) get(sizeof(_REAL)* Numang);
00274         Teq = (_REAL *) get(sizeof(_REAL)* Numang);
00275         Pk = (_REAL *) get(sizeof(_REAL)* Nptra);
00276         Pn = (_REAL *) get(sizeof(_REAL)* Nptra);
00277         Phase = (_REAL *) get(sizeof(_REAL)* Nptra);
00278         Solty = (_REAL *) get(sizeof(_REAL)* Natyp);
00279         Cn1 = (_REAL *) get(sizeof(_REAL)* Nttyp);
00280         Cn2 = (_REAL *) get(sizeof(_REAL)* Nttyp);
00281         BondHAt1 = (int *) get(sizeof(int)* Nbonh);
00282         BondHAt2 = (int *) get(sizeof(int)* Nbonh);
00283         BondHNum = (int *) get(sizeof(int)* Nbonh);
00284         BondAt1 = (int *) get(sizeof(int)* Nbona);
00285         BondAt2 = (int *) get(sizeof(int)* Nbona);
00286         BondNum = (int *) get(sizeof(int)* Nbona);
00287         AngleHAt1 = (int *) get(sizeof(int)* Ntheth);
00288         AngleHAt2 = (int *) get(sizeof(int)* Ntheth);
00289         AngleHAt3 = (int *) get(sizeof(int)* Ntheth);
00290         AngleHNum = (int *) get(sizeof(int)* Ntheth);
00291         AngleAt1 = (int *) get(sizeof(int)* Ntheta);
00292         AngleAt2 = (int *) get(sizeof(int)*Ntheta);
00293         AngleAt3 = (int *) get(sizeof(int)*Ntheta);
00294         AngleNum = (int *) get(sizeof(int)*Ntheta);
00295         DihHAt1 = (int *) get(sizeof(int)*Nphih);
00296         DihHAt2 = (int *) get(sizeof(int)*Nphih);
00297         DihHAt3 = (int *) get(sizeof(int)*Nphih);
00298         DihHAt4 = (int *) get(sizeof(int)*Nphih);
00299         DihHNum = (int *) get(sizeof(int)*Nphih);
00300         DihAt1 = (int *) get(sizeof(int)*Nphia);
00301         DihAt2 = (int *) get(sizeof(int)*Nphia);
00302         DihAt3 = (int *) get(sizeof(int)*Nphia);
00303         DihAt4 = (int *) get(sizeof(int)*Nphia);
00304         DihNum = (int *) get(sizeof(int)*Nphia);
00305         ExclAt = (int *) get(sizeof(int)*Nnb);
00306         HB12 = (_REAL *) get(sizeof(_REAL)*Nphb);
00307         HB6 = (_REAL *) get(sizeof(_REAL)*Nphb);
00308         AtomSym = (char *) get(4*Natom+81);
00309         AtomTree = (char *) get(4*Natom+81);
00310         TreeJoin = (int *) get(sizeof(int)*Natom);
00311         AtomRes = (int *) get(sizeof(int)*Natom);
00312 
00313         /* 
00314          * READ ATOM NAMES -IH(M04)
00315          */
00316 
00317         if (amber7_format)
00318                 if (!moveto(file,"ATOM_NAME"))
00319                 { genclose(file);
00320                   return 0;
00321                 }
00322         for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
00323                 preadln(file, "", &AtomNames[i*80]);
00324 
00325         /* 
00326          * READ ATOM CHARGES -X(L15)
00327          *      (pre-multiplied by an energy factor of 18.2223 == sqrt(332)
00328          *       for faster force field calculations)
00329          */
00330 
00331         if (amber7_format)
00332                 if (!moveto(file,"CHARGE"))
00333                 { genclose(file);
00334                   return 0;
00335                 }
00336         for (i=0; i<Natom; i++)
00337 #ifdef DOUBLE
00338                 fscanf(file, " %lf", &Charges[i]);
00339 #else
00340                 fscanf(file, " %f", &Charges[i]);
00341 #endif
00342 //      skipeoln(file);
00343         readtoeoln(file);
00344 
00345         /* 
00346          * READ ATOM MASSES -X(L20)
00347          */
00348 
00349         if (amber7_format)
00350                 if (!moveto(file,"MASS"))
00351                 { genclose(file);
00352                   return 0;
00353                 }
00354         for (i=0; i<Natom; i++)
00355 #ifdef DOUBLE
00356                 fscanf(file, " %le", &Masses[i]);
00357 #else
00358                 fscanf(file, " %e", &Masses[i]);
00359 #endif
00360 //      skipeoln(file);
00361         readtoeoln(file);
00362 
00363         /* 
00364          * READ ATOM L-J TYPES -IX(I04)
00365          */
00366 
00367         if (amber7_format)
00368         {       if (!moveto(file,"ATOM_TYPE_INDEX"))
00369                 { genclose(file);
00370                   return 0;
00371                 }
00372                 for (i=0; i<Natom; i++)
00373                         fscanf(file, " %d", &Iac[i]);
00374         }
00375         else
00376         {       if (!read_fortran_12I6(file,Iac,Natom))
00377                 { genclose(file);
00378                   return 0;
00379                 }
00380 //      skipeoln(file);
00381                 readtoeoln(file);
00382         }
00383 
00384         /* 
00385          * READ ATOM INDEX TO 1st IN EXCLUDED ATOM LIST "NATEX" -IX(I08)
00386          */
00387 
00388         if (amber7_format)
00389         {       if (!moveto(file,"NUMBER_EXCLUDED_ATOMS"))
00390                 { genclose(file);
00391                   return 0;
00392                 }
00393                 for (i=0; i<Natom; i++)
00394                         fscanf(file, " %d", &Iblo[i]);
00395         }
00396         else
00397         {       if (!read_fortran_12I6(file,Iblo,Natom))
00398                 { genclose(file);
00399                   return 0;
00400                 }
00401 //      skipeoln(file);
00402                 readtoeoln(file);
00403         }
00404 
00405         /* 
00406          * READ TYPE INDEX TO N-B TYPE -IX(I06)
00407          */
00408 
00409         if (amber7_format)
00410         {       if (!moveto(file,"NONBONDED_PARM_INDEX"))
00411                 { genclose(file);
00412                   return 0;
00413                 }
00414                 for (i=0; i<Ntype2d; i++)
00415                         fscanf(file, " %d", &Cno[i]);
00416         }
00417         else
00418         {       if (!read_fortran_12I6(file,Cno,Ntype2d))
00419                 { genclose(file);
00420                   return 0;
00421                 }
00422 //      skipeoln(file);
00423                 readtoeoln(file);
00424         }
00425 
00426         /* 
00427          * READ RES NAMES (4 chars each, 4th blank) -IH(M02)
00428          */
00429 
00430         if (amber7_format)
00431                 if (!moveto(file,"RESIDUE_LABEL"))
00432                 { genclose(file);
00433                   return 0;
00434                 }
00435         for (i=0; i<(Nres/20 + (Nres%20 ? 1 : 0)); i++)
00436                 preadln(file, "", &ResNames[i*80]);
00437 
00438         /* 
00439          * READ RES POINTERS TO 1st ATOM                -IX(I02)
00440          */
00441 
00442         if (amber7_format)
00443         {       if (!moveto(file,"RESIDUE_POINTER"))
00444                 { genclose(file);
00445                   return 0;
00446                 }
00447                 for (i=0; i<Nres; i++) 
00448                         fscanf(file, " %d", &Ipres[i]);
00449         }
00450         else
00451         {       if (!read_fortran_12I6(file,Ipres,Nres))
00452                 { genclose(file);
00453                   return 0;
00454                 }
00455 //      skipeoln(file);
00456                 readtoeoln(file);
00457         }
00458         Ipres[Nres] = Natom + 1;
00459 
00460         /* 
00461          * READ BOND FORCE CONSTANTS                    -RK()
00462          */
00463 
00464         if (amber7_format)
00465                 if (!moveto(file,"BOND_FORCE_CONSTANT"))
00466                 { genclose(file);
00467                   return 0;
00468                 }
00469         for (i=0; i< Numbnd; i++) 
00470 #ifdef DOUBLE
00471                 fscanf(file, " %lf", &Rk[i]);
00472 #else
00473                 fscanf(file, " %f", &Rk[i]);
00474 #endif
00475 //      skipeoln(file);
00476         readtoeoln(file);
00477 
00478         /* 
00479          * READ BOND LENGTH OF MINIMUM ENERGY           -REQ()
00480          */
00481 
00482         if (amber7_format)
00483                 if (!moveto(file,"BOND_EQUIL_VALUE"))
00484                 { genclose(file);
00485                   return 0;
00486                 }
00487         for (i=0; i< Numbnd; i++) 
00488 #ifdef DOUBLE
00489                 fscanf(file, " %lf", &Req[i]);
00490 #else
00491                 fscanf(file, " %f", &Req[i]);
00492 #endif
00493 //      skipeoln(file);
00494         readtoeoln(file);
00495 
00496         /* 
00497          * READ BOND ANGLE FORCE CONSTANTS (following Rk nomen) -TK()
00498          */
00499 
00500         if (amber7_format)
00501                 if (!moveto(file,"ANGLE_FORCE_CONSTANT"))
00502                 { genclose(file);
00503                   return 0;
00504                 }
00505         for (i=0; i< Numang; i++) 
00506 #ifdef DOUBLE
00507                 fscanf(file, " %lf", &Tk[i]);
00508 #else
00509                 fscanf(file, " %f", &Tk[i]);
00510 #endif
00511 //      skipeoln(file);
00512         readtoeoln(file);
00513 
00514         /* 
00515          * READ BOND ANGLE OF MINIMUM ENERGY (following Req nomen) -TEQ()
00516          */
00517 
00518         if (amber7_format)
00519                 if (!moveto(file,"ANGLE_EQUIL_VALUE"))
00520                 { genclose(file);
00521                   return 0;
00522                 }
00523         for (i=0; i< Numang; i++) 
00524 #ifdef DOUBLE
00525                 fscanf(file, " %lf", &Teq[i]);
00526 #else
00527                 fscanf(file, " %f", &Teq[i]);
00528 #endif
00529 //      skipeoln(file);
00530         readtoeoln(file);
00531 
00532         /* 
00533          * READ DIHEDRAL PEAK MAGNITUDE                 -PK()
00534          */
00535 
00536         if (amber7_format)
00537                 if (!moveto(file,"DIHEDRAL_FORCE_CONSTANT"))
00538                 { genclose(file);
00539                   return 0;
00540                 }
00541         for (i=0; i< Nptra; i++) 
00542 #ifdef DOUBLE
00543                 fscanf(file, " %lf", &Pk[i]);
00544 #else
00545                 fscanf(file, " %f", &Pk[i]);
00546 #endif
00547 //      skipeoln(file);
00548         readtoeoln(file);
00549 
00550         /* 
00551          * READ DIHEDRAL PERIODICITY                    -PN()
00552          */
00553 
00554         if (amber7_format)
00555                 if (!moveto(file,"DIHEDRAL_PERIODICITY"))
00556                 { genclose(file);
00557                   return 0;
00558                 }
00559         for (i=0; i< Nptra; i++) 
00560 #ifdef DOUBLE
00561                 fscanf(file, " %lf", &Pn[i]);
00562 #else
00563                 fscanf(file, " %f", &Pn[i]);
00564 #endif
00565 //      skipeoln(file);
00566         readtoeoln(file);
00567 
00568         /* 
00569          * READ DIHEDRAL PHASE                          -PHASE()
00570          */
00571 
00572         if (amber7_format)
00573                 if (!moveto(file,"DIHEDRAL_PHASE"))
00574                 { genclose(file);
00575                   return 0;
00576                 }
00577         for (i=0; i< Nptra; i++) 
00578 #ifdef DOUBLE
00579                 fscanf(file, " %lf", &Phase[i]);
00580 #else
00581                 fscanf(file, " %f", &Phase[i]);
00582 #endif
00583 //      skipeoln(file);
00584         readtoeoln(file);
00585 
00586         /* 
00587          * ?? "RESERVED"                                -SOLTY()
00588          */
00589 
00590         if (amber7_format)
00591                 if (!moveto(file,"SOLTY"))
00592                 { genclose(file);
00593                   return 0;
00594                 }
00595         for (i=0; i< Natyp; i++) 
00596 #ifdef DOUBLE
00597                 fscanf(file, " %lf", &Solty[i]);
00598 #else
00599                 fscanf(file, " %f", &Solty[i]);
00600 #endif
00601 //      skipeoln(file);
00602         readtoeoln(file);
00603 
00604         /* 
00605          * READ L-J R**12 FOR ALL PAIRS OF ATOM TYPES   -CN1()
00606          *      (SHOULD BE 0 WHERE H-BONDS)
00607          */
00608 
00609         if (amber7_format)
00610                 if (!moveto(file,"LENNARD_JONES_ACOEF"))
00611                 { genclose(file);
00612                   return 0;
00613                 }
00614         for (i=0; i< Nttyp; i++) 
00615 #ifdef DOUBLE
00616                 fscanf(file, " %lf", &Cn1[i]);
00617 #else
00618                 fscanf(file, " %f", &Cn1[i]);
00619 #endif
00620 //      skipeoln(file);
00621         readtoeoln(file);
00622 
00623         /* 
00624          * READ L-J R**6 FOR ALL PAIRS OF ATOM TYPES    -CN2()
00625          *      (SHOULD BE 0 WHERE H-BONDS)
00626          */
00627 
00628         if (amber7_format)
00629                 if (!moveto(file,"LENNARD_JONES_BCOEF"))
00630                 { genclose(file);
00631                   return 0;
00632                 }
00633         for (i=0; i< Nttyp; i++) 
00634 #ifdef DOUBLE
00635                 fscanf(file, " %lf", &Cn2[i]);
00636 #else
00637                 fscanf(file, " %f", &Cn2[i]);
00638 #endif
00639 //      skipeoln(file);
00640         readtoeoln(file);
00641 
00642         /* 
00643          * READ COVALENT BOND W/ HYDROGEN (3*(atnum-1)): 
00644          *      IBH = ATOM1             -IX(I12)
00645          *      JBH = ATOM2             -IX(I14)
00646          *      ICBH = BOND ARRAY PTR   -IX(I16)
00647          */
00648 
00649         if (amber7_format)
00650         {       if (!moveto(file,"BONDS_INC_HYDROGEN"))
00651                 { genclose(file);
00652                   return 0;
00653                 }
00654                 for (i=0; i<Nbonh; i++)
00655                         fscanf(file, " %d %d %d", 
00656                             &BondHAt1[i], &BondHAt2[i], &BondHNum[i]);
00657         }
00658         else
00659         {       buffer = new int[3*Nbonh];
00660                 if (!read_fortran_12I6(file,buffer,3*Nbonh))
00661                 { genclose(file);
00662                   return 0;
00663                 }
00664                 for (i=0; i<Nbonh; i++) 
00665                 { BondHAt1[i] = buffer[3*i];
00666                   BondHAt2[i] = buffer[3*i+1];
00667                   BondHNum[i] = buffer[3*i+2];
00668                 }
00669                 delete [] buffer;
00670 //      skipeoln(file);
00671                 readtoeoln(file);
00672         }
00673 
00674         /* 
00675          * READ COVALENT BOND W/OUT HYDROGEN (3*(atnum-1)):
00676          *      IB = ATOM1              -IX(I18)
00677          *      JB = ATOM2              -IX(I20)
00678          *      ICB = BOND ARRAY PTR    -IX(I22)
00679          */
00680 
00681         if (amber7_format)
00682         {       if (!moveto(file,"BONDS_WITHOUT_HYDROGEN"))
00683                 { genclose(file);
00684                   return 0;
00685                 }
00686                 for (i=0; i<Nbona; i++)
00687                         fscanf(file, " %d %d %d", 
00688                                 &BondAt1[i], &BondAt2[i], &BondNum[i]);
00689         }
00690         else
00691         {       buffer = new int[3*Nbona];
00692                 if (!read_fortran_12I6(file,buffer,3*Nbona))
00693                 { genclose(file);
00694                   return 0;
00695                 }
00696                 for (i=0; i<Nbona; i++)
00697                 { BondAt1[i] = buffer[3*i];
00698                   BondAt2[i] = buffer[3*i+1];
00699                   BondNum[i] = buffer[3*i+2];
00700                 }
00701                 delete [] buffer;
00702 //      skipeoln(file);
00703                 readtoeoln(file);
00704         }
00705 
00706         /* 
00707          * READ ANGLE W/ HYDROGEN: 
00708          *      ITH = ATOM1                     -IX(I24)
00709          *      JTH = ATOM2                     -IX(I26)
00710          *      KTH = ATOM3                     -IX(I28)
00711          *      ICTH = ANGLE ARRAY PTR          -IX(I30)
00712          */
00713 
00714         if (amber7_format)
00715         {       if (!moveto(file,"ANGLES_INC_HYDROGEN"))
00716                 { genclose(file);
00717                   return 0;
00718                 }
00719                 for (i=0; i<Ntheth; i++)
00720                         fscanf(file, " %d %d %d %d", 
00721                                 &AngleHAt1[i], &AngleHAt2[i], 
00722                                 &AngleHAt3[i], &AngleHNum[i]);
00723         }
00724         else
00725         {       buffer = new int[4*Ntheth];
00726                 if (!read_fortran_12I6(file,buffer,4*Ntheth))
00727                 { genclose(file);
00728                   return 0;
00729                 }
00730                 for (i=0; i<Ntheth; i++)
00731                 { AngleHAt1[i] = buffer[4*i];
00732                   AngleHAt2[i] = buffer[4*i+1];
00733                   AngleHAt3[i] = buffer[4*i+2];
00734                   AngleHNum[i] = buffer[4*i+3];
00735                 }
00736                 delete [] buffer;
00737 //      skipeoln(file);
00738                 readtoeoln(file);
00739         }
00740 
00741         /* 
00742          * READ ANGLE W/OUT HYDROGEN: 
00743          *      IT = ATOM1                      -IX(I32)
00744          *      JT = ATOM2                      -IX(I34)
00745          *      KT = ATOM3                      -IX(I36)
00746          *      ICT = ANGLE ARRAY PTR           -IX(I38)
00747          */
00748 
00749         if (amber7_format)
00750         {       if (!moveto(file,"ANGLES_WITHOUT_HYDROGEN"))
00751                 { genclose(file);
00752                   return 0;
00753                 }
00754                 for (i=0; i<Ntheta; i++)
00755                         fscanf(file, " %d %d %d %d", 
00756                                 &AngleAt1[i], &AngleAt2[i], 
00757                                 &AngleAt3[i], &AngleNum[i]);
00758         }
00759         else
00760         {       buffer = new int[4*Ntheta];
00761                 if (!read_fortran_12I6(file,buffer,4*Ntheta))
00762                 { genclose(file);
00763                   return 0;
00764                 }
00765                 for (i=0; i<Ntheta; i++)
00766                 { AngleAt1[i] = buffer[4*i];
00767                   AngleAt2[i] = buffer[4*i+1];
00768                   AngleAt3[i] = buffer[4*i+2];
00769                   AngleNum[i] = buffer[4*i+3];
00770                 }
00771                 delete [] buffer;
00772 //      skipeoln(file);
00773                 readtoeoln(file);
00774         }
00775 
00776         /* 
00777          * READ DIHEDRAL W/ HYDROGEN: 
00778          *      ITH = ATOM1                     -IX(40)
00779          *      JTH = ATOM2                     -IX(42)
00780          *      KTH = ATOM3                     -IX(44)
00781          *      LTH = ATOM4                     -IX(46)
00782          *      ICTH = DIHEDRAL ARRAY PTR       -IX(48)
00783          */
00784 
00785         if (amber7_format)
00786         {       if (!moveto(file,"DIHEDRALS_INC_HYDROGEN"))
00787                 { genclose(file);
00788                   return 0;
00789                 }
00790                 for (i=0; i<Nphih; i++)
00791                         fscanf(file, " %d %d %d %d %d", 
00792                                 &DihHAt1[i], &DihHAt2[i], &DihHAt3[i], 
00793                                 &DihHAt4[i], &DihHNum[i]);
00794         }
00795         else
00796         {       buffer = new int[5*Nphih];
00797                 if (!read_fortran_12I6(file,buffer,5*Nphih))
00798                 { genclose(file);
00799                   return 0;
00800                 }
00801                 for (i=0; i<Nphih; i++)
00802                 { DihHAt1[i] = buffer[5*i];
00803                   DihHAt2[i] = buffer[5*i+1];
00804                   DihHAt3[i] = buffer[5*i+2];
00805                   DihHAt4[i] = buffer[5*i+3];
00806                   DihHNum[i] = buffer[5*i+4];
00807                 }
00808                 delete [] buffer;
00809 //      skipeoln(file);
00810                 readtoeoln(file);
00811         }
00812 
00813         /* 
00814          * READ DIHEDRAL W/OUT HYDROGEN: 
00815          *      IT = ATOM1
00816          *      JT = ATOM2
00817          *      KT = ATOM3
00818          *      LT = ATOM4
00819          *      ICT = DIHEDRAL ARRAY PTR
00820          */
00821 
00822         if (amber7_format)
00823         {       if (!moveto(file,"DIHEDRALS_WITHOUT_HYDROGEN"))
00824                 { genclose(file);
00825                   return 0;
00826                 }
00827                 for (i=0; i<Nphia; i++)
00828                         fscanf(file, " %d %d %d %d %d", 
00829                                 &DihAt1[i], &DihAt2[i], &DihAt3[i], 
00830                                 &DihAt4[i], &DihNum[i]);
00831         }
00832         else
00833         {       buffer = new int[5*Nphia];
00834                 if (!read_fortran_12I6(file,buffer,5*Nphia))
00835                 { genclose(file);
00836                   return 0;
00837                 }
00838                 for (i=0; i<Nphia; i++) {
00839                   DihAt1[i] = buffer[5*i];
00840                   DihAt2[i] = buffer[5*i+1];
00841                   DihAt3[i] = buffer[5*i+2];
00842                   DihAt4[i] = buffer[5*i+3];
00843                   DihNum[i] = buffer[5*i+4];
00844                 }
00845                 delete [] buffer;
00846 //      skipeoln(file);
00847                 readtoeoln(file);
00848         }
00849 
00850         /*
00851          * READ EXCLUDED ATOM LIST      -IX(I10)
00852          */
00853 
00854         if (amber7_format)
00855         {       if (!moveto(file,"EXCLUDED_ATOMS_LIST"))
00856                 { genclose(file);
00857                   return 0;
00858                 }
00859                 for (i=0; i<Nnb; i++)
00860                         fscanf(file, " %d", &ExclAt[i]);
00861         }
00862         else
00863         {       if (!read_fortran_12I6(file,ExclAt,Nnb))
00864                 { genclose(file);
00865                   return 0;
00866                 }
00867 //      skipeoln(file);
00868                 readtoeoln(file);
00869         }
00870 
00871         /*
00872          * READ H-BOND R**12 TERM FOR ALL N-B TYPES     -ASOL()
00873          */
00874 
00875         if (amber7_format)
00876                 if (!moveto(file,"HBOND_ACOEF"))
00877                 { genclose(file);
00878                   return 0;
00879                 }
00880         for (i=0; i<Nphb; i++) 
00881 #ifdef DOUBLE
00882                 fscanf(file, " %lf", &HB12[i]);
00883 #else
00884                 fscanf(file, " %f", &HB12[i]);
00885 #endif
00886 //      skipeoln(file);
00887         readtoeoln(file);
00888 
00889         /*
00890          * READ H-BOND R**6 TERM FOR ALL N-B TYPES      -BSOL()
00891          */
00892 
00893         if (amber7_format)
00894                 if (!moveto(file,"HBOND_BCOEF"))
00895                 { genclose(file);
00896                   return 0;
00897                 }
00898         for (i=0; i<Nphb; i++) 
00899 #ifdef DOUBLE
00900                 fscanf(file, " %lf", &HB6[i]);
00901 #else
00902                 fscanf(file, " %f", &HB6[i]);
00903 #endif
00904 //      skipeoln(file);
00905         readtoeoln(file);
00906 
00907         /*
00908          * READ H-BOND CUTOFF (NOT USED) ??             -HBCUT()
00909          */
00910 
00911         if (amber7_format)
00912                 if (!moveto(file,"HBCUT"))
00913                 { genclose(file);
00914                   return 0;
00915                 }
00916         H = (_REAL *) get(Nphb * sizeof(_REAL));
00917         for (i=0; i<Nphb; i++) 
00918 #ifdef DOUBLE
00919                 fscanf(file, " %lf", &H[i]);
00920 #else
00921                 fscanf(file, " %f", &H[i]);
00922 #endif
00923         free((char *)H);
00924         
00925 //      skipeoln(file);
00926         readtoeoln(file);
00927 
00928         /*
00929          * READ ATOM SYMBOLS (FOR ANALYSIS PROGS)       -IH(M06)
00930          */
00931 
00932         if (amber7_format)
00933                 if (!moveto(file,"AMBER_ATOM_TYPE"))
00934                 { genclose(file);
00935                   return 0;
00936                 }
00937         for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
00938                 preadln(file, "", &AtomSym[i*80]);
00939 
00940         /*
00941          * READ TREE SYMBOLS (FOR ANALYSIS PROGS)       -IH(M08)
00942          */
00943 
00944         if (amber7_format)
00945                 if (!moveto(file,"TREE_CHAIN_CLASSIFICATION"))
00946                 { genclose(file);
00947                   return 0;
00948                 }
00949         for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
00950                 preadln(file, "", &AtomTree[i*80]);
00951 
00952         /*
00953          * READ TREE JOIN INFO (FOR ANALYSIS PROGS)     -IX(I64)
00954          */
00955 
00956         if (amber7_format)
00957         {       if (!moveto(file,"JOIN_ARRAY"))
00958                 { genclose(file);
00959                   return 0;
00960                 }
00961                 for (i=0; i<Natom; i++)
00962                         fscanf(file, " %d", &TreeJoin[i]);
00963         }
00964         else
00965         {       if (!read_fortran_12I6(file,TreeJoin,Natom))
00966                 { genclose(file);
00967                   return 0;
00968                 }
00969 //      skipeoln(file);
00970                 readtoeoln(file);
00971         }
00972 
00973         /*
00974          * READ PER-ATOM RES NUMBER                     -IX(I66)
00975          *      NOTE: this appears to be something entirely different
00976          *      NOTE: overwriting this with correct PER-ATOM RES NUMBERs
00977          */
00978 
00979         if (amber7_format)
00980         {       if (!moveto(file,"IROTAT"))
00981                 { genclose(file);
00982                   return 0;
00983                 }
00984                 for (i=0; i<Natom; i++)
00985                         fscanf(file, " %d", &AtomRes[i]);
00986         }
00987         else
00988                 if (!read_fortran_12I6(file,AtomRes,Natom))
00989                 { genclose(file);
00990                   return 0;
00991                 }
00992         res = 0;
00993         for (i=0; i<Natom; i++) {
00994                 if (i+1 == Ipres[res+1])        /* atom is 1st of next res */
00995                         res++;
00996                 AtomRes[i] = res;
00997         }
00998       
00999         /*
01000          * BOUNDARY CONDITION STUFF
01001          */
01002 
01003         if (!IfBox) {
01004                 Nspm = 1;
01005                 Boundary = (int *) get(sizeof(int)*Nspm);
01006                 Boundary[0] = Natom;
01007         } else {
01008                 if (amber7_format)
01009                 {       if (!moveto(file,"SOLVENT_POINTERS"))
01010                         { genclose(file);
01011                           return 0;
01012                         }
01013                         fscanf(file, " %d %d %d", &Iptres, &Nspm, 
01014                                                         &Nspsol);
01015                 }
01016                 else
01017                 {
01018 //              skipeoln(file);
01019                         readtoeoln(file);
01020                         buffer = new int[3];
01021                         if (!read_fortran_12I6(file,buffer,3))
01022                         { genclose(file);
01023                           return 0;
01024                         }
01025                         Iptres = buffer[0];
01026                         Nspm = buffer[1];
01027                         Nspsol = buffer[2];
01028                         delete [] buffer;
01029 //              skipeoln(file);
01030                         readtoeoln(file);
01031                 }
01032                 Boundary = (int *) get(sizeof(int)*Nspm);
01033                 if (amber7_format)
01034                 {       if (!moveto(file,"ATOMS_PER_MOLECULE"))
01035                         { genclose(file);
01036                           return 0;
01037                         }
01038                         for (i=0; i<Nspm; i++)
01039                                 fscanf(file, " %d", &Boundary[i]);
01040                 }
01041                 else
01042                 {       if (!read_fortran_12I6(file,Boundary,Nspm))
01043                         { genclose(file);
01044                           return 0;
01045                         }
01046 //              skipeoln(file);
01047                         readtoeoln(file);
01048                 }
01049                 if (amber7_format)
01050                         if (!moveto(file,"BOX_DIMENSIONS"))
01051                         { genclose(file);
01052                           return 0;
01053                         }
01054 #ifdef DOUBLE
01055                 fscanf(file, " %lf %lf %lf", 
01056 #else
01057                 fscanf(file, " %f %f %f", 
01058 #endif
01059                                 &Box[0], &Box[1], &Box[2]);
01060 //              skipeoln(file);
01061                 readtoeoln(file);
01062                 if (Iptres)
01063                         Ipatm = Ipres[Iptres] - 1; 
01064                 /* IF(IPTRES.GT.0) IPTATM = IX(I02+IPTRES-1+1)-1 */
01065         }
01066 
01067         /*
01068          * ----- LOAD THE CAP INFORMATION IF NEEDED -----
01069          */
01070 
01071         // I don't know the label for it, so I don't read it if
01072         // it's AMBER 7 format. It's not used in NAMD anyway.
01073 //      if (IfCap) {
01074         if (IfCap && !amber7_format) {
01075                 /* if (IfBox) 
01076                         skipeoln(file); */
01077 #ifdef DOUBLE
01078                 fscanf(file, " %d %lf %lf %lf %lf", 
01079 #else
01080                 fscanf(file, " %d %f %f %f %f", 
01081 #endif
01082                                 &Natcap, &Cutcap, 
01083                                 &Xcap, &Ycap, &Zcap);
01084         }
01085         genclose(file);
01086         if (debug) {
01087                 printf("rdprm done\n");
01088 #ifndef NAMD_NO_STDOUT_FLUSH
01089                 fflush(stdout);
01090 #endif
01091         }
01092         data_read = 1;
01093         return(1);
01094 }
01095 
01096 int Ambertoppar::firstwat()
01097 {
01098         char    *restr = ResNames; 
01099         char    *lastres = ResNames + Nres * 4 + 1;
01100         int     res = 0;
01101 
01102         /*
01103          *  find 1st water residue
01104          */
01105 
01106         for (; restr<lastres; restr+=4) {
01107                 if (!strncmp(restr, "WAT ", 4)) {
01108                   printf("first water: res = %d, atom = %d (%.4s)\n", 
01109                                         res+1, Ipres[res],
01110                                         &AtomNames[Ipres[res]]);
01111 #ifndef NAMD_NO_STDOUT_FLUSH
01112                   fflush(stdout);
01113 #endif
01114                   return(Ipres[res]-1);
01115                 }
01116                 res++;
01117         }
01118         return(0);
01119 }
01120 
01121 
01122 // Constructer: simply set all the pointers Null
01123 
01124 Ambertoppar::parm()
01125 {
01126   data_read = 0;        // No data are read yet
01127   AtomNames = ResNames = AtomSym = AtomTree = NULL;
01128   Charges = Masses = Rk = Req = Tk = Teq = Pk = Pn = Phase = NULL;
01129   Solty = Cn1 = Cn2 = HB12 = HB6 = NULL;
01130   Iac = Iblo = Cno = Ipres = ExclAt = TreeJoin = AtomRes = NULL;
01131   BondHAt1 = BondHAt2 = BondHNum = BondAt1 = BondAt2 = NULL;
01132   BondNum = AngleHAt1 = AngleHAt2 = AngleHAt3 = AngleHNum = NULL;
01133   AngleAt1 = AngleAt2 = AngleAt3 = AngleNum = DihHAt1 = NULL;
01134   DihHAt2 = DihHAt3 = DihHAt4 = DihHNum = DihAt1 = DihAt2 = NULL;
01135   DihAt3 = DihAt4 = DihNum = Boundary = NULL;
01136 }
01137 
01138 
01139 // Destructer: free all the allocated memory for arrays
01140 
01141 Ambertoppar::~parm()
01142 {
01143   free(AtomNames);
01144   free(Charges);
01145   free(Masses);
01146   free(Iac);
01147   free(Iblo);
01148   free(Cno);
01149   free(ResNames);
01150   free(Ipres);
01151   free(Rk);
01152   free(Req);
01153   free(Tk);
01154   free(Teq);
01155   free(Pk);
01156   free(Pn);
01157   free(Phase);
01158   free(Solty);
01159   free(Cn1);
01160   free(Cn2);
01161   free(BondHAt1);
01162   free(BondHAt2);
01163   free(BondHNum);
01164   free(BondAt1);
01165   free(BondAt2);
01166   free(BondNum);
01167   free(AngleHAt1);
01168   free(AngleHAt2);
01169   free(AngleHAt3);
01170   free(AngleHNum);
01171   free(AngleAt1);
01172   free(AngleAt2);
01173   free(AngleAt3);
01174   free(AngleNum);
01175   free(DihHAt1);
01176   free(DihHAt2);
01177   free(DihHAt3);
01178   free(DihHAt4);
01179   free(DihHNum);
01180   free(DihAt1);
01181   free(DihAt2);
01182   free(DihAt3);
01183   free(DihAt4);
01184   free(DihNum);
01185   free(ExclAt);
01186   free(HB12);
01187   free(HB6);
01188   free(AtomSym);
01189   free(AtomTree);
01190   free(TreeJoin);
01191   free(AtomRes);
01192 }
01193 
01194 
01195 // Read FORTRAN 12I6 format data (no space between adjacent data is assumed)
01196 // One needs to read the whole data block into the buffer here
01197 // fp - file pointer.
01198 // data - the buffer to hold the whole block. Should be allocated before the call.
01199 // count - number of ints in the block.
01200 
01201 int Ambertoppar::read_fortran_12I6(FILE *fp, int *data, int count)
01202 {
01203   int i,j;
01204   char buf[7];
01205   
01206   for (i=0; i<count; ++i)
01207   { for (j=0; j<6; ++j)
01208     { buf[j]=getc(fp);
01209       if (buf[j]=='\n' || buf[j]=='\0' || buf[j]==EOF)
01210         return 0;
01211     }
01212     buf[6] = '\0';
01213     if (sscanf(buf,"%d",data+i) != 1)
01214       return 0;
01215     if (i%12==11 && i<count-1)
01216       readtoeoln(fp);
01217   }
01218   
01219   return 1;
01220 }
01221 
01222 
01223 // Used for AMBER 7 format. The function moves the file postion to
01224 // the beginning of the data section labeled by "label"
01225 
01226 int Ambertoppar::moveto(FILE *fp, char *label)
01227 {
01228   char s[76],buf[81];
01229 
01230  while ( 1 ) {
01231   // Find the string "%FLAG"
01232   do preadln(fp, "parm file", buf);
01233   while (strncmp(buf,"%FLAG",5));
01234 
01235   // See if the label is what we expected
01236   sscanf(buf+5,"%s",s);
01237   if (strcasecmp(s,label))
01238     iout << iWARN << "Skipping " << s << " in parm file while seeking " << label << ".\n" << endi;
01239   else
01240     break;
01241  }
01242 
01243   // The next line should begin with "%FORMAT"
01244   preadln(fp, "parm file", buf);
01245   if (strncmp(buf,"%FORMAT",7))
01246     return 0;
01247 
01248   return 1;
01249 }

Generated on Fri Sep 22 01:17:14 2017 for NAMD by  doxygen 1.4.7