1 /***************************************************************************
2 * RCS INFORMATION:
3 *
4 * $RCSfile: ReadPARM7.h,v $
5 * $Author: johns $ $Locker: $ $State: Exp $
6 * $Revision: 1.32 $ $Date: 2016/11/06 17:54:56 $
7 *
8 ***************************************************************************
9 * DESCRIPTION:
10 * NOTE:: Significant modifications were made to the VMD version of
11 * Bill Ross's original code in order to make it easy to hook
12 * into VMD plugin structures.
13 * Further modifications were made to the VMD code to
14 * read amber 7 parm files, courtesy of Brian Bennion
15 * Here is what has changed:
16 * Functions became Class Methods, data became instance variables
17 * The Code to check for compressed files before opening was disabled
18 * Methods get_parm7_atom, get_parm7_bond, get_hydrogen_bond,
19 * get_parm7_natoms, get_parm7_nbonds, get_parm7_boxInfo were added in
20 * order to convert from prm.c parlance to VMD conventions.
21 ***************************************************************************/
22
23 /*
24 * COPYRIGHT 1992, REGENTS OF THE UNIVERSITY OF CALIFORNIA
25 *
26 * prm.c - read information from an amber PARM topology file:
27 * atom/residue/bond/charge info, plus force field data.
28 * This file and the accompanying prm.h may be distributed
29 * provided this notice is retained unmodified and provided
30 * that any modifications to the rest of the file are noted
31 * in comments.
32 *
33 * Bill Ross, UCSF 1994
34 */
35
36 #ifndef READPARM7_H
37 #define READPARM7_H
38
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <ctype.h>
42 #include <sys/types.h>
43 #include <sys/stat.h>
44 #include <errno.h>
45 #include <string.h>
46 #include "molfile_plugin.h" // needed for molfile return codes etc
47
48 #if defined(WIN32) || defined(WIN64)
49 #define strcasecmp stricmp
50 #define strncasecmp strnicmp
51 #endif
52
53 #if 0
54 #define _REAL double
55 #define DBLFMT "%lf"
56 #else
57 #define _REAL float
58 #define DBLFMT "%f"
59 #endif
60
61
62 typedef struct parm {
63 char title[85];
64 char version[85];
65 int IfBox, Nmxrs, IfCap,
66 Natom, Ntypes, Nbonds, Nbonh, Mbona, Ntheth, Mtheta,
67 Nphih, Mphia, Nhparm, Nparm, Nnb, Nres,Mptra,
68 Nbona, Ntheta, Nphia, Numbnd, Numang, Nptra,Jparm,
69 Natyp, Nphb, Nat3, Ntype2d, Nttyp, Nspm, Iptres, Nspsol,
70 Ipatm, Natcap,Ifpert,Nbper,Ngper,Ndper,Mbper,Mgper,Mdper,
71 Numextra;
72 _REAL Box[3], Cutcap, Xcap, Ycap, Zcap;
73 } parmstruct;
74
read_parm7_flag(FILE * file,const char * flag,const char * format)75 static int read_parm7_flag(FILE *file, const char *flag, const char *format) {
76 char buf[1024];
77
78 /* read the %FLAG text */
79 fscanf(file, "%s\n", buf);
80 if (strcmp("%FLAG", buf)) {
81 printf("AMBER 7 parm read error, at flag section %s,\n", flag);
82 printf(" expected %%FLAG but got %s\n", buf);
83 return 0; /* read of flag data failed */
84 }
85
86 /* read field name specifier */
87 fscanf(file, "%s\n", buf);
88 if (flag != NULL) {
89 if (strcmp(flag, buf)) {
90 printf("AMBER 7 parm read error at flag section %s,\n", flag);
91 printf(" expected flag field %s but got %s\n", flag, buf);
92 return 0; /* read of flag data failed */
93 }
94 }
95
96 /* read format string */
97 fscanf(file, "%s\n", buf);
98 if (format != NULL) {
99 if (strcmp(format, buf)) {
100 if ( (!strcmp(flag, "TITLE") || !strcmp(flag, "CTITLE") )
101 && !strcmp(format, "%FORMAT(20a4)")) {
102 if (!strcmp(buf, "%FORMAT(a80)"))
103 return 1; /* accept a80 as substitute for 20a4 */
104 }
105 printf("AMBER 7 parm read error at flag section %s,\n", flag);
106 printf(" expected format %s but got %s\n", format, buf);
107 return 0; /* read of flag data failed */
108 }
109 }
110
111 return 1; /* read of flag data succeeded */
112 }
113
114 /*
115 * open_parm7_file() - fopen regular or popen compressed file for reading
116 * Return FILE handle on success.
117 * set as_pipe to 1 if opened with popen, or 0 if opened with fopen.
118 */
119
open_parm7_file(const char * name,int * as_pipe)120 static FILE *open_parm7_file(const char *name, int *as_pipe) {
121 struct stat buf;
122 char cbuf[8192];
123 int length;
124 int &compressed = *as_pipe;
125 FILE *fp;
126
127 length = strlen(name);
128 compressed = 0; // Just to start
129 strcpy(cbuf, name);
130
131 /*
132 * if file doesn't exist, maybe it has been compressed/decompressed
133 */
134 if (stat(cbuf, &buf) == -1) {
135 switch (errno) {
136 case ENOENT:
137 if (!compressed) {
138 strcat(cbuf, ".Z");
139 if (stat(cbuf, &buf) == -1) {
140 printf("%s, %s: does not exist\n", name, cbuf);
141 return(NULL);
142 }
143 compressed++;
144
145 // Don't modify the filename
146 //strcat(name, ".Z"); /* TODO: add protection */
147 } else {
148 cbuf[length-2] = '\0';
149 if (stat(cbuf, &buf) == -1) {
150 printf("%s, %s: does not exist\n", name, cbuf);
151 return(NULL);
152 }
153 compressed = 0;
154 }
155 break;
156
157 default:
158 return(NULL);
159 }
160 }
161
162 /*
163 * open the file
164 */
165 #if defined(_MSC_VER) || defined(__MINGW32__)
166 if (compressed) {
167 /* NO "zcat" on Win32 */
168 printf("Cannot load compressed PARM files on Windows.\n");
169 return NULL;
170 }
171 #else
172 if (compressed) {
173 char pcmd[sizeof(cbuf) + 7];
174 sprintf(pcmd, "zcat '%s'", cbuf);
175 if ((fp = popen(pcmd, "r")) == NULL) {
176 perror(pcmd);
177 return NULL;
178 }
179 }
180 #endif
181 else {
182 if ((fp = fopen(cbuf, "r")) == NULL) {
183 perror(cbuf);
184 return NULL;
185 }
186 }
187
188 return(fp);
189 }
190
191
parse_parm7_atoms(const char * fmt,int natoms,molfile_atom_t * atoms,FILE * file)192 static int parse_parm7_atoms(const char *fmt,
193 int natoms, molfile_atom_t *atoms, FILE *file) {
194 char buf[85];
195
196 if (strncasecmp(fmt, "%FORMAT(20a4)", 13))
197 return 0;
198
199 int j=0;
200 for (int i=0; i<natoms; i++) {
201 molfile_atom_t *atom = atoms+i;
202 if (!(i%20)) {
203 j=0;
204 fgets(buf, 85, file);
205 }
206 strncpy(atom->name, buf+4*j, 4);
207 atom->name[4]='\0';
208 j++;
209 }
210 return 1;
211 }
212
213
parse_parm7_charge(const char * fmt,int natoms,molfile_atom_t * atoms,FILE * file)214 static int parse_parm7_charge(const char *fmt,
215 int natoms, molfile_atom_t *atoms, FILE *file) {
216 if (strncasecmp(fmt, "%FORMAT(5E16.8)", 15) &&
217 strncasecmp(fmt, "%FORMAT(3E24.16)", 16))
218 return 0;
219
220 for (int i=0; i<natoms; i++) {
221 double q=0;
222 if (fscanf(file, " %lf", &q) != 1) {
223 fprintf(stderr, "PARM7: error reading charge at index %d\n", i);
224 return 0;
225 }
226 atoms[i].charge = 0.0548778 * (float)q; /* convert to elementary charge units */
227 }
228
229 return 1;
230 }
231
232
parse_parm7_mass(const char * fmt,int natoms,molfile_atom_t * atoms,FILE * file)233 static int parse_parm7_mass(const char *fmt,
234 int natoms, molfile_atom_t *atoms, FILE *file) {
235 if (strncasecmp(fmt, "%FORMAT(5E16.8)", 15)) return 0;
236 for (int i=0; i<natoms; i++) {
237 double m=0;
238 if (fscanf(file, " %lf", &m) != 1) {
239 fprintf(stderr, "PARM7: error reading mass at index %d\n", i);
240 return 0;
241 }
242 atoms[i].mass = (float)m;
243 }
244 return 1;
245 }
246
247
parse_parm7_atype(const char * fmt,int natoms,molfile_atom_t * atoms,FILE * file)248 static int parse_parm7_atype(const char *fmt,
249 int natoms, molfile_atom_t *atoms, FILE *file) {
250 if (strncasecmp(fmt, "%FORMAT(20a4)", 13)) return 0;
251 char buf[85];
252 int j=0;
253 for (int i=0; i<natoms; i++) {
254 molfile_atom_t *atom = atoms+i;
255 if (!(i%20)) {
256 j=0;
257 fgets(buf, 85, file);
258 }
259 strncpy(atom->type, buf+4*j, 4);
260 atom->type[4]='\0';
261 j++;
262 }
263 return 1;
264 }
265
266
parse_parm7_resnames(const char * fmt,int nres,char * resnames,FILE * file)267 static int parse_parm7_resnames(const char *fmt,
268 int nres, char *resnames, FILE *file) {
269 if (strncasecmp(fmt, "%FORMAT(20a4)", 13)) return 0;
270 char buf[85];
271 int j=0;
272 for (int i=0; i<nres; i++) {
273 if (!(i%20)) {
274 j=0;
275 fgets(buf, 85, file);
276 }
277 strncpy(resnames, buf+4*j, 4);
278 resnames += 4;
279 j++;
280 }
281
282 return 1;
283 }
284
parse_parm7_respointers(const char * fmt,int natoms,molfile_atom_t * atoms,int nres,const char * resnames,FILE * file)285 static int parse_parm7_respointers(const char *fmt, int natoms,
286 molfile_atom_t *atoms, int nres, const char *resnames, FILE *file) {
287 if (strncasecmp(fmt, "%FORMAT(10I8)", 13)) return 0;
288 int cur, next;
289 fscanf(file, " %d", &cur);
290 for (int i=1; i<nres; i++) {
291 if (fscanf(file, " %d", &next) != 1) {
292 fprintf(stderr, "PARM7: error reading respointer records at residue %d\n", i);
293 return 0;
294 }
295 while (cur < next) {
296 if (cur > natoms) {
297 fprintf(stderr, "invalid atom index: %d\n", cur);
298 return 0;
299 }
300 strncpy(atoms[cur-1].resname, resnames, 4);
301 atoms[cur-1].resname[4] = '\0';
302 atoms[cur-1].resid = i;
303 cur++;
304 }
305 resnames += 4;
306 }
307 // store the last residue name
308 while (cur <= natoms) {
309 strncpy(atoms[cur-1].resname, resnames, 4);
310 atoms[cur-1].resname[4] = '\0';
311 atoms[cur-1].resid = nres;
312 cur++;
313 }
314
315 return 1;
316 }
317
318
parse_parm7_bonds(const char * fmt,int nbonds,int * from,int * to,FILE * file)319 static int parse_parm7_bonds(const char *fmt,
320 int nbonds, int *from, int *to, FILE *file) {
321 if (strncasecmp(fmt, "%FORMAT(10I8)", 13))
322 return 0;
323
324 int a, b, tmp;
325 for (int i=0; i<nbonds; i++) {
326 if (fscanf(file, " %d %d %d", &a, &b, &tmp) != 3) {
327 fprintf(stderr, "PARM7: error reading bond number %d\n", i);
328 return 0;
329 }
330 from[i] = a/3 + 1;
331 to[i] = b/3 + 1;
332 }
333
334 return 1;
335 }
336
337
338 /*
339 * close_parm7_file() - close fopened or popened file
340 */
close_parm7_file(FILE * fileptr,int popn)341 static void close_parm7_file(FILE *fileptr, int popn) {
342 #if defined(_MSC_VER) || defined(__MINGW32__)
343 if (popn) {
344 printf("pclose() no such function on win32!\n");
345 } else {
346 if (fclose(fileptr) == -1)
347 perror("fclose");
348 }
349 #else
350 if (popn) {
351 if (pclose(fileptr) == -1)
352 perror("pclose");
353 } else {
354 if (fclose(fileptr) == -1)
355 perror("fclose");
356 }
357 #endif
358 }
359
360 static const char *parm7 = "%8d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n";
361
read_parm7_header(FILE * file)362 static parmstruct *read_parm7_header(FILE *file) {
363 char sdum[512];
364 parmstruct *prm;
365 prm = new parmstruct;
366
367 /* READ VERSION */
368 fgets(sdum, 512, file);
369
370 /* READ TITLE */
371 fscanf(file, "%s\n", sdum); // "%FLAG"
372 if (strcmp("%FLAG", sdum)) {
373 printf("AMBER 7 parm read error, can't find TITLE flag.\n");
374 printf(" expected %%FLAG, got %s\n", sdum);
375 delete prm;
376 return NULL;
377 }
378 fscanf(file, "%s\n", sdum); // "TITLE" or "CTITLE"
379 if (strcmp("TITLE", sdum) && strcmp("CTITLE", sdum)) {
380 printf("AMBER 7 parm read error, at flag section TITLE,\n");
381 printf(" expected TITLE or CTITLE but got %s,\n", sdum);
382 delete prm;
383 return NULL;
384 }
385 fscanf(file, "%s\n", sdum); // "FORMAT (20a4)"
386 if (strcmp(sdum, "%FORMAT(20a4)") && strcmp(sdum, "%FORMAT(a80)")) {
387 printf("AMBER 7 parm read error, at flag section TITLE,\n");
388 printf(" expected %%FLAG but got %s,\n", sdum);
389 delete prm;
390 return NULL;
391 }
392
393 // read the title string itself, and handle empty lines
394 #if 0
395 // XXX this code fails with some AMBER 9 test files
396 fscanf(file, "%s\n", prm->title);
397 #else
398 // XXX this hack causes AMBER 9 prmtop files to load
399 fgets(prm->title, sizeof(prm->title), file);
400 #endif
401
402 if (strstr(prm->title, "%FLAG") == NULL) {
403 // Got a title string
404 if (!read_parm7_flag(file, "POINTERS", "%FORMAT(10I8)")) {
405 delete prm;
406 return NULL;
407 }
408 } else {
409 // NO title string, use a special method to pick up next flag
410 #if 0
411 fscanf(file,"%s\n", sdum);
412 if (strcmp("POINTERS", sdum)) {
413 printf("AMBER 7 parm read error at flag section POINTERS\n");
414 printf(" expected flag field POINTERS but got %s\n", sdum);
415 #else
416 if (strstr(prm->title, "POINTERS") == NULL) {
417 printf("AMBER 7 parm read error at flag section POINTERS\n");
418 printf(" expected flag field POINTERS but got %s\n", prm->title);
419 #endif
420 delete prm;
421 return NULL;
422 }
423 #if 0
424 fscanf(file,"%s\n", sdum);
425 if (strcasecmp("%FORMAT(10I8)", sdum)) {
426 #else
427 fgets(sdum, sizeof(sdum), file);
428 if ((strstr(sdum, "%FORMAT(10I8)") == NULL) &&
429 (strstr(sdum, "%FORMAT(10i8)") == NULL)) {
430 #endif
431 printf("AMBER 7 parm read error at flag section POINTERS,\n");
432 printf(" expected format %%FORMAT(10I8) but got %s\n", sdum);
433 delete prm;
434 return NULL;
435 }
436 }
437
438 /* READ POINTERS (CONTROL INTEGERS) */
439 fscanf(file,parm7,
440 &prm->Natom, &prm->Ntypes, &prm->Nbonh, &prm->Nbona,
441 &prm->Ntheth, &prm->Ntheta, &prm->Nphih, &prm->Nphia,
442 &prm->Jparm, &prm->Nparm);
443 fscanf(file, parm7,
444 &prm->Nnb, &prm->Nres, &prm->Mbona, &prm->Mtheta,
445 &prm->Mphia, &prm->Numbnd, &prm->Numang, &prm->Mptra,
446 &prm->Natyp, &prm->Nphb);
447 fscanf(file, parm7, &prm->Ifpert, &prm->Nbper, &prm->Ngper,
448 &prm->Ndper, &prm->Mbper, &prm->Mgper, &prm->Mdper,
449 &prm->IfBox, &prm->Nmxrs, &prm->IfCap);
450
451 fscanf(file,"%8d",&prm->Numextra); //BB
452 prm->Nptra=prm->Mptra; //BB new to amber 7 files...
453
454 prm->Nat3 = 3 * prm->Natom;
455 prm->Ntype2d = prm->Ntypes * prm->Ntypes;
456 prm->Nttyp = prm->Ntypes*(prm->Ntypes+1)/2;
457
458 return prm;
459 }
460
461
462 #endif
463