?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed bug 553 by checking whether an atom name is longer than 2 characters
[gromacs/rigid-bodies.git] / src / gmxlib / pdbio.c
blob33342dd1b5c1328dce2a30b56714d45b5f4f2777
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * GROningen Mixture of Alchemy and Childrens' Stories
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <ctype.h>
40
41 #include "sysstuff.h"
42 #include "string2.h"
43 #include "vec.h"
44 #include "smalloc.h"
45 #include "typedefs.h"
46 #include "symtab.h"
47 #include "pdbio.h"
48 #include "vec.h"
49 #include "copyrite.h"
50 #include "futil.h"
51 #include "atomprop.h"
52 #include "physics.h"
53 #include "pbc.h"
54 #include "gmxfio.h"
55
56 typedef struct {
57 int ai,aj;
58 } gmx_conection_t;
59
60 typedef struct gmx_conect_t {
61 int nconect;
62 gmx_bool bSorted;
63 gmx_conection_t *conect;
64 } gmx_conect_t;
65
66 static const char *pdbtp[epdbNR]={
67 "ATOM ","HETATM", "ANISOU", "CRYST1",
68 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
69 "CONECT"
70 };
71
72
73 /* this is not very good,
74 but these are only used in gmx_trjconv and gmx_editconv */
75 static gmx_bool bWideFormat=FALSE;
76 #define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
77
78 void set_pdb_wide_format(gmx_bool bSet)
79 {
80 bWideFormat = bSet;
81 }
82
83 static void xlate_atomname_pdb2gmx(char *name)
84 {
85 int i,length;
86 char temp;
87
88 length=strlen(name);
89 if (length>3 && isdigit(name[0])) {
90 temp=name[0];
91 for(i=1; i<length; i++)
92 name[i-1]=name[i];
93 name[length-1]=temp;
94 }
95 }
96
97 static void xlate_atomname_gmx2pdb(char *name)
98 {
99 int i,length;
100 char temp;
101
102 length=strlen(name);
103 if (length>3 && isdigit(name[length-1])) {
104 temp=name[length-1];
105 for(i=length-1; i>0; --i)
106 name[i]=name[i-1];
107 name[0]=temp;
108 }
109 }
110
111
112 void gmx_write_pdb_box(FILE *out,int ePBC,matrix box)
113 {
114 real alpha,beta,gamma;
115
116 if (ePBC == -1)
117 ePBC = guess_ePBC(box);
118
119 if (ePBC == epbcNONE)
120 return;
121
122 if (norm2(box[YY])*norm2(box[ZZ])!=0)
123 alpha = RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ]));
124 else
125 alpha = 90;
126 if (norm2(box[XX])*norm2(box[ZZ])!=0)
127 beta = RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ]));
128 else
129 beta = 90;
130 if (norm2(box[XX])*norm2(box[YY])!=0)
131 gamma = RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY]));
132 else
133 gamma = 90;
134 fprintf(out,"REMARK THIS IS A SIMULATION BOX\n");
135 if (ePBC != epbcSCREW) {
136 fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
137 10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
138 alpha,beta,gamma,"P 1",1);
139 } else {
140 /* Double the a-vector length and write the correct space group */
141 fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
142 20*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
143 alpha,beta,gamma,"P 21 1 1",1);
144
145 }
146 }
147
148 static void read_cryst1(char *line,int *ePBC,matrix box)
149 {
150 #define SG_SIZE 11
151 char sa[12],sb[12],sc[12],sg[SG_SIZE+1],ident;
152 double fa,fb,fc,alpha,beta,gamma,cosa,cosb,cosg,sing;
153 int syma,symb,symc;
154 int ePBC_file;
155
156 sscanf(line,"%*s%s%s%s%lf%lf%lf",sa,sb,sc,&alpha,&beta,&gamma);
157
158 ePBC_file = -1;
159 if (strlen(line) >= 55) {
160 strncpy(sg,line+55,SG_SIZE);
161 sg[SG_SIZE] = '\0';
162 ident = ' ';
163 syma = 0;
164 symb = 0;
165 symc = 0;
166 sscanf(sg,"%c %d %d %d",&ident,&syma,&symb,&symc);
167 if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1) {
168 fc = strtod(sc,NULL)*0.1;
169 ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
170 }
171 if (ident == 'P' && syma == 21 && symb == 1 && symc == 1) {
172 ePBC_file = epbcSCREW;
173 }
174 }
175 if (ePBC) {
176 *ePBC = ePBC_file;
177 }
178
179 if (box) {
180 fa = strtod(sa,NULL)*0.1;
181 fb = strtod(sb,NULL)*0.1;
182 fc = strtod(sc,NULL)*0.1;
183 if (ePBC_file == epbcSCREW) {
184 fa *= 0.5;
185 }
186 clear_mat(box);
187 box[XX][XX] = fa;
188 if ((alpha!=90.0) || (beta!=90.0) || (gamma!=90.0)) {
189 if (alpha != 90.0) {
190 cosa = cos(alpha*DEG2RAD);
191 } else {
192 cosa = 0;
193 }
194 if (beta != 90.0) {
195 cosb = cos(beta*DEG2RAD);
196 } else {
197 cosb = 0;
198 }
199 if (gamma != 90.0) {
200 cosg = cos(gamma*DEG2RAD);
201 sing = sin(gamma*DEG2RAD);
202 } else {
203 cosg = 0;
204 sing = 1;
205 }
206 box[YY][XX] = fb*cosg;
207 box[YY][YY] = fb*sing;
208 box[ZZ][XX] = fc*cosb;
209 box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
210 box[ZZ][ZZ] = sqrt(fc*fc
211 - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
212 } else {
213 box[YY][YY] = fb;
214 box[ZZ][ZZ] = fc;
215 }
216 }
217 }
218
219 void write_pdbfile_indexed(FILE *out,const char *title,
220 t_atoms *atoms,rvec x[],
221 int ePBC,matrix box,char chainid,
222 int model_nr, atom_id nindex, atom_id index[],
223 gmx_conect conect, gmx_bool bTerSepChains)
224 {
225 gmx_conect_t *gc = (gmx_conect_t *)conect;
226 char resnm[6],nm[6],pdbform[128],pukestring[100];
227 atom_id i,ii;
228 int resind,resnr,type;
229 unsigned char resic,ch;
230 real occup,bfac;
231 gmx_bool bOccup;
232 int nlongname=0;
233 int chainnum,lastchainnum;
234 int lastresind,lastchainresind;
235 gmx_residuetype_t rt;
236 const char *p_restype;
237 const char *p_lastrestype;
238
239 gmx_residuetype_init(&rt);
240
241 bromacs(pukestring,99);
242 fprintf(out,"TITLE %s\n",(title && title[0])?title:pukestring);
243 if (bWideFormat) {
244 fprintf(out,"REMARK This file does not adhere to the PDB standard\n");
245 fprintf(out,"REMARK As a result of, some programs may not like it\n");
246 }
247 if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) ) {
248 gmx_write_pdb_box(out,ePBC,box);
249 }
250 if (atoms->pdbinfo) {
251 /* Check whether any occupancies are set, in that case leave it as is,
252 * otherwise set them all to one
253 */
254 bOccup = TRUE;
255 for (ii=0; (ii<nindex) && bOccup; ii++) {
256 i = index[ii];
257 bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
258 }
259 }
260 else
261 bOccup = FALSE;
262
263 fprintf(out,"MODEL %8d\n",model_nr>=0 ? model_nr : 1);
264
265 lastchainresind = -1;
266 lastchainnum = -1;
267 resind = -1;
268 p_restype = NULL;
269
270 for (ii=0; ii<nindex; ii++) {
271 i=index[ii];
272 lastresind = resind;
273 resind = atoms->atom[i].resind;
274 chainnum = atoms->resinfo[resind].chainnum;
275 p_lastrestype = p_restype;
276 gmx_residuetype_get_type(rt,*atoms->resinfo[resind].name,&p_restype);
277
278 /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
279 if( bTerSepChains && ii>0 && chainnum != lastchainnum)
280 {
281 /* Only add TER if the previous chain contained protein/DNA/RNA. */
282 if(gmx_residuetype_is_protein(rt,p_lastrestype) || gmx_residuetype_is_dna(rt,p_lastrestype) || gmx_residuetype_is_rna(rt,p_lastrestype))
283 {
284 fprintf(out,"TER\n");
285 }
286 lastchainnum = chainnum;
287 lastchainresind = lastresind;
288 }
289
290 strncpy(resnm,*atoms->resinfo[resind].name,sizeof(resnm)-1);
291 strncpy(nm,*atoms->atomname[i],sizeof(nm)-1);
292 /* rename HG12 to 2HG1, etc. */
293 xlate_atomname_gmx2pdb(nm);
294 resnr = atoms->resinfo[resind].nr;
295 resic = atoms->resinfo[resind].ic;
296 if (chainid!=' ')
297 {
298 ch = chainid;
299 }
300 else
301 {
302 ch = atoms->resinfo[resind].chainid;
303
304 if (ch == 0)
305 {
306 ch = ' ';
307 }
308 }
309 if (resnr>=10000)
310 resnr = resnr % 10000;
311 if (atoms->pdbinfo) {
312 type = atoms->pdbinfo[i].type;
313 occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
314 bfac = atoms->pdbinfo[i].bfac;
315 }
316 else {
317 type = 0;
318 occup = 1.0;
319 bfac = 0.0;
320 }
321 if (bWideFormat)
322 strcpy(pdbform,
323 "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
324 else {
325 /* Check whether atomname is an element name */
326 if ((strlen(nm)<4) && (gmx_strcasecmp(nm,atoms->atom[i].elem) != 0))
327 strcpy(pdbform,pdbformat);
328 else {
329 strcpy(pdbform,pdbformat4);
330 if (strlen(nm) > 4) {
331 int maxwln=20;
332 if (nlongname < maxwln) {
333 fprintf(stderr,"WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n",nm);
334 } else if (nlongname == maxwln) {
335 fprintf(stderr,"WARNING: More than %d long atom names, will not write more warnings\n",maxwln);
336 }
337 nlongname++;
338 }
339 }
340 strcat(pdbform,"%6.2f%6.2f %2s\n");
341 }
342 fprintf(out,pdbform,pdbtp[type],(i+1)%100000,nm,resnm,ch,resnr,
343 (resic == '\0') ? ' ' : resic,
344 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],occup,bfac,atoms->atom[i].elem);
345 if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic) {
346 fprintf(out,"ANISOU%5u %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
347 (i+1)%100000,nm,resnm,ch,resnr,
348 (resic == '\0') ? ' ' : resic,
349 atoms->pdbinfo[i].uij[0],atoms->pdbinfo[i].uij[1],
350 atoms->pdbinfo[i].uij[2],atoms->pdbinfo[i].uij[3],
351 atoms->pdbinfo[i].uij[4],atoms->pdbinfo[i].uij[5]);
352 }
353 }
354
355 fprintf(out,"TER\n");
356 fprintf(out,"ENDMDL\n");
357
358 if (NULL != gc) {
359 /* Write conect records */
360 for(i=0; (i<gc->nconect); i++) {
361 fprintf(out,"CONECT%5d%5d\n",gc->conect[i].ai+1,gc->conect[i].aj+1);
362 }
363 }
364 }
365
366 void write_pdbfile(FILE *out,const char *title, t_atoms *atoms,rvec x[],
367 int ePBC,matrix box,char chainid,int model_nr,gmx_conect conect,gmx_bool bTerSepChains)
368 {
369 atom_id i,*index;
370
371 snew(index,atoms->nr);
372 for(i=0; i<atoms->nr; i++)
373 index[i]=i;
374 write_pdbfile_indexed(out,title,atoms,x,ePBC,box,chainid,model_nr,
375 atoms->nr,index,conect,bTerSepChains);
376 sfree(index);
377 }
378
379 int line2type(char *line)
380 {
381 int k;
382 char type[8];
383
384 for(k=0; (k<6); k++)
385 type[k]=line[k];
386 type[k]='\0';
387
388 for(k=0; (k<epdbNR); k++)
389 if (strncmp(type,pdbtp[k],strlen(pdbtp[k])) == 0)
390 break;
391
392 return k;
393 }
394
395 static void read_anisou(char line[],int natom,t_atoms *atoms)
396 {
397 int i,j,k,atomnr;
398 char nc='\0';
399 char anr[12],anm[12];
400
401 /* Skip over type */
402 j=6;
403 for(k=0; (k<5); k++,j++) anr[k]=line[j];
404 anr[k]=nc;
405 j++;
406 for(k=0; (k<4); k++,j++) anm[k]=line[j];
407 anm[k]=nc;
408 j++;
409
410 /* Strip off spaces */
411 trim(anm);
412
413 /* Search backwards for number and name only */
414 atomnr = strtol(anr, NULL, 10);
415 for(i=natom-1; (i>=0); i--)
416 if ((strcmp(anm,*(atoms->atomname[i])) == 0) &&
417 (atomnr == atoms->pdbinfo[i].atomnr))
418 break;
419 if (i < 0)
420 fprintf(stderr,"Skipping ANISOU record (atom %s %d not found)\n",
421 anm,atomnr);
422 else {
423 if (sscanf(line+29,"%d%d%d%d%d%d",
424 &atoms->pdbinfo[i].uij[U11],&atoms->pdbinfo[i].uij[U22],
425 &atoms->pdbinfo[i].uij[U33],&atoms->pdbinfo[i].uij[U12],
426 &atoms->pdbinfo[i].uij[U13],&atoms->pdbinfo[i].uij[U23])
427 == 6) {
428 atoms->pdbinfo[i].bAnisotropic = TRUE;
429 }
430 else {
431 fprintf(stderr,"Invalid ANISOU record for atom %d\n",i);
432 atoms->pdbinfo[i].bAnisotropic = FALSE;
433 }
434 }
435 }
436
437 void get_pdb_atomnumber(t_atoms *atoms,gmx_atomprop_t aps)
438 {
439 int i,atomnumber,len;
440 size_t k;
441 char anm[6],anm_copy[6];
442 char nc='\0';
443 real eval;
444
445 if (!atoms->pdbinfo) {
446 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
447 }
448 for(i=0; (i<atoms->nr); i++) {
449 strcpy(anm,atoms->pdbinfo[i].atomnm);
450 strcpy(anm_copy,atoms->pdbinfo[i].atomnm);
451 len = strlen(anm);
452 atomnumber = NOTSET;
453 if ((anm[0] != ' ') && ((len <=2) || ((len > 2) && !isdigit(anm[2])))) {
454 anm_copy[2] = nc;
455 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
456 atomnumber = gmx_nint(eval);
457 else {
458 anm_copy[1] = nc;
459 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
460 atomnumber = gmx_nint(eval);
461 }
462 }
463 if (atomnumber == NOTSET) {
464 k=0;
465 while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
466 k++;
467 anm_copy[0] = anm[k];
468 anm_copy[1] = nc;
469 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
470 atomnumber = gmx_nint(eval);
471 }
472 atoms->atom[i].atomnumber = atomnumber;
473 sprintf(atoms->atom[i].elem,"%2s",gmx_atomprop_element(aps,atomnumber));
474 if (debug)
475 fprintf(debug,"Atomnumber for atom '%s' is %d\n",anm,atomnumber);
476 }
477 }
478
479 static int read_atom(t_symtab *symtab,
480 char line[],int type,int natom,
481 t_atoms *atoms,rvec x[],int chainnum,gmx_bool bChange)
482 {
483 t_atom *atomn;
484 int j,k;
485 char nc='\0';
486 char anr[12],anm[12],anm_copy[12],altloc,resnm[12],rnr[12];
487 char xc[12],yc[12],zc[12],occup[12],bfac[12];
488 unsigned char resic;
489 char chainid;
490 int resnr,atomnumber;
491
492 if (natom>=atoms->nr)
493 gmx_fatal(FARGS,"\nFound more atoms (%d) in pdb file than expected (%d)",
494 natom+1,atoms->nr);
495
496 /* Skip over type */
497 j=6;
498 for(k=0; (k<5); k++,j++) anr[k]=line[j];
499 anr[k]=nc;
500 trim(anr);
501 j++;
502 for(k=0; (k<4); k++,j++) anm[k]=line[j];
503 anm[k]=nc;
504 strcpy(anm_copy,anm);
505 atomnumber = NOTSET;
506 trim(anm);
507 altloc=line[j];
508 j++;
509 for(k=0; (k<4); k++,j++)
510 resnm[k]=line[j];
511 resnm[k]=nc;
512 trim(resnm);
513
514 chainid = line[j];
515 j++;
516
517 for(k=0; (k<4); k++,j++) {
518 rnr[k] = line[j];
519 }
520 rnr[k] = nc;
521 trim(rnr);
522 resnr = strtol(rnr, NULL, 10);
523 resic = line[j];
524 j+=4;
525
526 /* X,Y,Z Coordinate */
527 for(k=0; (k<8); k++,j++) xc[k]=line[j];
528 xc[k]=nc;
529 for(k=0; (k<8); k++,j++) yc[k]=line[j];
530 yc[k]=nc;
531 for(k=0; (k<8); k++,j++) zc[k]=line[j];
532 zc[k]=nc;
533
534 /* Weight */
535 for(k=0; (k<6); k++,j++) occup[k]=line[j];
536 occup[k]=nc;
537
538 /* B-Factor */
539 for(k=0; (k<7); k++,j++) bfac[k]=line[j];
540 bfac[k]=nc;
541
542 if (atoms->atom) {
543 atomn=&(atoms->atom[natom]);
544 if ((natom==0) ||
545 atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
546 atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
547 (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name,resnm) != 0))
548 {
549 if (natom == 0) {
550 atomn->resind = 0;
551 } else {
552 atomn->resind = atoms->atom[natom-1].resind + 1;
553 }
554 atoms->nres = atomn->resind + 1;
555 t_atoms_set_resinfo(atoms,natom,symtab,resnm,resnr,resic,chainnum,chainid);
556 }
557 else
558 {
559 atomn->resind = atoms->atom[natom-1].resind;
560 }
561 if (bChange) {
562 xlate_atomname_pdb2gmx(anm);
563 }
564 atoms->atomname[natom]=put_symtab(symtab,anm);
565 atomn->m = 0.0;
566 atomn->q = 0.0;
567 atomn->atomnumber = atomnumber;
568 atomn->elem[0] = '\0';
569 }
570 x[natom][XX]=strtod(xc,NULL)*0.1;
571 x[natom][YY]=strtod(yc,NULL)*0.1;
572 x[natom][ZZ]=strtod(zc,NULL)*0.1;
573 if (atoms->pdbinfo) {
574 atoms->pdbinfo[natom].type=type;
575 atoms->pdbinfo[natom].atomnr=strtol(anr, NULL, 10);
576 atoms->pdbinfo[natom].altloc=altloc;
577 strcpy(atoms->pdbinfo[natom].atomnm,anm_copy);
578 atoms->pdbinfo[natom].bfac=strtod(bfac,NULL);
579 atoms->pdbinfo[natom].occup=strtod(occup,NULL);
580 }
581 natom++;
582
583 return natom;
584 }
585
586 gmx_bool is_hydrogen(const char *nm)
587 {
588 char buf[30];
589
590 strcpy(buf,nm);
591 trim(buf);
592
593 if (buf[0] == 'H')
594 return TRUE;
595 else if ((isdigit(buf[0])) && (buf[1] == 'H'))
596 return TRUE;
597 return FALSE;
598 }
599
600 gmx_bool is_dummymass(const char *nm)
601 {
602 char buf[30];
603
604 strcpy(buf,nm);
605 trim(buf);
606
607 if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
608 return TRUE;
609
610 return FALSE;
611 }
612
613 static void gmx_conect_addline(gmx_conect_t *con,char *line)
614 {
615 int n,ai,aj;
616 char format[32],form2[32];
617
618 sprintf(form2,"%%*s");
619 sprintf(format,"%s%%d",form2);
620 if (sscanf(line,format,&ai) == 1) {
621 do {
622 strcat(form2,"%*s");
623 sprintf(format,"%s%%d",form2);
624 n = sscanf(line,format,&aj);
625 if (n == 1) {
626 srenew(con->conect,++con->nconect);
627 con->conect[con->nconect-1].ai = ai-1;
628 con->conect[con->nconect-1].aj = aj-1;
629 }
630 } while (n == 1);
631 }
632 }
633
634 void gmx_conect_dump(FILE *fp,gmx_conect conect)
635 {
636 gmx_conect_t *gc = (gmx_conect_t *)conect;
637 int i;
638
639 for(i=0; (i<gc->nconect); i++)
640 fprintf(fp,"%6s%5d%5d\n","CONECT",
641 gc->conect[i].ai+1,gc->conect[i].aj+1);
642 }
643
644 gmx_conect gmx_conect_init()
645 {
646 gmx_conect_t *gc;
647
648 snew(gc,1);
649
650 return (gmx_conect) gc;
651 }
652
653 void gmx_conect_done(gmx_conect conect)
654 {
655 gmx_conect_t *gc = (gmx_conect_t *)conect;
656
657 sfree(gc->conect);
658 }
659
660 gmx_bool gmx_conect_exist(gmx_conect conect,int ai,int aj)
661 {
662 gmx_conect_t *gc = (gmx_conect_t *)conect;
663 int i;
664
665 /* if (!gc->bSorted)
666 sort_conect(gc);*/
667
668 for(i=0; (i<gc->nconect); i++)
669 if (((gc->conect[i].ai == ai) &&
670 (gc->conect[i].aj == aj)) ||
671 ((gc->conect[i].aj == ai) &&
672 (gc->conect[i].ai == aj)))
673 return TRUE;
674 return FALSE;
675 }
676
677 void gmx_conect_add(gmx_conect conect,int ai,int aj)
678 {
679 gmx_conect_t *gc = (gmx_conect_t *)conect;
680 int i;
681
682 /* if (!gc->bSorted)
683 sort_conect(gc);*/
684
685 if (!gmx_conect_exist(conect,ai,aj)) {
686 srenew(gc->conect,++gc->nconect);
687 gc->conect[gc->nconect-1].ai = ai;
688 gc->conect[gc->nconect-1].aj = aj;
689 }
690 }
691
692 int read_pdbfile(FILE *in,char *title,int *model_nr,
693 t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
694 gmx_conect conect)
695 {
696 gmx_conect_t *gc = (gmx_conect_t *)conect;
697 t_symtab symtab;
698 gmx_bool bCOMPND;
699 gmx_bool bConnWarn = FALSE;
700 char line[STRLEN+1];
701 int line_type;
702 char *c,*d;
703 int natom,chainnum,nres_ter_prev=0;
704 char chidmax=' ';
705 gmx_bool bStop=FALSE;
706
707 if (ePBC)
708 {
709 /* Only assume pbc when there is a CRYST1 entry */
710 *ePBC = epbcNONE;
711 }
712 if (box != NULL)
713 {
714 clear_mat(box);
715 }
716
717 open_symtab(&symtab);
718
719 bCOMPND=FALSE;
720 title[0]='\0';
721 natom=0;
722 chainnum=0;
723 while (!bStop && (fgets2(line,STRLEN,in) != NULL))
724 {
725 line_type = line2type(line);
726
727 switch(line_type)
728 {
729 case epdbATOM:
730 case epdbHETATM:
731 natom = read_atom(&symtab,line,line_type,natom,atoms,x,chainnum,bChange);
732 break;
733
734 case epdbANISOU:
735 if (atoms->pdbinfo)
736 {
737 read_anisou(line,natom,atoms);
738 }
739 break;
740
741 case epdbCRYST1:
742 read_cryst1(line,ePBC,box);
743 break;
744
745 case epdbTITLE:
746 case epdbHEADER:
747 if (strlen(line) > 6)
748 {
749 c=line+6;
750 /* skip HEADER or TITLE and spaces */
751 while (c && (c[0]!=' ')) c++;
752 while (c && (c[0]==' ')) c++;
753 /* truncate after title */
754 d=strstr(c," ");
755 if (d)
756 {
757 d[0]='\0';
758 }
759 if (strlen(c)>0)
760 {
761 strcpy(title,c);
762 }
763 }
764 break;
765
766 case epdbCOMPND:
767 if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:")))
768 {
769 if ( !(c=strstr(line+6,"MOLECULE:")) )
770 {
771 c=line;
772 }
773 /* skip 'MOLECULE:' and spaces */
774 while (c && (c[0]!=' ')) c++;
775 while (c && (c[0]==' ')) c++;
776 /* truncate after title */
777 d=strstr(c," ");
778 if (d)
779 {
780 while ( (d[-1]==';') && d>c) d--;
781 d[0]='\0';
782 }
783 if (strlen(c) > 0)
784 {
785 if (bCOMPND)
786 {
787 strcat(title,"; ");
788 strcat(title,c);
789 }
790 else
791 {
792 strcpy(title,c);
793 }
794 }
795 bCOMPND=TRUE;
796 }
797 break;
798
799 case epdbTER:
800 chainnum++;
801 break;
802
803 case epdbMODEL:
804 if(model_nr)
805 {
806 sscanf(line,"%*s%d",model_nr);
807 }
808 break;
809
810 case epdbENDMDL:
811 bStop=TRUE;
812 break;
813 case epdbCONECT:
814 if (gc)
815 {
816 gmx_conect_addline(gc,line);
817 }
818 else if (!bConnWarn)
819 {
820 fprintf(stderr,"WARNING: all CONECT records are ignored\n");
821 bConnWarn = TRUE;
822 }
823 break;
824
825 default:
826 break;
827 }
828 }
829
830 free_symtab(&symtab);
831 return natom;
832 }
833
834 void get_pdb_coordnum(FILE *in,int *natoms)
835 {
836 char line[STRLEN];
837
838 *natoms=0;
839 while (fgets2(line,STRLEN,in))
840 {
841 if ( strncmp(line,"ENDMDL",6) == 0 )
842 {
843 break;
844 }
845 if ((strncmp(line,"ATOM ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
846 {
847 (*natoms)++;
848 }
849 }
850 }
851
852 void read_pdb_conf(const char *infile,char *title,
853 t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
854 gmx_conect conect)
855 {
856 FILE *in;
857
858 in = gmx_fio_fopen(infile,"r");
859 read_pdbfile(in,title,NULL,atoms,x,ePBC,box,bChange,conect);
860 gmx_fio_fclose(in);
861 }
862
863 gmx_conect gmx_conect_generate(t_topology *top)
864 {
865 int f,i;
866 gmx_conect gc;
867
868 /* Fill the conect records */
869 gc = gmx_conect_init();
870
871 for(f=0; (f<F_NRE); f++) {
872 if (IS_CHEMBOND(f))
873 for(i=0; (i<top->idef.il[f].nr); i+=interaction_function[f].nratoms+1) {
874 gmx_conect_add(gc,top->idef.il[f].iatoms[i+1],
875 top->idef.il[f].iatoms[i+2]);
876 }
877 }
878 return gc;
879 }