src/contrib/scripts/amber2gmxrtp

   1 #!/bin/awk -f
   2
   3 # reads Amber building block definition file(s) and converts to Gromacs
   4 # Version 1.1
   5 # Copyright (c) 2002
   6 # Anton Feenstra - Vrije Universiteit Amsterdam - The Netherlands
   7
   8 function fill_with_underscores(s, n) {
   9   ss=s;
  10   while (length(ss)<n)
  11     ss=ss"_";
  12   return ss;
  13 }
  14
  15 # keep track of which section ('card') of the residue entry we are reading
  16 # if tp<0, 'type' is incremented
  17 function next_type(tp) {
  18   if (tp>=0)
  19     type=tp;
  20   else
  21     type ++;
  22   if (debug) printf("Start reading section: %d = %s\n", type, section[type]);
  23 }
  24
  25 # read list of dihedrals as defined in the ff files (written by amber2gmxff)
  26 function read_dihs(file, dihs) {
  27   ndihs=0;
  28   printf("Reading dihs from '%s'... ", file);
  29   while ( 0<getline <file ) {
  30     n=split($0,d);
  31     for(i=0; i<n; i++)
  32       dihs[ndihs,i]=d[i+1];
  33     ndihs++;
  34   }
  35   close(file);
  36   if (ndihs)
  37     printf("Read %d dihs\n", ndihs);
  38   else
  39     warning("could not read dihs");
  40   return ndihs;
  41 }
  42
  43 # locate dihedrals in list and return number of matches
  44 function find_dih(ai, aj, ak, al, ndihs, dihs, found_dihs) {
  45   nmatch=0;
  46   dih[0]=ai; dih[1]=aj; dih[2]=ak; dih[3]=al;
  47   for(fi=0; fi<ndihs; fi++) {
  48     bMatch=1;
  49     if (debug) printf(".");
  50     for(fd=0; fd<=1; fd++) { # forward and backward comparison
  51       for(fj=0; fj<4 && bMatch; fj++)
  52         if (dihs[fi,fj]!="X")
  53           if (fd==0) # forward comparison
  54             bMatch = (dihs[fi,fj]==dih[fj]);
  55           else       # backward comparison
  56             bMatch = (dihs[fi,fj]==dih[3-fj]);
  57       if (bMatch) {
  58         for(fj=0; fj<4; fj++)
  59           found_dihs[nmatch,fj]=dihs[fi,fj];
  60         nmatch++;
  61       }
  62     }
  63   }
  64   return nmatch;
  65 }
  66
  67 # write list of dihedrals
  68 function dump_dihs(ndihs, dihs) {
  69   for(i=0; i<ndihs; i++) {
  70     printf("dih %d:", i);
  71     for(j=0; j<4; j++)
  72       printf(" %s", dihs[i,j]);
  73     printf("\n");
  74   }
  75 }
  76
  77 function anm2atp(anmi, nat, anms, atps) {
  78   anm=anmi;
  79   gsub("-","",anm);
  80   gsub("\\+","",anm);
  81   for(ii=0; ii<nat; ii++)
  82     if (anm==anms[ii])
  83       return atps[ii];
  84 }
  85
  86 function is_hydrogen(anr) {
  87   return index(isymbl[anr], "H")==1;
  88 }
  89
  90 function swap_atoms(atnr1, atnr2) {
  91   # swap array entries:
  92   dum=igraph[atnr1]; igraph[atnr1]=igraph[atnr2]; igraph[atnr2]=dum;
  93   dum=isymbl[atnr1]; isymbl[atnr1]=isymbl[atnr2]; isymbl[atnr2]=dum;
  94   for(_j=0; _j<3; _j++)
  95     dum=zz[atnr1,_j]; zz[atnr1,_j]=zz[atnr2,_j];  zz[atnr2,_j] =dum;
  96   dum=theta[atnr1];  theta[atnr1] =theta[atnr2];  theta[atnr2] =dum;
  97   dum=chg[atnr1];    chg[atnr1]   =chg[atnr2];    chg[atnr2]   =dum;
  98   dum=allM[atnr1];   allM[atnr1]  =allM[atnr2];   allM[atnr2]  =dum;
  99   dum=allMtp[atnr1]; allMtp[atnr1]=allMtp[atnr2]; allMtp[atnr2]=dum;
 100   # swap atom number references in z-matrix definition:
 101   for(_i=0; _i<nat; _i++)
 102     for(_j=0; _j<nat; _j++)
 103       if      (zz[_i,_j]==atnr1) zz[_i,_j]=atnr2;
 104       else if (zz[_i,_j]==atnr2) zz[_i,_j]=atnr1;
 105 }
 106
 107 function warning(s) {
 108   printf("WARNING: %s\n", s);
 109   nwarn++;
 110 }
 111
 112 function print_warn() {
 113   if (nwarn)
 114     printf("\nThere were %d warnings\n\n", nwarn);
 115 }
 116
 117 function fatal_error(s) {
 118   if (nwarn) printf("(%d warnings pending)\n", nwarn);
 119   printf("\nFATAL ERROR: %s\n\n", s);
 120   xxit=1;
 121   exit -1;
 122 }
 123
 124 BEGIN {
 125   if (ARGC < 2) {
 126     print("Usage:");
 127     print("amber2gmxrtp [ffname=<name>] [debug=1] all##.in");
 128     print("");
 129     print("Reads Amber residue database file (all##.in) and");
 130     print("writes Gromacs residue topology file (ff<name>.rtp),");
 131     print("Gromacs hydrogen database file (ff<name>.hdb),");
 132 #    print("and Gromacs termini database files (ff<name>-[cn].tdb).");
 133     print("NOTE: For the time being, termini are implemented as");
 134     print("additional residues in the ff<name>.rtp as entries with");
 135     print("residue names '<res>-C' or '<res>-N'.");
 136     print("Default for <name> in output files is 'amber'.");
 137     print("");
 138     print("Use 'debug=1' for extremely verbose output.");
 139     print("");
 140     xxit=1;
 141     exit;
 142   }
 143   if ( !ffname ) ffname="amber";
 144   ffrtp = "ff" ffname ".rtp";
 145   ffhdb = "ff" ffname ".hdb";
 146   ffctdb = "ff" ffname "-c.tdb";
 147   ffntdb = "ff" ffname "-n.tdb";
 148   ffpdihs = ffname "-dihedrals.txt";
 149   ffidihs = ffname "-impropers.txt";
 150   nptps=read_dihs(ffpdihs, pdihtps);
 151   nitps=read_dihs(ffidihs, idihtps);
 152   if (debug) dump_dihs(nptps, pdihtps);
 153   if (debug) dump_dihs(nitps, idihtps);
 154   warning("NOTE: .tdb not implemented yet (dummy files will be written).");
 155   printf("Sending output to: %s\n", ffrtp);
 156   printf("; Amber forcefield converted to Gromacs\n") > ffrtp;
 157   printf("; \n") > ffrtp;
 158   type=0;
 159   nres=0; # residues
 160   nM=0;
 161   bM=0;
 162   # declare sections and names (order corresponds to amber file format)
 163   section[tp_RES_TIT=0]="residue title";
 164   section[tp_OUT_FIL=1]="output file";
 165   section[tp_RES_NAM=2]="residue name";
 166   section[tp_GEO_TYP=3]="geometry type";
 167   section[tp_LOO_CUT=4]="loop cutoff";
 168   section[tp_ATOMS  =5]="atoms";
 169   section[tp_EXTRAS =6]="extras";
 170   section[tp_CHARGE =7]="charge";
 171   section[tp_LOOP   =8]="loop";
 172   section[tp_IMP_DIH=9]="improper";
 173
 174   # database types (third number on first line):
 175   db_tp[db_tp_UNI    =1  ]="United-atom";
 176   db_tp[db_tp_UNI_NT =100]="United-atom, NH3+ terminal";
 177   db_tp[db_tp_UNI_CT =101]="United-atom, COO- terminal";
 178   db_tp[db_tp_ALL    =2  ]="All-atom";
 179   db_tp[db_tp_ALL_NT =200]="All-atom, NH3+ terminal";
 180   db_tp[db_tp_ALL_CT =201]="All-atom, COO- terminal";
 181   db_tp[db_tp_OPLS   =3  ]="OPLS United-atom";
 182   db_tp[db_tp_OPLS_NT=300]="OPLS United-atom, NH3+ terminal";
 183   db_tp[db_tp_OPLS_CT=301]="OPLS United-atom, COO- terminal";
 184   db_suf[db_tp_UNI]=db_suf[db_tp_ALL]=db_suf[db_tp_OPLS]="";
 185   db_suf[db_tp_UNI_NT]=db_suf[db_tp_ALL_NT]=db_suf[db_tp_OPLS_NT]="-N";
 186   db_suf[db_tp_UNI_CT]=db_suf[db_tp_ALL_CT]=db_suf[db_tp_OPLS_CT]="-C";
 187   db_pre[db_tp_UNI_NT]=db_pre[db_tp_ALL_NT]=db_pre[db_tp_OPLS_NT]=1;
 188   # marker for z-matrix entries that should be ignored:
 189   zz_NONE=-12345;
 190
 191   # declare hydrogen database hydrogen addition types (Gromacs manual p 95)
 192   # atom orders:
 193   hotp[htp_o_IJK =0,0]=0; hotp[htp_o_IJK,1]=1; hotp[htp_o_IJK,2]=2;
 194   hotp[htp_o_JIK =1,0]=0; hotp[htp_o_JIK,1]=1; hotp[htp_o_JIK,2]=-1;
 195   hotp[htp_o_JKI =2,0]=0; hotp[htp_o_JKI,1]=-1; hotp[htp_o_JKI,2]=1;
 196   hotp[htp_o_JIKL=3,0]=0; hotp[htp_o_JIKL,1]=1;
 197   hotp[htp_o_JIKL,2]=-2; hotp[htp_o_JIKL,3]=-1;
 198   hotp[htp_o_I   =4,0]=0;
 199   hotp[htp_o_ACE =5,0]=0; hotp[htp_o_ACE,1]=-1; hotp[htp_o_ACE,2]=-2;
 200   # the actual hydrogen types:
 201   htp_n[htp_NONE   =-12345]=1;
 202   htp_n[htp_PLANE  =1]=3; htp_o[htp_PLANE  ]=htp_o_JIK;
 203   htp_n[htp_SINGLE =2]=3; htp_o[htp_SINGLE ]=htp_o_IJK;
 204   htp_n[htp_2PLANE =3]=3; htp_o[htp_2PLANE ]=htp_o_IJK;
 205   htp_n[htp_TTETRA =4]=3; htp_o[htp_TTETRA ]=htp_o_IJK;
 206   htp_n[htp_1TETRA =5]=4; htp_o[htp_1TETRA ]=htp_o_JIKL;
 207   htp_n[htp_2TETRA =6]=3; htp_o[htp_2TETRA ]=htp_o_JKI;
 208   htp_n[htp_WATER  =7]=1; htp_o[htp_WATER  ]=htp_o_I;
 209   htp_n[htp_CARBOX =8]=3; htp_o[htp_CARBOX ]=htp_o_IJK;
 210   htp_n[htp_CARBOXH=9]=3; htp_o[htp_CARBOXH]=htp_o_IJK;
 211
 212   # write dummy .tdb files:
 213   print "[ null ]" > ffctdb;
 214   print "[ null ]" > ffntdb;
 215
 216   # write empty xlateat.dat:
 217   print "0" > "xlateat.dat";
 218 }
 219
 220 NR==1 {
 221   printf("Start reading file %s\n", FILENAME);
 222   printf("; from file %s\n", FILENAME) > ffrtp;
 223   idbgen = $1+0;
 224   irest  = $2+0;
 225   itypf  = $3+0;
 226   printf("; \n") > ffrtp;
 227   printf("; reading database: %s\n", db_tp[itypf]) > ffrtp;
 228   printf("; \n") > ffrtp;
 229   if (debug) print "dbgen:", idbgen, irest, itypf, db_suf[itypf];
 230   if (db_pre[itypf]) {
 231     warning("terminus: will ignore all bonds to preceding residue in this file");
 232   }
 233   next;
 234 }
 235
 236 NR==2 {
 237   namdbf = $1;
 238   if (debug) print "dbname:", namdbf;
 239   next;
 240 }
 241
 242 $1=="STOP" { # end of database file
 243   printf("Finished reading file %s\n", FILENAME);
 244   printf("; end of file %s\n", FILENAME) > ffrtp;
 245   NR=0;
 246   next;
 247 }
 248
 249 $1=="DONE" {
 250   if (nat==0) fatal_error("no atoms for residue " resnm);
 251   if (nM) bM=1;
 252   # fix some exceptions:
 253   swap_first=0;
 254   kill_bond=0;
 255   if (resnm=="ACE" || resnm=="HOH") {
 256     # in 'ACE', HH31 comes before CH3, then HH32 and HH33 follow.
 257     # we need to change this to CH3, HH31-3.
 258     # in 'HOH', H1 comes before O. we need O, H1-2
 259     warning("EXCEPTION: swapping atoms "igraph[0]" and "igraph[1]" in residue "resnm);
 260     swap_first=1;
 261     # also, the bond to 'preceding' residue must go:
 262     warning("EXCEPTION: ignoring bond from "bonds[0,1]" to preceding residue in residue "resnm);
 263     kill_bond=1;
 264   }
 265   if (db_pre[itypf]) { # do we have no preceding residue (e.g. N-ter does not)
 266     kill_bond=1;
 267   }
 268   if (swap_first) {
 269     if (debug) print("SWAP FIRST");
 270     swap_atoms(0, 1);
 271   }
 272   if (kill_bond) {
 273     if (debug) print("KILL BOND");
 274     bonds[0,0]=-1;bonds[0,1]=-1;
 275     zz[0,0]=zz[0,1]=zz[0,2]=zz_NONE;
 276     zz[1,1]=zz[1,2]=zz_NONE;
 277     zz[2,2]=zz_NONE;
 278   }
 279   # write rtp
 280   printf("[ %s ]\n", resnmt) > ffrtp;
 281   printf(" [ atoms ]\n") > ffrtp;
 282   printf(";%4s %5s %8s %5s ; %s\n", "name", "type", "charge", "cg", "qtot") > ffrtp;
 283   cg=0;
 284   qt=0;
 285   for(i=0; i<nat; i++) {
 286     qt+=chg[i];
 287     printf("%5s %5s %8.3f %5d ; %g\n", igraph[i],isymbl[i], chg[i],cg++, qt) > ffrtp;
 288   }
 289
 290   # write hdb
 291   nha=0;
 292   for(i=0; i<nat; i++) {
 293     if (is_hydrogen(i)) {
 294       hi[nha]=i-1;
 295       nhh[nha]=0;
 296       while (is_hydrogen(i) && i<nat) {
 297         i++;
 298         nhh[nha]++;
 299       }
 300       nha++;
 301     }
 302   }
 303   printf("%s\t%d\n", resnmt, nha) > ffhdb;
 304   for(i=0; i<nha; i++) {
 305     if (nhh[i]==3) { htp=htp_TTETRA; # three must be tetraeder
 306     } else if (theta[hi[i]+1]==109.5 ||
 307                theta[hi[i]+1]==109.47 ) { # tetraedrical or OH
 308       if (nhh[i]==2) {
 309         if (r[hi[i]+1]==1.01) htp=htp_TTETRA; # -NH2: (-NH3+ after deprotonation)
 310         else htp=htp_2TETRA; # C-CH2-C
 311       } else if (nhh[i]==1) { # either (tert)C-H or OH
 312         if (r[hi[i]+1]==0.96) htp=htp_SINGLE; # OH
 313         else htp=htp_1TETRA; # (tert)C-H
 314       }
 315     } else if (theta[hi[i]+1]==119.8 || # backbone N-H / arg N-H2
 316                theta[hi[i]+1]==118.5 || # arginine N-H
 317                ( theta[hi[i]+1]>=120.0 &&
 318                  theta[hi[i]+1]<=126.0 ) ) { # his/trp N-H's
 319       if (nhh[i]==1) htp=htp_PLANE;
 320       else if (nhh[i]==2) htp=htp_2PLANE;
 321     } else if (theta[hi[i]+1]==109.47 || # OH
 322                theta[hi[i]+1]==118.5  ||
 323                theta[hi[i]+1]==127.0  || # nucl ac.
 324                theta[hi[i]+1]==117.7  || # nucl ac.
 325                theta[hi[i]+1]==117.36 || # nucl ac.
 326                theta[hi[i]+1]==116.77 || # nucl ac.
 327                theta[hi[i]+1]==114.97 || # nucl ac.
 328                theta[hi[i]+1]==107.0  || # nucl ac.
 329                theta[hi[i]+1]==113.0  ||
 330                theta[hi[i]+1]== 96.0 ) { # backbone N-H
 331       htp=htp_SINGLE;
 332     } else if (theta[hi[i]+1]==101.43) {
 333       htp=htp_WATER;
 334     } else if (theta[hi[i]+1]== 60.0) { # PRO-N (-NH2-)
 335       htp=htp_2TETRA;
 336     } else
 337       htp=htp_NONE;
 338     if (debug) printf("i:%d hi:%d(%s) nhh:%d htp:%d htpn:%d theta:%g r:%g\n",
 339                       i, hi[i], igraph[hi[i]], nhh[i], htp, htp_n[htp],
 340                       theta[hi[i]+1], r[hi[i]+1]);
 341     printf("\t%d\t%d", nhh[i], htp) > ffhdb;
 342     # more exceptions:
 343     if (resnm=="ACE" && htp==htp_TTETRA)
 344       htpo=htp_o_ACE;
 345     else
 346       htpo=htp_o[htp];
 347     for(j=0; j<htp_n[htp]; j++) {
 348       k=hotp[htpo,j];
 349       if (debug) printf("j:%d htp:%d htp_o:%d k:%d",
 350                         j, htp, htpo, k);
 351       if (k==0) # atom 'i': bonded atom
 352         l = igraph[hi[i]];
 353       else if (k<0) { # atoms 'ijk': reference atoms upstream
 354         l = "";
 355         # search other non-hydrogen z-matrix entries for 'this' atom number:
 356         for(jj=0; jj<3 && k<0; jj++)
 357           for(ii=hi[i]+1; ii<nat && k<0; ii++) {
 358             if (debug) printf("[ii:%d H:%d jj:%d zz:%d]",
 359                               ii, is_hydrogen(ii), jj, zz[ii,jj]);
 360             if (!is_hydrogen(ii) && zz[ii,jj]!=zz_NONE && zz[ii,jj]==hi[i]) {
 361               if (debug) printf(":%s", igraph[ii]);
 362               k++; # k<0, we count up to zero
 363             }
 364           }
 365         if (k==0)
 366           l = igraph[ii-1];
 367         else {
 368           if (nl) { # finally, check 'loops':
 369             if (debug) printf(" nl:%d", nl);
 370             for(ii=0; ii<nl; ii++)
 371               for(jj=0; jj<2; jj++) {
 372                 if (debug) printf(" %s", loops[ii,jj]);
 373                 if (loops[ii,jj]==igraph[hi[i]])
 374                   l = loops[ii,1-jj];
 375               }
 376             if (debug) printf("\n");
 377           }
 378           if (l=="")
 379             fatal_error("no bonded atom found for atom "igraph[hi[i]]" ("hi[i]+1") in residue "resnmt);
 380         }
 381       } else { # atoms 'ijk': reference atoms downstream (take from z-matrix)
 382         if (zz[hi[i],k-1]<0) # atom in previous residue, use 'allM':
 383           l = "-"allM[nM+zz[hi[i],k-1]];
 384         else # atom in this residue:
 385           l = igraph[zz[hi[i],k-1]];
 386       }
 387       if (debug) printf(" l:%s\n",l);
 388       printf("\t%s", l) > ffhdb;
 389     }
 390     printf("\n") > ffhdb;
 391   }
 392   # done hdb
 393
 394   if (nat>1) { # bonds etc. do not make sense for a single atom:
 395
 396     if (nb==0) fatal_error("no bonds for residue " resnmt);
 397     printf(" [ bonds ]\n") > ffrtp;
 398     printf(";%4s %5s\n", "ai", "aj") > ffrtp;
 399     ignb=0;
 400     for(i=0; i<nb; i++) {
 401       if (bonds[i,0]<0 && bonds[i,1]<0)
 402         ignb++;
 403       else {
 404         if (bonds[i,0]=="") bonds[i,0]="-"allM[nM-1];
 405         printf("%5s %5s\n", bonds[i,0], bonds[i,1]) > ffrtp;
 406       }
 407     }
 408     nb-=ignb;
 409     if (np>0) {
 410       bFirstPdih=1;
 411       for(i=0; i<np; i++) {
 412         # expand 'negative' atom id's:
 413         for(j=0; j<4; j++) {
 414           atp = pdihsatp[i,j]+0;
 415           if (atp<0) {
 416             # fix atom name:
 417             pdihs[i,j] = "-"allM[nM+atp];
 418             # fix atom type:
 419             pdihsatp[i,j] = allMtp[nM+atp];
 420           }
 421         }
 422         # find number of pre-defined 4-atom based dihedral defines, if zero
 423         # we use implicit (2-atom based) dihedrals (which need not be listed):
 424         n = find_dih(pdihsatp[i,0],pdihsatp[i,1],pdihsatp[i,2],pdihsatp[i,3],
 425                      nptps, pdihtps, null);
 426         for(k=0; k<n || (n==0 && k==0); k++)
 427           if ( n ) {
 428             if (bFirstPdih) {
 429               bFirstPdih=0;
 430               printf(" [ dihedrals ]\n") > ffrtp;
 431               printf(";%4s %5s %5s %5s\n", "ai", "aj", "ak", "al") > ffrtp;
 432             }
 433             for(j=0; j<4; j++)
 434               printf("%5s ", pdihs[i,j]) > ffrtp;
 435             printf("    ad") > ffrtp;
 436             for(j=0; j<4; j++)
 437               printf("_%2s", fill_with_underscores(pdihsatp[i,j],2) ) > ffrtp;
 438             if ( n>0 ) printf("_%d", k+1) > ffrtp;
 439             printf("\n") > ffrtp;
 440           }
 441       }
 442     }
 443     if (ni) {
 444       printf(" [ impropers ]\n") > ffrtp;
 445       printf(";%4s %5s %5s %5s\n", "ai", "aj", "ak", "al") > ffrtp;
 446       for(i=0; i<ni; i++) {
 447         for(j=0; j<4; j++) {
 448           if (idihs[i,j]=="-M") idihs[i,j]="-"allM[nM-1];
 449           if (idihs[i,j]=="+M") idihs[i,j]="+"allM[0];
 450           atps[j] = anm2atp(idihs[i,j], nat, igraph, isymbl);
 451           printf("%5s ", idihs[i,j]) > ffrtp;
 452         }
 453         if (debug) {
 454           printf("%s %s %s %s ", atps[0], atps[1], atps[2], atps[3]);
 455           printf("%s %s %s %s ", idihs[i,0], idihs[i,1], idihs[i,2], idihs[i,3]);
 456         }
 457         n = find_dih(atps[0], atps[1], atps[2], atps[3], nitps, idihtps, fdihs);
 458         if (debug) printf("\n");
 459         if (n==1) {
 460           printf("    ai") > ffrtp;
 461           for(j=0; j<4; j++)
 462             printf("_%2s", fill_with_underscores(fdihs[0,j],2) ) > ffrtp;
 463         } else
 464           warning((n?"multiple ("n")":"no")" idih definition for "atps[0]" "atps[1]" "atps[2]" "atps[3]);
 465         printf("\n") > ffrtp;
 466
 467       }
 468     }
 469   }
 470   printf("; end residue %s: %d atoms, %d bonds, %d pdihs, %d idihs\n",
 471          resnmt, nat, nb, np, ni) > ffrtp;
 472   printf("\n") > ffrtp;
 473
 474   # clean up
 475   for(i in xlate) delete xlate[i];
 476
 477   if (debug) print "done.";
 478   next_type(tp_RES_TIT);
 479   next;
 480 }
 481
 482 type==tp_RES_TIT { # title
 483   # initialize per-residue variables
 484   nat=0;
 485   nb=0;
 486   nl=0;
 487   ni=0;
 488   np=0;
 489   # read this section
 490   title = $0;
 491   gsub("   *", " ", title);
 492   if (debug) print "title:", title;
 493   printf("; %s\n", title) > ffrtp;
 494   printf("; \n") > ffrtp;
 495   next_type(-1);
 496   next;
 497 }
 498
 499 type==tp_OUT_FIL { # file
 500   file = $0;
 501   if (debug) print "file:", file;
 502   next_type(-1);
 503   next;
 504 }
 505
 506 type==tp_RES_NAM { # residue name
 507   resnm = $1;
 508   intx  = $2;
 509   kform = $3+0;
 510   if (debug) print "residue:", nres, resnm, intx, kform;
 511   if (length(resnm)>3) {
 512     warning("Residue name '"resnm"' more than three characters long");
 513   }
 514   # add possible suffix for terminus:
 515   resnmt=resnm db_suf[itypf];
 516
 517   printf("Residue: %d %s (%s)\n", nres, title, resnmt);
 518
 519   resnms[nres]=resnm;
 520   nres++;
 521
 522   next_type(-1);
 523   next;
 524 }
 525
 526 type==tp_GEO_TYP { # geometry
 527   ifixc  = $1;
 528   iomit  = $2;
 529   isymdu = $3;
 530   ipos   = $4;
 531   if (debug) print "geometry:", ifixc, iomit, isymdu, ipos;
 532   next_type(-1);
 533   next;
 534 }
 535
 536 type==tp_LOO_CUT { # cut
 537   cut=$1+0;
 538   if (debug) print "cut:", cut;
 539   next_type(-1);
 540   next;
 541 }
 542
 543 type==tp_ATOMS && NF==0 { # end of atoms
 544   ndum=0;
 545   next_type(-1);
 546   next;
 547 }
 548
 549 type==tp_ATOMS { # atoms
 550   if ($2=="DUMM") {
 551     ndum++;
 552   } else {
 553     i         =  $1-ndum-1;
 554     igraph[i] =  $2;
 555     isymbl[i] =  $3;
 556     itree     =  $4;
 557     zz[i,0]   =  $5-ndum-1;
 558     zz[i,1]   =  $6-ndum-1;
 559     zz[i,2]   =  $7-ndum-1;
 560     r[i]      =  $8+0;
 561     theta[i]  =  $9+0;
 562     phi       = $10+0;
 563     chg[i]    = $11+0;
 564     nat++;
 565     # fix hydrogen atom numbering:
 566     if (is_hydrogen(i) && length(igraph[i])>length_heavy) {
 567       hnum++;
 568       igr=substr(igraph[i], 1, length_heavy) hnum;
 569       if (igr!=igraph[i]) {
 570         # store translated atom names:
 571         if (debug) print "Xlate:", i, "'"igraph[i]"'", "'"isymbl[i]"'",
 572                      length_heavy, "'"igr"'";
 573         xlate[igraph[i]]=igr;
 574         igraph[i]=igr;
 575       }
 576     } else {
 577       length_heavy=length(igraph[i])
 578       hnum=0;
 579     }
 580     if (itree=="M" && !bM) {
 581       allM[nM]   = igraph[i];
 582       allMtp[nM] = isymbl[i];
 583       nM++;
 584     }
 585     # add bonds
 586     bonds[nb, 0]=igraph[zz[i,0]];
 587     bonds[nb, 1]=igraph[i];
 588     nb++;
 589     # angles are automatic
 590     # add proper dihs
 591     pdihs[np, 0]=igraph[i];
 592     pdihs[np, 1]=igraph[zz[i,0]];
 593     pdihs[np, 2]=igraph[zz[i,1]];
 594     pdihs[np, 3]=igraph[zz[i,2]];
 595     pdihsatp[np, 0]=isymbl[i];
 596     pdihsatp[np, 1]= zz[i,0]<0 ? zz[i,0] : isymbl[zz[i,0]];
 597     pdihsatp[np, 2]= zz[i,1]<0 ? zz[i,1] : isymbl[zz[i,1]];
 598     pdihsatp[np, 3]= zz[i,2]<0 ? zz[i,2] : isymbl[zz[i,2]];
 599     np++;
 600     if (debug) printf("atom: %d %s %s %s %d %d %d %g %g %g %g\n",
 601                       i+1, igraph[i], isymbl[i], itree,
 602                       zz[i,0], zz[i,1], zz[i,2], r[i], theta[i], phi, chg[i]);
 603   }
 604 }
 605
 606 type==tp_EXTRAS { # extra's
 607   if (debug) print "extra:", $0;
 608   for(i=tp_CHARGE; i<=tp_IMP_DIH; i++) {
 609     if (toupper(section[i]) == $1) {
 610       next_type(i);
 611       next;
 612     }
 613   }
 614 }
 615
 616 type>tp_EXTRAS && NF==0 { # end of an extra's section
 617   if (iq && iq!=nat)
 618     warning("ERROR: # charges ("iq") and # atoms ("nat") differ");
 619   iq=0;
 620   next_type(tp_EXTRAS);
 621   next;
 622 }
 623
 624 type==tp_CHARGE { # charge
 625   for(i=1; i<=NF; i++) {
 626     chg[iq] = $i;
 627     iq++;
 628   }
 629   if (debug) print "charges:", iq;
 630 }
 631
 632 type==tp_LOOP { # loop
 633   for(i=1; i<=2; i++) bonds[nb, i-1] = xlate[$i] ? xlate[$i] : $i;
 634   nb++;
 635   for(i=1; i<=2; i++) loops[nl, i-1] = xlate[$i] ? xlate[$i] : $i;
 636   nl++;
 637 }
 638
 639 type==tp_IMP_DIH { # improper
 640   for(i=1; i<=4; i++) idihs[ni, i-1] = xlate[$i] ? xlate[$i] : $i;
 641   ni++;
 642 }
 643
 644 END {
 645   if (xxit) exit;
 646   printf("Found %d residues\n", nres);
 647   print_warn();
 648   printf("; Found %d residues\n", nres) > ffrtp;
 649 }
 650
 651 #last line