3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
60 typedef struct gmx_conect_t
{
63 gmx_conection_t
*conect
;
66 static const char *pdbtp
[epdbNR
]={
67 "ATOM ","HETATM", "ANISOU", "CRYST1",
68 "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
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"
78 void set_pdb_wide_format(gmx_bool bSet
)
83 static void xlate_atomname_pdb2gmx(char *name
)
89 if (length
>3 && isdigit(name
[0])) {
91 for(i
=1; i
<length
; i
++)
97 static void xlate_atomname_gmx2pdb(char *name
)
103 if (length
>3 && isdigit(name
[length
-1])) {
105 for(i
=length
-1; i
>0; --i
)
112 void gmx_write_pdb_box(FILE *out
,int ePBC
,matrix box
)
114 real alpha
,beta
,gamma
;
117 ePBC
= guess_ePBC(box
);
119 if (ePBC
== epbcNONE
)
122 if (norm2(box
[YY
])*norm2(box
[ZZ
])!=0)
123 alpha
= RAD2DEG
*acos(cos_angle_no_table(box
[YY
],box
[ZZ
]));
126 if (norm2(box
[XX
])*norm2(box
[ZZ
])!=0)
127 beta
= RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[ZZ
]));
130 if (norm2(box
[XX
])*norm2(box
[YY
])!=0)
131 gamma
= RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[YY
]));
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);
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);
148 static void read_cryst1(char *line
,int *ePBC
,matrix box
)
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
;
156 sscanf(line
,"%*s%s%s%s%lf%lf%lf",sa
,sb
,sc
,&alpha
,&beta
,&gamma
);
159 if (strlen(line
) >= 55) {
160 strncpy(sg
,line
+55,SG_SIZE
);
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
);
171 if (ident
== 'P' && syma
== 21 && symb
== 1 && symc
== 1) {
172 ePBC_file
= epbcSCREW
;
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
) {
188 if ((alpha
!=90.0) || (beta
!=90.0) || (gamma
!=90.0)) {
190 cosa
= cos(alpha
*DEG2RAD
);
195 cosb
= cos(beta
*DEG2RAD
);
200 cosg
= cos(gamma
*DEG2RAD
);
201 sing
= sin(gamma
*DEG2RAD
);
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
]);
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
)
225 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
226 char resnm
[6],nm
[6],pdbform
[128],pukestring
[100];
228 int resind
,resnr
,type
;
229 unsigned char resic
,ch
;
233 int chainnum
,lastchainnum
;
234 int lastresind
,lastchainresind
;
235 gmx_residuetype_t rt
;
236 const char *p_restype
;
237 const char *p_lastrestype
;
239 gmx_residuetype_init(&rt
);
241 bromacs(pukestring
,99);
242 fprintf(out
,"TITLE %s\n",(title
&& title
[0])?title
:pukestring
);
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");
247 if (box
&& ( norm2(box
[XX
]) || norm2(box
[YY
]) || norm2(box
[ZZ
]) ) ) {
248 gmx_write_pdb_box(out
,ePBC
,box
);
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
255 for (ii
=0; (ii
<nindex
) && bOccup
; ii
++) {
257 bOccup
= bOccup
&& (atoms
->pdbinfo
[i
].occup
== 0.0);
263 fprintf(out
,"MODEL %8d\n",model_nr
>=0 ? model_nr
: 1);
265 lastchainresind
= -1;
270 for (ii
=0; ii
<nindex
; ii
++) {
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
);
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
)
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
))
284 fprintf(out
,"TER\n");
286 lastchainnum
= chainnum
;
287 lastchainresind
= lastresind
;
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
;
302 ch
= atoms
->resinfo
[resind
].chainid
;
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
;
323 "%-6s%5u %-4.4s %3.3s %c%4d%c %10.5f%10.5f%10.5f%8.4f%8.4f %2s\n");
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
);
329 strcpy(pdbform
,pdbformat4
);
330 if (strlen(nm
) > 4) {
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
);
340 strcat(pdbform
,"%6.2f%6.2f %2s\n");
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]);
355 fprintf(out
,"TER\n");
356 fprintf(out
,"ENDMDL\n");
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);
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
)
371 snew(index
,atoms
->nr
);
372 for(i
=0; i
<atoms
->nr
; i
++)
374 write_pdbfile_indexed(out
,title
,atoms
,x
,ePBC
,box
,chainid
,model_nr
,
375 atoms
->nr
,index
,conect
,bTerSepChains
);
379 int line2type(char *line
)
388 for(k
=0; (k
<epdbNR
); k
++)
389 if (strncmp(type
,pdbtp
[k
],strlen(pdbtp
[k
])) == 0)
395 static void read_anisou(char line
[],int natom
,t_atoms
*atoms
)
399 char anr
[12],anm
[12];
403 for(k
=0; (k
<5); k
++,j
++) anr
[k
]=line
[j
];
406 for(k
=0; (k
<4); k
++,j
++) anm
[k
]=line
[j
];
410 /* Strip off spaces */
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
))
420 fprintf(stderr
,"Skipping ANISOU record (atom %s %d not found)\n",
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
])
428 atoms
->pdbinfo
[i
].bAnisotropic
= TRUE
;
431 fprintf(stderr
,"Invalid ANISOU record for atom %d\n",i
);
432 atoms
->pdbinfo
[i
].bAnisotropic
= FALSE
;
437 void get_pdb_atomnumber(t_atoms
*atoms
,gmx_atomprop_t aps
)
439 int i
,atomnumber
,len
;
441 char anm
[6],anm_copy
[6];
445 if (!atoms
->pdbinfo
) {
446 gmx_incons("Trying to deduce atomnumbers when no pdb information is present");
448 for(i
=0; (i
<atoms
->nr
); i
++) {
449 strcpy(anm
,atoms
->pdbinfo
[i
].atomnm
);
450 strcpy(anm_copy
,atoms
->pdbinfo
[i
].atomnm
);
453 if ((anm
[0] != ' ') && ((len
<=2) || ((len
> 2) && !isdigit(anm
[2])))) {
455 if (gmx_atomprop_query(aps
,epropElement
,"???",anm_copy
,&eval
))
456 atomnumber
= gmx_nint(eval
);
459 if (gmx_atomprop_query(aps
,epropElement
,"???",anm_copy
,&eval
))
460 atomnumber
= gmx_nint(eval
);
463 if (atomnumber
== NOTSET
) {
465 while ((k
< strlen(anm
)) && (isspace(anm
[k
]) || isdigit(anm
[k
])))
467 anm_copy
[0] = anm
[k
];
469 if (gmx_atomprop_query(aps
,epropElement
,"???",anm_copy
,&eval
))
470 atomnumber
= gmx_nint(eval
);
472 atoms
->atom
[i
].atomnumber
= atomnumber
;
473 sprintf(atoms
->atom
[i
].elem
,"%2s",gmx_atomprop_element(aps
,atomnumber
));
475 fprintf(debug
,"Atomnumber for atom '%s' is %d\n",anm
,atomnumber
);
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
)
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];
490 int resnr
,atomnumber
;
492 if (natom
>=atoms
->nr
)
493 gmx_fatal(FARGS
,"\nFound more atoms (%d) in pdb file than expected (%d)",
498 for(k
=0; (k
<5); k
++,j
++) anr
[k
]=line
[j
];
502 for(k
=0; (k
<4); k
++,j
++) anm
[k
]=line
[j
];
504 strcpy(anm_copy
,anm
);
509 for(k
=0; (k
<4); k
++,j
++)
517 for(k
=0; (k
<4); k
++,j
++) {
522 resnr
= strtol(rnr
, NULL
, 10);
526 /* X,Y,Z Coordinate */
527 for(k
=0; (k
<8); k
++,j
++) xc
[k
]=line
[j
];
529 for(k
=0; (k
<8); k
++,j
++) yc
[k
]=line
[j
];
531 for(k
=0; (k
<8); k
++,j
++) zc
[k
]=line
[j
];
535 for(k
=0; (k
<6); k
++,j
++) occup
[k
]=line
[j
];
539 for(k
=0; (k
<7); k
++,j
++) bfac
[k
]=line
[j
];
543 atomn
=&(atoms
->atom
[natom
]);
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))
552 atomn
->resind
= atoms
->atom
[natom
-1].resind
+ 1;
554 atoms
->nres
= atomn
->resind
+ 1;
555 t_atoms_set_resinfo(atoms
,natom
,symtab
,resnm
,resnr
,resic
,chainnum
,chainid
);
559 atomn
->resind
= atoms
->atom
[natom
-1].resind
;
562 xlate_atomname_pdb2gmx(anm
);
564 atoms
->atomname
[natom
]=put_symtab(symtab
,anm
);
567 atomn
->atomnumber
= atomnumber
;
568 atomn
->elem
[0] = '\0';
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
);
586 gmx_bool
is_hydrogen(const char *nm
)
595 else if ((isdigit(buf
[0])) && (buf
[1] == 'H'))
600 gmx_bool
is_dummymass(const char *nm
)
607 if ((buf
[0] == 'M') && isdigit(buf
[strlen(buf
)-1]))
613 static void gmx_conect_addline(gmx_conect_t
*con
,char *line
)
616 char format
[32],form2
[32];
618 sprintf(form2
,"%%*s");
619 sprintf(format
,"%s%%d",form2
);
620 if (sscanf(line
,format
,&ai
) == 1) {
623 sprintf(format
,"%s%%d",form2
);
624 n
= sscanf(line
,format
,&aj
);
626 srenew(con
->conect
,++con
->nconect
);
627 con
->conect
[con
->nconect
-1].ai
= ai
-1;
628 con
->conect
[con
->nconect
-1].aj
= aj
-1;
634 void gmx_conect_dump(FILE *fp
,gmx_conect conect
)
636 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
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);
644 gmx_conect
gmx_conect_init()
650 return (gmx_conect
) gc
;
653 void gmx_conect_done(gmx_conect conect
)
655 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
660 gmx_bool
gmx_conect_exist(gmx_conect conect
,int ai
,int aj
)
662 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
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
)))
677 void gmx_conect_add(gmx_conect conect
,int ai
,int aj
)
679 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
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
;
692 int read_pdbfile(FILE *in
,char *title
,int *model_nr
,
693 t_atoms
*atoms
,rvec x
[],int *ePBC
,matrix box
,gmx_bool bChange
,
696 gmx_conect_t
*gc
= (gmx_conect_t
*)conect
;
699 gmx_bool bConnWarn
= FALSE
;
703 int natom
,chainnum
,nres_ter_prev
=0;
705 gmx_bool bStop
=FALSE
;
709 /* Only assume pbc when there is a CRYST1 entry */
717 open_symtab(&symtab
);
723 while (!bStop
&& (fgets2(line
,STRLEN
,in
) != NULL
))
725 line_type
= line2type(line
);
731 natom
= read_atom(&symtab
,line
,line_type
,natom
,atoms
,x
,chainnum
,bChange
);
737 read_anisou(line
,natom
,atoms
);
742 read_cryst1(line
,ePBC
,box
);
747 if (strlen(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 */
767 if ((!strstr(line
,": ")) || (strstr(line
+6,"MOLECULE:")))
769 if ( !(c
=strstr(line
+6,"MOLECULE:")) )
773 /* skip 'MOLECULE:' and spaces */
774 while (c
&& (c
[0]!=' ')) c
++;
775 while (c
&& (c
[0]==' ')) c
++;
776 /* truncate after title */
780 while ( (d
[-1]==';') && d
>c
) d
--;
806 sscanf(line
,"%*s%d",model_nr
);
816 gmx_conect_addline(gc
,line
);
820 fprintf(stderr
,"WARNING: all CONECT records are ignored\n");
830 free_symtab(&symtab
);
834 void get_pdb_coordnum(FILE *in
,int *natoms
)
839 while (fgets2(line
,STRLEN
,in
))
841 if ( strncmp(line
,"ENDMDL",6) == 0 )
845 if ((strncmp(line
,"ATOM ",6) == 0) || (strncmp(line
,"HETATM",6) == 0))
852 void read_pdb_conf(const char *infile
,char *title
,
853 t_atoms
*atoms
,rvec x
[],int *ePBC
,matrix box
,gmx_bool bChange
,
858 in
= gmx_fio_fopen(infile
,"r");
859 read_pdbfile(in
,title
,NULL
,atoms
,x
,ePBC
,box
,bChange
,conect
);
863 gmx_conect
gmx_conect_generate(t_topology
*top
)
868 /* Fill the conect records */
869 gc
= gmx_conect_init();
871 for(f
=0; (f
<F_NRE
); 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]);