2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
63 typedef struct gmx_conect_t {
66 gmx_conection_t *conect;
69 static const char *pdbtp[epdbNR]={
70 "ATOM ","HETATM", "ANISOU", "CRYST1",
71 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
76 /* this is not very good,
77 but these are only used in gmx_trjconv and gmx_editconv */
78 static gmx_bool bWideFormat=FALSE;
79 #define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
81 void set_pdb_wide_format(gmx_bool bSet)
86 static void xlate_atomname_pdb2gmx(char *name)
92 if (length>3 && isdigit(name[0])) {
94 for(i=1; i<length; i++)
100 static void xlate_atomname_gmx2pdb(char *name)
106 if (length>3 && isdigit(name[length-1])) {
108 for(i=length-1; i>0; --i)
115 void gmx_write_pdb_box(FILE *out,int ePBC,matrix box)
117 real alpha,beta,gamma;
120 ePBC = guess_ePBC(box);
122 if (ePBC == epbcNONE)
125 if (norm2(box[YY])*norm2(box[ZZ])!=0)
126 alpha = RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ]));
129 if (norm2(box[XX])*norm2(box[ZZ])!=0)
130 beta = RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ]));
133 if (norm2(box[XX])*norm2(box[YY])!=0)
134 gamma = RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY]));
137 fprintf(out,"REMARK THIS IS A SIMULATION BOX\n");
138 if (ePBC != epbcSCREW) {
139 fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
140 10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
141 alpha,beta,gamma,"P 1",1);
143 /* Double the a-vector length and write the correct space group */
144 fprintf(out,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
145 20*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
146 alpha,beta,gamma,"P 21 1 1",1);
151 static void read_cryst1(char *line,int *ePBC,matrix box)
154 char sa[12],sb[12],sc[12],sg[SG_SIZE+1],ident;
155 double fa,fb,fc,alpha,beta,gamma,cosa,cosb,cosg,sing;
159 sscanf(line,"%*s%s%s%s%lf%lf%lf",sa,sb,sc,&alpha,&beta,&gamma);
162 if (strlen(line) >= 55) {
163 strncpy(sg,line+55,SG_SIZE);
169 sscanf(sg,"%c %d %d %d",&ident,&syma,&symb,&symc);
170 if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1) {
171 fc = strtod(sc,NULL)*0.1;
172 ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
174 if (ident == 'P' && syma == 21 && symb == 1 && symc == 1) {
175 ePBC_file = epbcSCREW;
183 fa = strtod(sa,NULL)*0.1;
184 fb = strtod(sb,NULL)*0.1;
185 fc = strtod(sc,NULL)*0.1;
186 if (ePBC_file == epbcSCREW) {
191 if ((alpha!=90.0) || (beta!=90.0) || (gamma!=90.0)) {
193 cosa = cos(alpha*DEG2RAD);
198 cosb = cos(beta*DEG2RAD);
203 cosg = cos(gamma*DEG2RAD);
204 sing = sin(gamma*DEG2RAD);
209 box[YY][XX] = fb*cosg;
210 box[YY][YY] = fb*sing;
211 box[ZZ][XX] = fc*cosb;
212 box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
213 box[ZZ][ZZ] = sqrt(fc*fc
214 - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
222 void write_pdbfile_indexed(FILE *out,const char *title,
223 t_atoms *atoms,rvec x[],
224 int ePBC,matrix box,char chainid,
225 int model_nr, atom_id nindex, atom_id index[],
226 gmx_conect conect, gmx_bool bTerSepChains)
228 gmx_conect_t *gc = (gmx_conect_t *)conect;
229 char resnm[6],nm[6],pdbform[128],pukestring[100];
231 int resind,resnr,type;
232 unsigned char resic,ch;
236 int chainnum,lastchainnum;
237 int lastresind,lastchainresind;
238 gmx_residuetype_t rt;
239 const char *p_restype;
240 const char *p_lastrestype;
242 gmx_residuetype_init(&rt);
244 bromacs(pukestring,99);
245 fprintf(out,"TITLE %s\n",(title && title[0])?title:pukestring);
247 fprintf(out,"REMARK This file does not adhere to the PDB standard\n");
248 fprintf(out,"REMARK As a result of, some programs may not like it\n");
250 if (box && ( norm2(box[XX]) || norm2(box[YY]) || norm2(box[ZZ]) ) ) {
251 gmx_write_pdb_box(out,ePBC,box);
253 if (atoms->pdbinfo) {
254 /* Check whether any occupancies are set, in that case leave it as is,
255 * otherwise set them all to one
258 for (ii=0; (ii<nindex) && bOccup; ii++) {
260 bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
266 fprintf(out,"MODEL %8d\n",model_nr>0 ? model_nr : 1);
268 lastchainresind = -1;
273 for (ii=0; ii<nindex; ii++) {
276 resind = atoms->atom[i].resind;
277 chainnum = atoms->resinfo[resind].chainnum;
278 p_lastrestype = p_restype;
279 gmx_residuetype_get_type(rt,*atoms->resinfo[resind].name,&p_restype);
281 /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
282 if( bTerSepChains && ii>0 && chainnum != lastchainnum)
284 /* Only add TER if the previous chain contained protein/DNA/RNA. */
285 if(gmx_residuetype_is_protein(rt,p_lastrestype) || gmx_residuetype_is_dna(rt,p_lastrestype) || gmx_residuetype_is_rna(rt,p_lastrestype))
287 fprintf(out,"TER\n");
289 lastchainnum = chainnum;
290 lastchainresind = lastresind;
293 strncpy(resnm,*atoms->resinfo[resind].name,sizeof(resnm)-1);
294 strncpy(nm,*atoms->atomname[i],sizeof(nm)-1);
295 /* rename HG12 to 2HG1, etc. */
296 xlate_atomname_gmx2pdb(nm);
297 resnr = atoms->resinfo[resind].nr;
298 resic = atoms->resinfo[resind].ic;
305 ch = atoms->resinfo[resind].chainid;
313 resnr = resnr % 10000;
314 if (atoms->pdbinfo) {
315 type = atoms->pdbinfo[i].type;
316 occup = bOccup ? 1.0 : atoms->pdbinfo[i].occup;
317 bfac = atoms->pdbinfo[i].bfac;
326 "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
328 /* Check whether atomname is an element name */
329 if ((strlen(nm)<4) && (gmx_strcasecmp(nm,atoms->atom[i].elem) != 0))
330 strcpy(pdbform,pdbformat);
332 strcpy(pdbform,pdbformat4);
333 if (strlen(nm) > 4) {
335 if (nlongname < maxwln) {
336 fprintf(stderr,"WARNING: Writing out atom name (%s) longer than 4 characters to .pdb file\n",nm);
337 } else if (nlongname == maxwln) {
338 fprintf(stderr,"WARNING: More than %d long atom names, will not write more warnings\n",maxwln);
343 strcat(pdbform,"%6.2f%6.2f %2s\n");
345 fprintf(out,pdbform,pdbtp[type],(i+1)%100000,nm,resnm,ch,resnr,
346 (resic == '\0') ? ' ' : resic,
347 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],occup,bfac,atoms->atom[i].elem);
348 if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic) {
349 fprintf(out,"ANISOU%5u %-4.4s%3.3s %c%4d%c %7d%7d%7d%7d%7d%7d\n",
350 (i+1)%100000,nm,resnm,ch,resnr,
351 (resic == '\0') ? ' ' : resic,
352 atoms->pdbinfo[i].uij[0],atoms->pdbinfo[i].uij[1],
353 atoms->pdbinfo[i].uij[2],atoms->pdbinfo[i].uij[3],
354 atoms->pdbinfo[i].uij[4],atoms->pdbinfo[i].uij[5]);
358 fprintf(out,"TER\n");
359 fprintf(out,"ENDMDL\n");
362 /* Write conect records */
363 for(i=0; (i<gc->nconect); i++) {
364 fprintf(out,"CONECT%5d%5d\n",gc->conect[i].ai+1,gc->conect[i].aj+1);
368 gmx_residuetype_destroy(rt);
371 void write_pdbfile(FILE *out,const char *title, t_atoms *atoms,rvec x[],
372 int ePBC,matrix box,char chainid,int model_nr,gmx_conect conect,gmx_bool bTerSepChains)
376 snew(index,atoms->nr);
377 for(i=0; i<atoms->nr; i++)
379 write_pdbfile_indexed(out,title,atoms,x,ePBC,box,chainid,model_nr,
380 atoms->nr,index,conect,bTerSepChains);
384 int line2type(char *line)
393 for(k=0; (k<epdbNR); k++)
394 if (strncmp(type,pdbtp[k],strlen(pdbtp[k])) == 0)
400 static void read_anisou(char line[],int natom,t_atoms *atoms)
404 char anr[12],anm[12];
408 for(k=0; (k<5); k++,j++) anr[k]=line[j];
411 for(k=0; (k<4); k++,j++) anm[k]=line[j];
415 /* Strip off spaces */
418 /* Search backwards for number and name only */
419 atomnr = strtol(anr, NULL, 10);
420 for(i=natom-1; (i>=0); i--)
421 if ((strcmp(anm,*(atoms->atomname[i])) == 0) &&
422 (atomnr == atoms->pdbinfo[i].atomnr))
425 fprintf(stderr,"Skipping ANISOU record (atom %s %d not found)\n",
428 if (sscanf(line+29,"%d%d%d%d%d%d",
429 &atoms->pdbinfo[i].uij[U11],&atoms->pdbinfo[i].uij[U22],
430 &atoms->pdbinfo[i].uij[U33],&atoms->pdbinfo[i].uij[U12],
431 &atoms->pdbinfo[i].uij[U13],&atoms->pdbinfo[i].uij[U23])
433 atoms->pdbinfo[i].bAnisotropic = TRUE;
436 fprintf(stderr,"Invalid ANISOU record for atom %d\n",i);
437 atoms->pdbinfo[i].bAnisotropic = FALSE;
442 void get_pdb_atomnumber(t_atoms *atoms,gmx_atomprop_t aps)
444 int i,atomnumber,len;
446 char anm[6],anm_copy[6],*ptr;
450 if (!atoms->pdbinfo) {
451 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
453 for(i=0; (i<atoms->nr); i++) {
454 strcpy(anm,atoms->pdbinfo[i].atomnm);
455 strcpy(anm_copy,atoms->pdbinfo[i].atomnm);
458 if ((anm[0] != ' ') && ((len <=2) || ((len > 2) && !isdigit(anm[2])))) {
460 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
461 atomnumber = gmx_nint(eval);
464 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
465 atomnumber = gmx_nint(eval);
468 if (atomnumber == NOTSET) {
470 while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
472 anm_copy[0] = anm[k];
474 if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
475 atomnumber = gmx_nint(eval);
477 atoms->atom[i].atomnumber = atomnumber;
478 ptr = gmx_atomprop_element(aps,atomnumber);
479 strncpy(atoms->atom[i].elem,ptr==NULL ? "" : ptr,4);
481 fprintf(debug,"Atomnumber for atom '%s' is %d\n",anm,atomnumber);
485 static int read_atom(t_symtab *symtab,
486 char line[],int type,int natom,
487 t_atoms *atoms,rvec x[],int chainnum,gmx_bool bChange)
492 char anr[12],anm[12],anm_copy[12],altloc,resnm[12],rnr[12];
493 char xc[12],yc[12],zc[12],occup[12],bfac[12];
496 int resnr,atomnumber;
498 if (natom>=atoms->nr)
499 gmx_fatal(FARGS,"\nFound more atoms (%d) in pdb file than expected (%d)",
504 for(k=0; (k<5); k++,j++) anr[k]=line[j];
508 for(k=0; (k<4); k++,j++) anm[k]=line[j];
510 strcpy(anm_copy,anm);
516 for(k=0; (k<4); k++,j++)
524 for(k=0; (k<4); k++,j++) {
529 resnr = strtol(rnr, NULL, 10);
533 /* X,Y,Z Coordinate */
534 for(k=0; (k<8); k++,j++) xc[k]=line[j];
536 for(k=0; (k<8); k++,j++) yc[k]=line[j];
538 for(k=0; (k<8); k++,j++) zc[k]=line[j];
542 for(k=0; (k<6); k++,j++) occup[k]=line[j];
546 for(k=0; (k<7); k++,j++) bfac[k]=line[j];
550 atomn=&(atoms->atom[natom]);
552 atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
553 atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
554 (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name,resnm) != 0))
559 atomn->resind = atoms->atom[natom-1].resind + 1;
561 atoms->nres = atomn->resind + 1;
562 t_atoms_set_resinfo(atoms,natom,symtab,resnm,resnr,resic,chainnum,chainid);
566 atomn->resind = atoms->atom[natom-1].resind;
569 xlate_atomname_pdb2gmx(anm);
571 atoms->atomname[natom]=put_symtab(symtab,anm);
574 atomn->atomnumber = atomnumber;
575 atomn->elem[0] = '\0';
577 x[natom][XX]=strtod(xc,NULL)*0.1;
578 x[natom][YY]=strtod(yc,NULL)*0.1;
579 x[natom][ZZ]=strtod(zc,NULL)*0.1;
580 if (atoms->pdbinfo) {
581 atoms->pdbinfo[natom].type=type;
582 atoms->pdbinfo[natom].atomnr=strtol(anr, NULL, 10);
583 atoms->pdbinfo[natom].altloc=altloc;
584 strcpy(atoms->pdbinfo[natom].atomnm,anm_copy);
585 atoms->pdbinfo[natom].bfac=strtod(bfac,NULL);
586 atoms->pdbinfo[natom].occup=strtod(occup,NULL);
593 gmx_bool is_hydrogen(const char *nm)
602 else if ((isdigit(buf[0])) && (buf[1] == 'H'))
607 gmx_bool is_dummymass(const char *nm)
614 if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
620 static void gmx_conect_addline(gmx_conect_t *con,char *line)
623 char format[32],form2[32];
625 sprintf(form2,"%%*s");
626 sprintf(format,"%s%%d",form2);
627 if (sscanf(line,format,&ai) == 1) {
630 sprintf(format,"%s%%d",form2);
631 n = sscanf(line,format,&aj);
633 srenew(con->conect,++con->nconect);
634 con->conect[con->nconect-1].ai = ai-1;
635 con->conect[con->nconect-1].aj = aj-1;
641 void gmx_conect_dump(FILE *fp,gmx_conect conect)
643 gmx_conect_t *gc = (gmx_conect_t *)conect;
646 for(i=0; (i<gc->nconect); i++)
647 fprintf(fp,"%6s%5d%5d\n","CONECT",
648 gc->conect[i].ai+1,gc->conect[i].aj+1);
651 gmx_conect gmx_conect_init()
657 return (gmx_conect) gc;
660 void gmx_conect_done(gmx_conect conect)
662 gmx_conect_t *gc = (gmx_conect_t *)conect;
667 gmx_bool gmx_conect_exist(gmx_conect conect,int ai,int aj)
669 gmx_conect_t *gc = (gmx_conect_t *)conect;
675 for(i=0; (i<gc->nconect); i++)
676 if (((gc->conect[i].ai == ai) &&
677 (gc->conect[i].aj == aj)) ||
678 ((gc->conect[i].aj == ai) &&
679 (gc->conect[i].ai == aj)))
684 void gmx_conect_add(gmx_conect conect,int ai,int aj)
686 gmx_conect_t *gc = (gmx_conect_t *)conect;
692 if (!gmx_conect_exist(conect,ai,aj)) {
693 srenew(gc->conect,++gc->nconect);
694 gc->conect[gc->nconect-1].ai = ai;
695 gc->conect[gc->nconect-1].aj = aj;
699 int read_pdbfile(FILE *in,char *title,int *model_nr,
700 t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
703 gmx_conect_t *gc = (gmx_conect_t *)conect;
706 gmx_bool bConnWarn = FALSE;
710 int natom,chainnum,nres_ter_prev=0;
712 gmx_bool bStop=FALSE;
716 /* Only assume pbc when there is a CRYST1 entry */
724 open_symtab(&symtab);
730 while (!bStop && (fgets2(line,STRLEN,in) != NULL))
732 line_type = line2type(line);
738 natom = read_atom(&symtab,line,line_type,natom,atoms,x,chainnum,bChange);
744 read_anisou(line,natom,atoms);
749 read_cryst1(line,ePBC,box);
754 if (strlen(line) > 6)
757 /* skip HEADER or TITLE and spaces */
758 while (c[0]!=' ') c++;
759 while (c[0]==' ') c++;
760 /* truncate after title */
774 if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:")))
776 if ( !(c=strstr(line+6,"MOLECULE:")) )
780 /* skip 'MOLECULE:' and spaces */
781 while (c[0]!=' ') c++;
782 while (c[0]==' ') c++;
783 /* truncate after title */
787 while ( (d[-1]==';') && d>c) d--;
813 sscanf(line,"%*s%d",model_nr);
823 gmx_conect_addline(gc,line);
827 fprintf(stderr,"WARNING: all CONECT records are ignored\n");
837 free_symtab(&symtab);
841 void get_pdb_coordnum(FILE *in,int *natoms)
846 while (fgets2(line,STRLEN,in))
848 if ( strncmp(line,"ENDMDL",6) == 0 )
852 if ((strncmp(line,"ATOM ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
859 void read_pdb_conf(const char *infile,char *title,
860 t_atoms *atoms,rvec x[],int *ePBC,matrix box,gmx_bool bChange,
865 in = gmx_fio_fopen(infile,"r");
866 read_pdbfile(in,title,NULL,atoms,x,ePBC,box,bChange,conect);
870 gmx_conect gmx_conect_generate(t_topology *top)
875 /* Fill the conect records */
876 gc = gmx_conect_init();
878 for(f=0; (f<F_NRE); f++) {
880 for(i=0; (i<top->idef.il[f].nr); i+=interaction_function[f].nratoms+1) {
881 gmx_conect_add(gc,top->idef.il[f].iatoms[i+1],
882 top->idef.il[f].iatoms[i+2]);