NAMD
parm.C
Go to the documentation of this file.
1 
2 /*
3  * COPYRIGHT 1992, REGENTS OF THE UNIVERSITY OF CALIFORNIA
4  *
5  * prm.c - read information from an amber PARM topology file:
6  * atom/residue/bond/charge info, plus force field data.
7  * This file and the accompanying prm.h may be distributed
8  * provided this notice is retained unmodified and provided
9  * that any modifications to the rest of the file are noted
10  * in comments.
11  *
12  * Bill Ross, UCSF 1994
13  */
14 
15 // The original implementation requires space between adjacent
16 // data. This is not true for AMBER 6 format files if it's a
17 // large system. So a large part of the code is modified. All
18 // reads of integer in 12I6 format is now implemented using the
19 // new method, which doesn't assume space between entries.
20 // Also, code has been added to read AMBER 7 format files.
21 // Apr. 9, 2002
22 
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <ctype.h>
26 #include <sys/types.h>
27 #include <sys/stat.h>
28 #include <errno.h>
29 #include <string.h>
30 #include "strlib.h" // For strcasecmp and strncasecmp
31 
32 #include "common.h"
33 #include "InfoStream.h"
34 #include "parm.h"
35 
36 static int debug = 0; /* set it if you want */
37 
38 /* fortran formats
39  * 9118 FORMAT(12I6)
40  * 9128 FORMAT(5E16.8)
41  */
42 
43 // This trick is misleading and still needs to have space between
44 // data thus won't work for 6 digit integers. Don't use it.
45 //const char *f9118 = "%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d%6d\n";
46 
47 
48 // This function was borrowed from VMD code in "ReadPARM.C". Here it
49 // is used to replace the "skipeoln()" in the original code, which
50 // seems to be undefined.
51 static int readtoeoln(FILE *f) {
52  char c;
53 
54  /* skip to eoln */
55  while((c = getc(f)) != '\n') {
56  if (c == EOF)
57  return -1;
58  }
59 
60  return 0;
61 }
62 
63 /***********************************************************************
64  GENOPEN()
65 ************************************************************************/
66 
67 /*
68  * genopen() - fopen regular or popen compressed file for reading
69  */
70 // replaced with NAMD implementation
71 
72 FILE *Ambertoppar::genopen(const char *name)
73 {
74  return(Fopen(name,"r"));
75 }
76 
77 /***********************************************************************
78  GENCLOSE()
79 ************************************************************************/
80 
81 /*
82  * genclose() - close fopened or popened file
83  */
84 // replaced with NAMD implementation
85 
86 void Ambertoppar::genclose(FILE *fileptr)
87 {
88  Fclose(fileptr);
89 }
90 
91 
92 /***********************************************************************
93  GET()
94 ************************************************************************/
95 
96 char *Ambertoppar::get(int size)
97 {
98  char *ptr;
99 
100 #ifdef DEBUG
101  printf("malloc %d\n", size);
102 #ifndef NAMD_NO_STDOUT_FLUSH
103  fflush(stdout);
104 #endif
105 #endif
106  if (size ==0)
107  return((char *) NULL);
108 
109  if ((ptr = (char *) malloc((unsigned)size)) == NULL) {
110  printf("malloc %d", size);
111 #ifndef NAMD_NO_STDOUT_FLUSH
112  fflush(stdout);
113 #endif
114  NAMD_die("Memory allocation error in Ambertoppar::get()");
115  }
116  return(ptr);
117 }
118 
119 /***********************************************************************
120  PREADLN()
121 ************************************************************************/
122 
123 void Ambertoppar::preadln(FILE *file, const char *name, char *string)
124 {
125  int i, j;
126 
127  for (i=0; i<81; i++) {
128  if ((j = getc(file)) == EOF) {
129  NAMD_die("Unexpected EOF in Amber parm file");
130  printf("Error: unexpected EOF in %s\n", name);
131  }
132  string[i] = (char) j;
133  if (string[i] == '\n') {
134  break;
135  }
136  }
137  if (i == 80 && string[i] != '\n') {
138  NAMD_die("Line too long in Amber parm file");
139  printf("Error: line too long in %s:\n%.80s", name, string);
140  }
141 }
142 
143 /***************************************************************************
144  READPARM()
145 ****************************************************************************/
146 
147 /*
148  * readparm() - instantiate a given Ambertoppar
149  */
150 
151 int Ambertoppar::readparm(char *name)
152 {
153  _REAL *H;
154  int i, idum, res, ifpert;
155  int *buffer, amber7_format;
156  FILE *file;
157 
158  if (data_read)
159  { iout << "Duplicate parm data in one object!\n" << endi;
160  return(0);
161  }
162 
163 // printf("Reading parm file (%s)\n", name);
164  iout << "Reading parm file (" << name << ") ...\n" << endi;
165 
166 // if ((file = genopen(name, "parm")) == NULL)
167  if ((file = genopen(name)) == NULL)
168  return(0);
169 
170  /* READ TITLE */
171 
172  preadln(file, name, ititl);
173 // "ititle" doesn't guarantee to have '\0' (as the end of a string),
174 // so the following is disabled in order to avoid strange output
175 // printf("%s title:\n%s", name, ititl);
176 
177  // Check whether it's in AMBER 7 format or in old format
178 
179  if (strncmp(ititl,"%VERSION",8))
180  amber7_format = 0; // old format
181  else
182  { amber7_format = 1; // AMBER 7 format
183  iout << "PARM file in AMBER 7 format\n" << endi;
184  if (!moveto(file,"TITLE"))
185  { genclose(file);
186  return 0;
187  }
188  preadln(file, name, ititl);
189  }
190 
191  /* READ CONTROL INTEGERS */
192 
193  if (amber7_format)
194  { if (!moveto(file,"POINTERS"))
195  { genclose(file);
196  return 0;
197  }
198 // fscanf(file, f9118,
199  fscanf(file, "%d%d%d%d%d%d%d%d%d%d%d%d",
200  &Natom, &Ntypes, &Nbonh, &Mbona,
201  &Ntheth, &Mtheta, &Nphih, &Mphia,
202  &Nhparm, &Nparm, &Nnb, &Nres);
203 
204 // fscanf(file, f9118,
205  fscanf(file, "%d%d%d%d%d%d%d%d%d%d%d%d",
206  &Nbona, &Ntheta, &Nphia, &Numbnd,
207  &Numang, &Nptra, &Natyp, &Nphb,
208  &ifpert, &idum, &idum, &idum);
209 
210  fscanf(file, " %d %d %d %d %d %d",
211  &idum, &idum,&idum,&IfBox,&Nmxrs,&IfCap);
212  }
213  else
214  { buffer = new int[30];
215  if (!read_fortran_12I6(file,buffer,30))
216  { genclose(file);
217  return 0;
218  }
219  Natom = buffer[0];
220  Ntypes = buffer[1];
221  Nbonh = buffer[2];
222  Mbona = buffer[3];
223  Ntheth = buffer[4];
224  Mtheta = buffer[5];
225  Nphih = buffer[6];
226  Mphia = buffer[7];
227  Nhparm = buffer[8];
228  Nparm = buffer[9];
229  Nnb = buffer[10];
230  Nres = buffer[11];
231  Nbona = buffer[12];
232  Ntheta = buffer[13];
233  Nphia = buffer[14];
234  Numbnd = buffer[15];
235  Numang = buffer[16];
236  Nptra = buffer[17];
237  Natyp = buffer[18];
238  Nphb = buffer[19];
239  ifpert = buffer[20];
240  IfBox = buffer[27];
241  Nmxrs = buffer[28];
242  IfCap = buffer[29];
243  delete [] buffer;
244 // skipeoln(file);
245  readtoeoln(file);
246  }
247 
248  if (ifpert) {
249  printf("not equipped to read perturbation prmtop\n");
250  return(0);
251  }
252 
253  /* ALLOCATE MEMORY */
254 
255  Nat3 = 3 * Natom;
256  Ntype2d = Ntypes * Ntypes;
257  Nttyp = Ntypes*(Ntypes+1)/2;
258 
259  /*
260  * get most of the indirect stuff; some extra allowed for char arrays
261  */
262 
263  AtomNames = (char *) get(4*Natom+81);
264  Charges = (_REAL *) get(sizeof(_REAL)*Natom);
265  Masses = (_REAL *) get(sizeof(_REAL)*Natom);
266  Iac = (int *) get(sizeof(int)*Natom);
267  Iblo = (int *) get(sizeof(int)*Natom);
268  Cno = (int *) get(sizeof(int)* Ntype2d);
269  ResNames = (char *) get(4* Nres+81);
270  Ipres = (int *) get(sizeof(int)*( Nres+1));
271  Rk = (_REAL *) get(sizeof(_REAL)* Numbnd);
272  Req = (_REAL *) get(sizeof(_REAL)* Numbnd);
273  Tk = (_REAL *) get(sizeof(_REAL)* Numang);
274  Teq = (_REAL *) get(sizeof(_REAL)* Numang);
275  Pk = (_REAL *) get(sizeof(_REAL)* Nptra);
276  Pn = (_REAL *) get(sizeof(_REAL)* Nptra);
277  Phase = (_REAL *) get(sizeof(_REAL)* Nptra);
278  Solty = (_REAL *) get(sizeof(_REAL)* Natyp);
279  Cn1 = (_REAL *) get(sizeof(_REAL)* Nttyp);
280  Cn2 = (_REAL *) get(sizeof(_REAL)* Nttyp);
281  BondHAt1 = (int *) get(sizeof(int)* Nbonh);
282  BondHAt2 = (int *) get(sizeof(int)* Nbonh);
283  BondHNum = (int *) get(sizeof(int)* Nbonh);
284  BondAt1 = (int *) get(sizeof(int)* Nbona);
285  BondAt2 = (int *) get(sizeof(int)* Nbona);
286  BondNum = (int *) get(sizeof(int)* Nbona);
287  AngleHAt1 = (int *) get(sizeof(int)* Ntheth);
288  AngleHAt2 = (int *) get(sizeof(int)* Ntheth);
289  AngleHAt3 = (int *) get(sizeof(int)* Ntheth);
290  AngleHNum = (int *) get(sizeof(int)* Ntheth);
291  AngleAt1 = (int *) get(sizeof(int)* Ntheta);
292  AngleAt2 = (int *) get(sizeof(int)*Ntheta);
293  AngleAt3 = (int *) get(sizeof(int)*Ntheta);
294  AngleNum = (int *) get(sizeof(int)*Ntheta);
295  DihHAt1 = (int *) get(sizeof(int)*Nphih);
296  DihHAt2 = (int *) get(sizeof(int)*Nphih);
297  DihHAt3 = (int *) get(sizeof(int)*Nphih);
298  DihHAt4 = (int *) get(sizeof(int)*Nphih);
299  DihHNum = (int *) get(sizeof(int)*Nphih);
300  DihAt1 = (int *) get(sizeof(int)*Nphia);
301  DihAt2 = (int *) get(sizeof(int)*Nphia);
302  DihAt3 = (int *) get(sizeof(int)*Nphia);
303  DihAt4 = (int *) get(sizeof(int)*Nphia);
304  DihNum = (int *) get(sizeof(int)*Nphia);
305  ExclAt = (int *) get(sizeof(int)*Nnb);
306  HB12 = (_REAL *) get(sizeof(_REAL)*Nphb);
307  HB6 = (_REAL *) get(sizeof(_REAL)*Nphb);
308  AtomSym = (char *) get(4*Natom+81);
309  AtomTree = (char *) get(4*Natom+81);
310  TreeJoin = (int *) get(sizeof(int)*Natom);
311  AtomRes = (int *) get(sizeof(int)*Natom);
312 
313  /*
314  * READ ATOM NAMES -IH(M04)
315  */
316 
317  if (amber7_format)
318  if (!moveto(file,"ATOM_NAME"))
319  { genclose(file);
320  return 0;
321  }
322  for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
323  preadln(file, "", &AtomNames[i*80]);
324 
325  /*
326  * READ ATOM CHARGES -X(L15)
327  * (pre-multiplied by an energy factor of 18.2223 == sqrt(332)
328  * for faster force field calculations)
329  */
330 
331  if (amber7_format)
332  if (!moveto(file,"CHARGE"))
333  { genclose(file);
334  return 0;
335  }
336  for (i=0; i<Natom; i++)
337 #ifdef DOUBLE
338  fscanf(file, " %lf", &Charges[i]);
339 #else
340  fscanf(file, " %f", &Charges[i]);
341 #endif
342 // skipeoln(file);
343  readtoeoln(file);
344 
345  /*
346  * READ ATOM MASSES -X(L20)
347  */
348 
349  if (amber7_format)
350  if (!moveto(file,"MASS"))
351  { genclose(file);
352  return 0;
353  }
354  for (i=0; i<Natom; i++)
355 #ifdef DOUBLE
356  fscanf(file, " %le", &Masses[i]);
357 #else
358  fscanf(file, " %e", &Masses[i]);
359 #endif
360 // skipeoln(file);
361  readtoeoln(file);
362 
363  /*
364  * READ ATOM L-J TYPES -IX(I04)
365  */
366 
367  if (amber7_format)
368  { if (!moveto(file,"ATOM_TYPE_INDEX"))
369  { genclose(file);
370  return 0;
371  }
372  for (i=0; i<Natom; i++)
373  fscanf(file, " %d", &Iac[i]);
374  }
375  else
376  { if (!read_fortran_12I6(file,Iac,Natom))
377  { genclose(file);
378  return 0;
379  }
380 // skipeoln(file);
381  readtoeoln(file);
382  }
383 
384  /*
385  * READ ATOM INDEX TO 1st IN EXCLUDED ATOM LIST "NATEX" -IX(I08)
386  */
387 
388  if (amber7_format)
389  { if (!moveto(file,"NUMBER_EXCLUDED_ATOMS"))
390  { genclose(file);
391  return 0;
392  }
393  for (i=0; i<Natom; i++)
394  fscanf(file, " %d", &Iblo[i]);
395  }
396  else
397  { if (!read_fortran_12I6(file,Iblo,Natom))
398  { genclose(file);
399  return 0;
400  }
401 // skipeoln(file);
402  readtoeoln(file);
403  }
404 
405  /*
406  * READ TYPE INDEX TO N-B TYPE -IX(I06)
407  */
408 
409  if (amber7_format)
410  { if (!moveto(file,"NONBONDED_PARM_INDEX"))
411  { genclose(file);
412  return 0;
413  }
414  for (i=0; i<Ntype2d; i++)
415  fscanf(file, " %d", &Cno[i]);
416  }
417  else
418  { if (!read_fortran_12I6(file,Cno,Ntype2d))
419  { genclose(file);
420  return 0;
421  }
422 // skipeoln(file);
423  readtoeoln(file);
424  }
425 
426  /*
427  * READ RES NAMES (4 chars each, 4th blank) -IH(M02)
428  */
429 
430  if (amber7_format)
431  if (!moveto(file,"RESIDUE_LABEL"))
432  { genclose(file);
433  return 0;
434  }
435  for (i=0; i<(Nres/20 + (Nres%20 ? 1 : 0)); i++)
436  preadln(file, "", &ResNames[i*80]);
437 
438  /*
439  * READ RES POINTERS TO 1st ATOM -IX(I02)
440  */
441 
442  if (amber7_format)
443  { if (!moveto(file,"RESIDUE_POINTER"))
444  { genclose(file);
445  return 0;
446  }
447  for (i=0; i<Nres; i++)
448  fscanf(file, " %d", &Ipres[i]);
449  }
450  else
451  { if (!read_fortran_12I6(file,Ipres,Nres))
452  { genclose(file);
453  return 0;
454  }
455 // skipeoln(file);
456  readtoeoln(file);
457  }
458  Ipres[Nres] = Natom + 1;
459 
460  /*
461  * READ BOND FORCE CONSTANTS -RK()
462  */
463 
464  if (amber7_format)
465  if (!moveto(file,"BOND_FORCE_CONSTANT"))
466  { genclose(file);
467  return 0;
468  }
469  for (i=0; i< Numbnd; i++)
470 #ifdef DOUBLE
471  fscanf(file, " %lf", &Rk[i]);
472 #else
473  fscanf(file, " %f", &Rk[i]);
474 #endif
475 // skipeoln(file);
476  readtoeoln(file);
477 
478  /*
479  * READ BOND LENGTH OF MINIMUM ENERGY -REQ()
480  */
481 
482  if (amber7_format)
483  if (!moveto(file,"BOND_EQUIL_VALUE"))
484  { genclose(file);
485  return 0;
486  }
487  for (i=0; i< Numbnd; i++)
488 #ifdef DOUBLE
489  fscanf(file, " %lf", &Req[i]);
490 #else
491  fscanf(file, " %f", &Req[i]);
492 #endif
493 // skipeoln(file);
494  readtoeoln(file);
495 
496  /*
497  * READ BOND ANGLE FORCE CONSTANTS (following Rk nomen) -TK()
498  */
499 
500  if (amber7_format)
501  if (!moveto(file,"ANGLE_FORCE_CONSTANT"))
502  { genclose(file);
503  return 0;
504  }
505  for (i=0; i< Numang; i++)
506 #ifdef DOUBLE
507  fscanf(file, " %lf", &Tk[i]);
508 #else
509  fscanf(file, " %f", &Tk[i]);
510 #endif
511 // skipeoln(file);
512  readtoeoln(file);
513 
514  /*
515  * READ BOND ANGLE OF MINIMUM ENERGY (following Req nomen) -TEQ()
516  */
517 
518  if (amber7_format)
519  if (!moveto(file,"ANGLE_EQUIL_VALUE"))
520  { genclose(file);
521  return 0;
522  }
523  for (i=0; i< Numang; i++)
524 #ifdef DOUBLE
525  fscanf(file, " %lf", &Teq[i]);
526 #else
527  fscanf(file, " %f", &Teq[i]);
528 #endif
529 // skipeoln(file);
530  readtoeoln(file);
531 
532  /*
533  * READ DIHEDRAL PEAK MAGNITUDE -PK()
534  */
535 
536  if (amber7_format)
537  if (!moveto(file,"DIHEDRAL_FORCE_CONSTANT"))
538  { genclose(file);
539  return 0;
540  }
541  for (i=0; i< Nptra; i++)
542 #ifdef DOUBLE
543  fscanf(file, " %lf", &Pk[i]);
544 #else
545  fscanf(file, " %f", &Pk[i]);
546 #endif
547 // skipeoln(file);
548  readtoeoln(file);
549 
550  /*
551  * READ DIHEDRAL PERIODICITY -PN()
552  */
553 
554  if (amber7_format)
555  if (!moveto(file,"DIHEDRAL_PERIODICITY"))
556  { genclose(file);
557  return 0;
558  }
559  for (i=0; i< Nptra; i++)
560 #ifdef DOUBLE
561  fscanf(file, " %lf", &Pn[i]);
562 #else
563  fscanf(file, " %f", &Pn[i]);
564 #endif
565 // skipeoln(file);
566  readtoeoln(file);
567 
568  /*
569  * READ DIHEDRAL PHASE -PHASE()
570  */
571 
572  if (amber7_format)
573  if (!moveto(file,"DIHEDRAL_PHASE"))
574  { genclose(file);
575  return 0;
576  }
577  for (i=0; i< Nptra; i++)
578 #ifdef DOUBLE
579  fscanf(file, " %lf", &Phase[i]);
580 #else
581  fscanf(file, " %f", &Phase[i]);
582 #endif
583 // skipeoln(file);
584  readtoeoln(file);
585 
586  /*
587  * ?? "RESERVED" -SOLTY()
588  */
589 
590  if (amber7_format)
591  if (!moveto(file,"SOLTY"))
592  { genclose(file);
593  return 0;
594  }
595  for (i=0; i< Natyp; i++)
596 #ifdef DOUBLE
597  fscanf(file, " %lf", &Solty[i]);
598 #else
599  fscanf(file, " %f", &Solty[i]);
600 #endif
601 // skipeoln(file);
602  readtoeoln(file);
603 
604  /*
605  * READ L-J R**12 FOR ALL PAIRS OF ATOM TYPES -CN1()
606  * (SHOULD BE 0 WHERE H-BONDS)
607  */
608 
609  if (amber7_format)
610  if (!moveto(file,"LENNARD_JONES_ACOEF"))
611  { genclose(file);
612  return 0;
613  }
614  for (i=0; i< Nttyp; i++)
615 #ifdef DOUBLE
616  fscanf(file, " %lf", &Cn1[i]);
617 #else
618  fscanf(file, " %f", &Cn1[i]);
619 #endif
620 // skipeoln(file);
621  readtoeoln(file);
622 
623  /*
624  * READ L-J R**6 FOR ALL PAIRS OF ATOM TYPES -CN2()
625  * (SHOULD BE 0 WHERE H-BONDS)
626  */
627 
628  if (amber7_format)
629  if (!moveto(file,"LENNARD_JONES_BCOEF"))
630  { genclose(file);
631  return 0;
632  }
633  for (i=0; i< Nttyp; i++)
634 #ifdef DOUBLE
635  fscanf(file, " %lf", &Cn2[i]);
636 #else
637  fscanf(file, " %f", &Cn2[i]);
638 #endif
639 // skipeoln(file);
640  readtoeoln(file);
641 
642  /*
643  * READ COVALENT BOND W/ HYDROGEN (3*(atnum-1)):
644  * IBH = ATOM1 -IX(I12)
645  * JBH = ATOM2 -IX(I14)
646  * ICBH = BOND ARRAY PTR -IX(I16)
647  */
648 
649  if (amber7_format)
650  { if (!moveto(file,"BONDS_INC_HYDROGEN"))
651  { genclose(file);
652  return 0;
653  }
654  for (i=0; i<Nbonh; i++)
655  fscanf(file, " %d %d %d",
656  &BondHAt1[i], &BondHAt2[i], &BondHNum[i]);
657  }
658  else
659  { buffer = new int[3*Nbonh];
660  if (!read_fortran_12I6(file,buffer,3*Nbonh))
661  { genclose(file);
662  return 0;
663  }
664  for (i=0; i<Nbonh; i++)
665  { BondHAt1[i] = buffer[3*i];
666  BondHAt2[i] = buffer[3*i+1];
667  BondHNum[i] = buffer[3*i+2];
668  }
669  delete [] buffer;
670 // skipeoln(file);
671  readtoeoln(file);
672  }
673 
674  /*
675  * READ COVALENT BOND W/OUT HYDROGEN (3*(atnum-1)):
676  * IB = ATOM1 -IX(I18)
677  * JB = ATOM2 -IX(I20)
678  * ICB = BOND ARRAY PTR -IX(I22)
679  */
680 
681  if (amber7_format)
682  { if (!moveto(file,"BONDS_WITHOUT_HYDROGEN"))
683  { genclose(file);
684  return 0;
685  }
686  for (i=0; i<Nbona; i++)
687  fscanf(file, " %d %d %d",
688  &BondAt1[i], &BondAt2[i], &BondNum[i]);
689  }
690  else
691  { buffer = new int[3*Nbona];
692  if (!read_fortran_12I6(file,buffer,3*Nbona))
693  { genclose(file);
694  return 0;
695  }
696  for (i=0; i<Nbona; i++)
697  { BondAt1[i] = buffer[3*i];
698  BondAt2[i] = buffer[3*i+1];
699  BondNum[i] = buffer[3*i+2];
700  }
701  delete [] buffer;
702 // skipeoln(file);
703  readtoeoln(file);
704  }
705 
706  /*
707  * READ ANGLE W/ HYDROGEN:
708  * ITH = ATOM1 -IX(I24)
709  * JTH = ATOM2 -IX(I26)
710  * KTH = ATOM3 -IX(I28)
711  * ICTH = ANGLE ARRAY PTR -IX(I30)
712  */
713 
714  if (amber7_format)
715  { if (!moveto(file,"ANGLES_INC_HYDROGEN"))
716  { genclose(file);
717  return 0;
718  }
719  for (i=0; i<Ntheth; i++)
720  fscanf(file, " %d %d %d %d",
721  &AngleHAt1[i], &AngleHAt2[i],
722  &AngleHAt3[i], &AngleHNum[i]);
723  }
724  else
725  { buffer = new int[4*Ntheth];
726  if (!read_fortran_12I6(file,buffer,4*Ntheth))
727  { genclose(file);
728  return 0;
729  }
730  for (i=0; i<Ntheth; i++)
731  { AngleHAt1[i] = buffer[4*i];
732  AngleHAt2[i] = buffer[4*i+1];
733  AngleHAt3[i] = buffer[4*i+2];
734  AngleHNum[i] = buffer[4*i+3];
735  }
736  delete [] buffer;
737 // skipeoln(file);
738  readtoeoln(file);
739  }
740 
741  /*
742  * READ ANGLE W/OUT HYDROGEN:
743  * IT = ATOM1 -IX(I32)
744  * JT = ATOM2 -IX(I34)
745  * KT = ATOM3 -IX(I36)
746  * ICT = ANGLE ARRAY PTR -IX(I38)
747  */
748 
749  if (amber7_format)
750  { if (!moveto(file,"ANGLES_WITHOUT_HYDROGEN"))
751  { genclose(file);
752  return 0;
753  }
754  for (i=0; i<Ntheta; i++)
755  fscanf(file, " %d %d %d %d",
756  &AngleAt1[i], &AngleAt2[i],
757  &AngleAt3[i], &AngleNum[i]);
758  }
759  else
760  { buffer = new int[4*Ntheta];
761  if (!read_fortran_12I6(file,buffer,4*Ntheta))
762  { genclose(file);
763  return 0;
764  }
765  for (i=0; i<Ntheta; i++)
766  { AngleAt1[i] = buffer[4*i];
767  AngleAt2[i] = buffer[4*i+1];
768  AngleAt3[i] = buffer[4*i+2];
769  AngleNum[i] = buffer[4*i+3];
770  }
771  delete [] buffer;
772 // skipeoln(file);
773  readtoeoln(file);
774  }
775 
776  /*
777  * READ DIHEDRAL W/ HYDROGEN:
778  * ITH = ATOM1 -IX(40)
779  * JTH = ATOM2 -IX(42)
780  * KTH = ATOM3 -IX(44)
781  * LTH = ATOM4 -IX(46)
782  * ICTH = DIHEDRAL ARRAY PTR -IX(48)
783  */
784 
785  if (amber7_format)
786  { if (!moveto(file,"DIHEDRALS_INC_HYDROGEN"))
787  { genclose(file);
788  return 0;
789  }
790  for (i=0; i<Nphih; i++)
791  fscanf(file, " %d %d %d %d %d",
792  &DihHAt1[i], &DihHAt2[i], &DihHAt3[i],
793  &DihHAt4[i], &DihHNum[i]);
794  }
795  else
796  { buffer = new int[5*Nphih];
797  if (!read_fortran_12I6(file,buffer,5*Nphih))
798  { genclose(file);
799  return 0;
800  }
801  for (i=0; i<Nphih; i++)
802  { DihHAt1[i] = buffer[5*i];
803  DihHAt2[i] = buffer[5*i+1];
804  DihHAt3[i] = buffer[5*i+2];
805  DihHAt4[i] = buffer[5*i+3];
806  DihHNum[i] = buffer[5*i+4];
807  }
808  delete [] buffer;
809 // skipeoln(file);
810  readtoeoln(file);
811  }
812 
813  /*
814  * READ DIHEDRAL W/OUT HYDROGEN:
815  * IT = ATOM1
816  * JT = ATOM2
817  * KT = ATOM3
818  * LT = ATOM4
819  * ICT = DIHEDRAL ARRAY PTR
820  */
821 
822  if (amber7_format)
823  { if (!moveto(file,"DIHEDRALS_WITHOUT_HYDROGEN"))
824  { genclose(file);
825  return 0;
826  }
827  for (i=0; i<Nphia; i++)
828  fscanf(file, " %d %d %d %d %d",
829  &DihAt1[i], &DihAt2[i], &DihAt3[i],
830  &DihAt4[i], &DihNum[i]);
831  }
832  else
833  { buffer = new int[5*Nphia];
834  if (!read_fortran_12I6(file,buffer,5*Nphia))
835  { genclose(file);
836  return 0;
837  }
838  for (i=0; i<Nphia; i++) {
839  DihAt1[i] = buffer[5*i];
840  DihAt2[i] = buffer[5*i+1];
841  DihAt3[i] = buffer[5*i+2];
842  DihAt4[i] = buffer[5*i+3];
843  DihNum[i] = buffer[5*i+4];
844  }
845  delete [] buffer;
846 // skipeoln(file);
847  readtoeoln(file);
848  }
849 
850  /*
851  * READ EXCLUDED ATOM LIST -IX(I10)
852  */
853 
854  if (amber7_format)
855  { if (!moveto(file,"EXCLUDED_ATOMS_LIST"))
856  { genclose(file);
857  return 0;
858  }
859  for (i=0; i<Nnb; i++)
860  fscanf(file, " %d", &ExclAt[i]);
861  }
862  else
863  { if (!read_fortran_12I6(file,ExclAt,Nnb))
864  { genclose(file);
865  return 0;
866  }
867 // skipeoln(file);
868  readtoeoln(file);
869  }
870 
871  /*
872  * READ H-BOND R**12 TERM FOR ALL N-B TYPES -ASOL()
873  */
874 
875  if (amber7_format)
876  if (!moveto(file,"HBOND_ACOEF"))
877  { genclose(file);
878  return 0;
879  }
880  for (i=0; i<Nphb; i++)
881 #ifdef DOUBLE
882  fscanf(file, " %lf", &HB12[i]);
883 #else
884  fscanf(file, " %f", &HB12[i]);
885 #endif
886 // skipeoln(file);
887  readtoeoln(file);
888 
889  /*
890  * READ H-BOND R**6 TERM FOR ALL N-B TYPES -BSOL()
891  */
892 
893  if (amber7_format)
894  if (!moveto(file,"HBOND_BCOEF"))
895  { genclose(file);
896  return 0;
897  }
898  for (i=0; i<Nphb; i++)
899 #ifdef DOUBLE
900  fscanf(file, " %lf", &HB6[i]);
901 #else
902  fscanf(file, " %f", &HB6[i]);
903 #endif
904 // skipeoln(file);
905  readtoeoln(file);
906 
907  /*
908  * READ H-BOND CUTOFF (NOT USED) ?? -HBCUT()
909  */
910 
911  if (amber7_format)
912  if (!moveto(file,"HBCUT"))
913  { genclose(file);
914  return 0;
915  }
916  H = (_REAL *) get(Nphb * sizeof(_REAL));
917  for (i=0; i<Nphb; i++)
918 #ifdef DOUBLE
919  fscanf(file, " %lf", &H[i]);
920 #else
921  fscanf(file, " %f", &H[i]);
922 #endif
923  free((char *)H);
924 
925 // skipeoln(file);
926  readtoeoln(file);
927 
928  /*
929  * READ ATOM SYMBOLS (FOR ANALYSIS PROGS) -IH(M06)
930  */
931 
932  if (amber7_format)
933  if (!moveto(file,"AMBER_ATOM_TYPE"))
934  { genclose(file);
935  return 0;
936  }
937  for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
938  preadln(file, "", &AtomSym[i*80]);
939 
940  /*
941  * READ TREE SYMBOLS (FOR ANALYSIS PROGS) -IH(M08)
942  */
943 
944  if (amber7_format)
945  if (!moveto(file,"TREE_CHAIN_CLASSIFICATION"))
946  { genclose(file);
947  return 0;
948  }
949  for (i=0; i<(Natom/20 + (Natom%20 ? 1 : 0)); i++)
950  preadln(file, "", &AtomTree[i*80]);
951 
952  /*
953  * READ TREE JOIN INFO (FOR ANALYSIS PROGS) -IX(I64)
954  */
955 
956  if (amber7_format)
957  { if (!moveto(file,"JOIN_ARRAY"))
958  { genclose(file);
959  return 0;
960  }
961  for (i=0; i<Natom; i++)
962  fscanf(file, " %d", &TreeJoin[i]);
963  }
964  else
965  { if (!read_fortran_12I6(file,TreeJoin,Natom))
966  { genclose(file);
967  return 0;
968  }
969 // skipeoln(file);
970  readtoeoln(file);
971  }
972 
973  /*
974  * READ PER-ATOM RES NUMBER -IX(I66)
975  * NOTE: this appears to be something entirely different
976  * NOTE: overwriting this with correct PER-ATOM RES NUMBERs
977  */
978 
979  if (amber7_format)
980  { if (!moveto(file,"IROTAT"))
981  { genclose(file);
982  return 0;
983  }
984  for (i=0; i<Natom; i++)
985  fscanf(file, " %d", &AtomRes[i]);
986  }
987  else
988  if (!read_fortran_12I6(file,AtomRes,Natom))
989  { genclose(file);
990  return 0;
991  }
992  res = 0;
993  for (i=0; i<Natom; i++) {
994  if (i+1 == Ipres[res+1]) /* atom is 1st of next res */
995  res++;
996  AtomRes[i] = res;
997  }
998 
999  /*
1000  * BOUNDARY CONDITION STUFF
1001  */
1002 
1003  if (!IfBox) {
1004  Nspm = 1;
1005  Boundary = (int *) get(sizeof(int)*Nspm);
1006  Boundary[0] = Natom;
1007  } else {
1008  if (amber7_format)
1009  { if (!moveto(file,"SOLVENT_POINTERS"))
1010  { genclose(file);
1011  return 0;
1012  }
1013  fscanf(file, " %d %d %d", &Iptres, &Nspm,
1014  &Nspsol);
1015  }
1016  else
1017  {
1018 // skipeoln(file);
1019  readtoeoln(file);
1020  buffer = new int[3];
1021  if (!read_fortran_12I6(file,buffer,3))
1022  { genclose(file);
1023  return 0;
1024  }
1025  Iptres = buffer[0];
1026  Nspm = buffer[1];
1027  Nspsol = buffer[2];
1028  delete [] buffer;
1029 // skipeoln(file);
1030  readtoeoln(file);
1031  }
1032  Boundary = (int *) get(sizeof(int)*Nspm);
1033  if (amber7_format)
1034  { if (!moveto(file,"ATOMS_PER_MOLECULE"))
1035  { genclose(file);
1036  return 0;
1037  }
1038  for (i=0; i<Nspm; i++)
1039  fscanf(file, " %d", &Boundary[i]);
1040  }
1041  else
1042  { if (!read_fortran_12I6(file,Boundary,Nspm))
1043  { genclose(file);
1044  return 0;
1045  }
1046 // skipeoln(file);
1047  readtoeoln(file);
1048  }
1049  if (amber7_format)
1050  if (!moveto(file,"BOX_DIMENSIONS"))
1051  { genclose(file);
1052  return 0;
1053  }
1054 #ifdef DOUBLE
1055  fscanf(file, " %lf %lf %lf",
1056 #else
1057  fscanf(file, " %f %f %f",
1058 #endif
1059  &Box[0], &Box[1], &Box[2]);
1060 // skipeoln(file);
1061  readtoeoln(file);
1062  if (Iptres)
1063  Ipatm = Ipres[Iptres] - 1;
1064  /* IF(IPTRES.GT.0) IPTATM = IX(I02+IPTRES-1+1)-1 */
1065  }
1066 
1067  /*
1068  * ----- LOAD THE CAP INFORMATION IF NEEDED -----
1069  */
1070 
1071  // I don't know the label for it, so I don't read it if
1072  // it's AMBER 7 format. It's not used in NAMD anyway.
1073 // if (IfCap) {
1074  if (IfCap && !amber7_format) {
1075  /* if (IfBox)
1076  skipeoln(file); */
1077 #ifdef DOUBLE
1078  fscanf(file, " %d %lf %lf %lf %lf",
1079 #else
1080  fscanf(file, " %d %f %f %f %f",
1081 #endif
1082  &Natcap, &Cutcap,
1083  &Xcap, &Ycap, &Zcap);
1084  }
1085  genclose(file);
1086  if (debug) {
1087  printf("rdprm done\n");
1088 #ifndef NAMD_NO_STDOUT_FLUSH
1089  fflush(stdout);
1090 #endif
1091  }
1092  data_read = 1;
1093  return(1);
1094 }
1095 
1097 {
1098  char *restr = ResNames;
1099  char *lastres = ResNames + Nres * 4 + 1;
1100  int res = 0;
1101 
1102  /*
1103  * find 1st water residue
1104  */
1105 
1106  for (; restr<lastres; restr+=4) {
1107  if (!strncmp(restr, "WAT ", 4)) {
1108  printf("first water: res = %d, atom = %d (%.4s)\n",
1109  res+1, Ipres[res],
1110  &AtomNames[Ipres[res]]);
1111 #ifndef NAMD_NO_STDOUT_FLUSH
1112  fflush(stdout);
1113 #endif
1114  return(Ipres[res]-1);
1115  }
1116  res++;
1117  }
1118  return(0);
1119 }
1120 
1121 
1122 // Constructer: simply set all the pointers Null
1123 
1125 {
1126  data_read = 0; // No data are read yet
1127  AtomNames = ResNames = AtomSym = AtomTree = NULL;
1128  Charges = Masses = Rk = Req = Tk = Teq = Pk = Pn = Phase = NULL;
1129  Solty = Cn1 = Cn2 = HB12 = HB6 = NULL;
1130  Iac = Iblo = Cno = Ipres = ExclAt = TreeJoin = AtomRes = NULL;
1131  BondHAt1 = BondHAt2 = BondHNum = BondAt1 = BondAt2 = NULL;
1133  AngleAt1 = AngleAt2 = AngleAt3 = AngleNum = DihHAt1 = NULL;
1134  DihHAt2 = DihHAt3 = DihHAt4 = DihHNum = DihAt1 = DihAt2 = NULL;
1135  DihAt3 = DihAt4 = DihNum = Boundary = NULL;
1136 }
1137 
1138 
1139 // Destructer: free all the allocated memory for arrays
1140 
1142 {
1143  free(AtomNames);
1144  free(Charges);
1145  free(Masses);
1146  free(Iac);
1147  free(Iblo);
1148  free(Cno);
1149  free(ResNames);
1150  free(Ipres);
1151  free(Rk);
1152  free(Req);
1153  free(Tk);
1154  free(Teq);
1155  free(Pk);
1156  free(Pn);
1157  free(Phase);
1158  free(Solty);
1159  free(Cn1);
1160  free(Cn2);
1161  free(BondHAt1);
1162  free(BondHAt2);
1163  free(BondHNum);
1164  free(BondAt1);
1165  free(BondAt2);
1166  free(BondNum);
1167  free(AngleHAt1);
1168  free(AngleHAt2);
1169  free(AngleHAt3);
1170  free(AngleHNum);
1171  free(AngleAt1);
1172  free(AngleAt2);
1173  free(AngleAt3);
1174  free(AngleNum);
1175  free(DihHAt1);
1176  free(DihHAt2);
1177  free(DihHAt3);
1178  free(DihHAt4);
1179  free(DihHNum);
1180  free(DihAt1);
1181  free(DihAt2);
1182  free(DihAt3);
1183  free(DihAt4);
1184  free(DihNum);
1185  free(ExclAt);
1186  free(HB12);
1187  free(HB6);
1188  free(AtomSym);
1189  free(AtomTree);
1190  free(TreeJoin);
1191  free(AtomRes);
1192 }
1193 
1194 
1195 // Read FORTRAN 12I6 format data (no space between adjacent data is assumed)
1196 // One needs to read the whole data block into the buffer here
1197 // fp - file pointer.
1198 // data - the buffer to hold the whole block. Should be allocated before the call.
1199 // count - number of ints in the block.
1200 
1201 int Ambertoppar::read_fortran_12I6(FILE *fp, int *data, int count)
1202 {
1203  int i,j;
1204  char buf[7];
1205 
1206  for (i=0; i<count; ++i)
1207  { for (j=0; j<6; ++j)
1208  { buf[j]=getc(fp);
1209  if (buf[j]=='\n' || buf[j]=='\0' || buf[j]==EOF)
1210  return 0;
1211  }
1212  buf[6] = '\0';
1213  if (sscanf(buf,"%d",data+i) != 1)
1214  return 0;
1215  if (i%12==11 && i<count-1)
1216  readtoeoln(fp);
1217  }
1218 
1219  return 1;
1220 }
1221 
1222 
1223 // Used for AMBER 7 format. The function moves the file postion to
1224 // the beginning of the data section labeled by "label"
1225 
1226 int Ambertoppar::moveto(FILE *fp, const char *label)
1227 {
1228  char s[76],buf[81];
1229 
1230  while ( 1 ) {
1231  // Find the string "%FLAG"
1232  do preadln(fp, "parm file", buf);
1233  while (strncmp(buf,"%FLAG",5));
1234 
1235  // See if the label is what we expected
1236  sscanf(buf+5,"%s",s);
1237  if (strcasecmp(s,label))
1238  iout << iWARN << "Skipping " << s << " in parm file while seeking " << label << ".\n" << endi;
1239  else
1240  break;
1241  }
1242 
1243  // The next line should begin with "%FORMAT"
1244  preadln(fp, "parm file", buf);
1245  if (strncmp(buf,"%FORMAT",7))
1246  return 0;
1247 
1248  return 1;
1249 }
_REAL * Pk
Definition: parm.h:24
int * DihHAt1
Definition: parm.h:27
_REAL Zcap
Definition: parm.h:26
int Nttyp
Definition: parm.h:17
_REAL Ycap
Definition: parm.h:26
_REAL * Teq
Definition: parm.h:24
char * AtomTree
Definition: parm.h:23
int * BondHNum
Definition: parm.h:27
void preadln(FILE *, const char *, char *)
Definition: parm.C:123
int * DihAt3
Definition: parm.h:27
char ititl[81]
Definition: parm.h:16
int Natom
Definition: parm.h:17
int IfBox
Definition: parm.h:17
_REAL * Phase
Definition: parm.h:24
int Natyp
Definition: parm.h:17
int Nphih
Definition: parm.h:17
int * AngleHNum
Definition: parm.h:27
int Iptres
Definition: parm.h:17
~parm()
Definition: parm.C:1141
int Mtheta
Definition: parm.h:17
void genclose(FILE *)
Definition: parm.C:86
int * BondAt2
Definition: parm.h:27
int Nbonh
Definition: parm.h:17
static int readtoeoln(FILE *f)
Definition: parm.C:51
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
static int debug
Definition: parm.C:36
int * DihHAt4
Definition: parm.h:27
_REAL * Tk
Definition: parm.h:24
#define iout
Definition: InfoStream.h:87
int Natcap
Definition: parm.h:17
int Mbona
Definition: parm.h:17
int Nhparm
Definition: parm.h:17
int readparm(char *)
Definition: parm.C:151
int * Iblo
Definition: parm.h:27
char * ResNames
Definition: parm.h:23
int * DihAt4
Definition: parm.h:27
int Mphia
Definition: parm.h:17
_REAL * Masses
Definition: parm.h:24
int * AngleHAt1
Definition: parm.h:27
int Numbnd
Definition: parm.h:17
_REAL * HB6
Definition: parm.h:24
_REAL * HB12
Definition: parm.h:24
int * AngleNum
Definition: parm.h:27
int * DihHNum
Definition: parm.h:27
int * Iac
Definition: parm.h:27
_REAL * Charges
Definition: parm.h:24
int Ntypes
Definition: parm.h:17
int * DihAt1
Definition: parm.h:27
int * BondNum
Definition: parm.h:27
int * DihNum
Definition: parm.h:27
_REAL Xcap
Definition: parm.h:26
int Nbona
Definition: parm.h:17
int Nmxrs
Definition: parm.h:17
int * BondHAt1
Definition: parm.h:27
int * AngleHAt3
Definition: parm.h:27
int firstwat()
Definition: parm.C:1096
char * AtomNames
Definition: parm.h:23
_REAL * Cn2
Definition: parm.h:24
FILE * Fopen(const char *filename, const char *mode)
Definition: common.C:265
int * Boundary
Definition: parm.h:27
int * AtomRes
Definition: parm.h:27
FILE * genopen(const char *name)
Definition: parm.C:72
int Ntheta
Definition: parm.h:17
int * AngleAt3
Definition: parm.h:27
int * AngleAt1
Definition: parm.h:27
int Nparm
Definition: parm.h:17
void NAMD_die(const char *err_msg)
Definition: common.C:83
char * AtomSym
Definition: parm.h:23
int * BondAt1
Definition: parm.h:27
int Ntype2d
Definition: parm.h:17
int Nspm
Definition: parm.h:17
int Ntheth
Definition: parm.h:17
_REAL * Rk
Definition: parm.h:24
Definition: Box.h:14
int data_read
Definition: parm.h:34
int Fclose(FILE *fout)
Definition: common.C:359
int Nnb
Definition: parm.h:17
int Nat3
Definition: parm.h:17
int IfCap
Definition: parm.h:17
int * DihHAt2
Definition: parm.h:27
int Nphb
Definition: parm.h:17
_REAL * Cn1
Definition: parm.h:24
int * TreeJoin
Definition: parm.h:27
parm()
Definition: parm.C:1124
int * DihAt2
Definition: parm.h:27
char * get(int)
Definition: parm.C:96
int Numang
Definition: parm.h:17
_REAL * Req
Definition: parm.h:24
int read_fortran_12I6(FILE *, int *, int)
Definition: parm.C:1201
int Ipatm
Definition: parm.h:17
_REAL * Solty
Definition: parm.h:24
int Nphia
Definition: parm.h:17
int * AngleAt2
Definition: parm.h:27
int moveto(FILE *, const char *)
Definition: parm.C:1226
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int * ExclAt
Definition: parm.h:27
int * AngleHAt2
Definition: parm.h:27
int Nres
Definition: parm.h:17
int * Ipres
Definition: parm.h:27
int Nptra
Definition: parm.h:17
#define _REAL
Definition: parm.h:11
int Nspsol
Definition: parm.h:17
int * BondHAt2
Definition: parm.h:27
int * Cno
Definition: parm.h:27
_REAL * Pn
Definition: parm.h:24
int * DihHAt3
Definition: parm.h:27
_REAL Cutcap
Definition: parm.h:26