Project

General

Profile

pdbio.c

fixed pdbio.c - Kyle Beauchamp, 08/22/2010 10:22 PM

 
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
  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 bool bWideFormat=FALSE;
76
#define REMARK_SIM_BOX "REMARK    THIS IS A SIMULATION BOX"
77

    
78
void set_pdb_wide_format(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, 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
  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,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;
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
    atomnumber = NOTSET;
452
    if (anm[0] != ' ') {
453
      anm_copy[2] = nc;
454
      if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
455
        atomnumber = gmx_nint(eval);
456
      else {
457
        anm_copy[1] = nc;
458
        if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
459
          atomnumber = gmx_nint(eval);
460
      }
461
    }
462
    if (atomnumber == NOTSET) {
463
      k=0;
464
      while ((k < strlen(anm)) && (isspace(anm[k]) || isdigit(anm[k])))
465
        k++;
466
      anm_copy[0] = anm[k];
467
      anm_copy[1] = nc;
468
      if (gmx_atomprop_query(aps,epropElement,"???",anm_copy,&eval))
469
        atomnumber = gmx_nint(eval);
470
    }
471
    atoms->atom[i].atomnumber = atomnumber;
472
    sprintf(atoms->atom[i].elem,"%2s",gmx_atomprop_element(aps,atomnumber));
473
    if (debug)
474
      fprintf(debug,"Atomnumber for atom '%s' is %d\n",anm,atomnumber);
475
  }
476
}
477

    
478
static int read_atom(t_symtab *symtab,
479
                     char line[],int type,int natom,
480
                     t_atoms *atoms,rvec x[],int chainnum,bool bChange)
481
{
482
  t_atom *atomn;
483
  int  j,k;
484
  char nc='\0';
485
  char anr[12],anm[12],anm_copy[12],altloc,resnm[12],rnr[12];
486
  char xc[12],yc[12],zc[12],occup[12],bfac[12];
487
  unsigned char resic;
488
  char chainid;
489
  int  resnr,atomnumber;
490

    
491
  if (natom>=atoms->nr)
492
    gmx_fatal(FARGS,"\nFound more atoms (%d) in pdb file than expected (%d)",
493
              natom+1,atoms->nr);
494

    
495
  /* Skip over type */  
496
  j=6;
497
  for(k=0; (k<5); k++,j++) anr[k]=line[j];
498
  anr[k]=nc;
499
  trim(anr);
500
  j++;
501
  for(k=0; (k<4); k++,j++) anm[k]=line[j];
502
  anm[k]=nc;
503
  strcpy(anm_copy,anm);
504
  atomnumber = NOTSET;
505
  trim(anm);
506
  altloc=line[j];
507
  j++;
508
  for(k=0; (k<4); k++,j++) 
509
    resnm[k]=line[j];
510
  resnm[k]=nc;
511
  trim(resnm);
512

    
513
  chainid = line[j];
514
  j++;
515
  
516
  for(k=0; (k<4); k++,j++) {
517
    rnr[k] = line[j];
518
  }
519
  rnr[k] = nc;
520
  trim(rnr);
521
  resnr = strtol(rnr, NULL, 10); 
522
  resic = line[j];
523
  j+=4;
524

    
525
  /* X,Y,Z Coordinate */
526
  for(k=0; (k<8); k++,j++) xc[k]=line[j];
527
  xc[k]=nc;
528
  for(k=0; (k<8); k++,j++) yc[k]=line[j];
529
  yc[k]=nc;
530
  for(k=0; (k<8); k++,j++) zc[k]=line[j];
531
  zc[k]=nc;
532
  
533
  /* Weight */
534
  for(k=0; (k<6); k++,j++) occup[k]=line[j];
535
  occup[k]=nc;
536
  
537
  /* B-Factor */
538
  for(k=0; (k<7); k++,j++) bfac[k]=line[j];
539
  bfac[k]=nc;
540

    
541
  if (atoms->atom) {
542
    atomn=&(atoms->atom[natom]);
543
    if ((natom==0) ||
544
        atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
545
        atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
546
        (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name,resnm) != 0))
547
    {
548
      if (natom == 0) {
549
        atomn->resind = 0;
550
      } else {
551
        atomn->resind = atoms->atom[natom-1].resind + 1;
552
      }
553
      atoms->nres = atomn->resind + 1;
554
      t_atoms_set_resinfo(atoms,natom,symtab,resnm,resnr,resic,chainnum,chainid);
555
    }
556
    else
557
    {
558
      atomn->resind = atoms->atom[natom-1].resind;
559
    }
560
    if (bChange) {
561
      xlate_atomname_pdb2gmx(anm); 
562
    }
563
    atoms->atomname[natom]=put_symtab(symtab,anm);
564
    atomn->m = 0.0;
565
    atomn->q = 0.0;
566
    atomn->atomnumber = atomnumber;
567
    atomn->elem[0] = '\0';
568
  }
569
  x[natom][XX]=strtod(xc,NULL)*0.1;
570
  x[natom][YY]=strtod(yc,NULL)*0.1;
571
  x[natom][ZZ]=strtod(zc,NULL)*0.1;
572
  if (atoms->pdbinfo) {
573
    atoms->pdbinfo[natom].type=type;
574
    atoms->pdbinfo[natom].atomnr=strtol(anr, NULL, 10); 
575
    atoms->pdbinfo[natom].altloc=altloc;
576
    strcpy(atoms->pdbinfo[natom].atomnm,anm_copy);
577
    atoms->pdbinfo[natom].bfac=strtod(bfac,NULL);
578
    atoms->pdbinfo[natom].occup=strtod(occup,NULL);
579
  }
580
  natom++;
581
  
582
  return natom;
583
}
584

    
585
bool is_hydrogen(const char *nm)
586
{
587
  char buf[30];
588
  
589
  strcpy(buf,nm);
590
  trim(buf);
591
  
592
  if (buf[0] == 'H')
593
    return TRUE;
594
  else if ((isdigit(buf[0])) && (buf[1] == 'H'))
595
    return TRUE;
596
  return FALSE;
597
}
598

    
599
bool is_dummymass(const char *nm)
600
{
601
  char buf[30];
602
  
603
  strcpy(buf,nm);
604
  trim(buf);
605
  
606
  if ((buf[0] == 'M') && isdigit(buf[strlen(buf)-1]))
607
    return TRUE;
608
      
609
  return FALSE;
610
}
611

    
612
static void gmx_conect_addline(gmx_conect_t *con,char *line)
613
{
614
  int n,ai,aj;
615
  char format[32],form2[32];
616
  
617
  sprintf(form2,"%%*s");
618
  sprintf(format,"%s%%d",form2);
619
  if (sscanf(line,format,&ai) == 1) {
620
    do {
621
      strcat(form2,"%*s");
622
      sprintf(format,"%s%%d",form2);
623
      n = sscanf(line,format,&aj);
624
      if (n == 1) {
625
        srenew(con->conect,++con->nconect);
626
        con->conect[con->nconect-1].ai = ai-1;
627
        con->conect[con->nconect-1].aj = aj-1;
628
      }
629
    } while (n == 1);
630
  }
631
}
632

    
633
void gmx_conect_dump(FILE *fp,gmx_conect conect)
634
{
635
  gmx_conect_t *gc = (gmx_conect_t *)conect;
636
  int i;
637
  
638
  for(i=0; (i<gc->nconect); i++)
639
    fprintf(fp,"%6s%5d%5d\n","CONECT",
640
            gc->conect[i].ai+1,gc->conect[i].aj+1);
641
}
642

    
643
gmx_conect gmx_conect_init()
644
{
645
  gmx_conect_t *gc;
646
  
647
  snew(gc,1);
648
  
649
  return (gmx_conect) gc;
650
}
651

    
652
void gmx_conect_done(gmx_conect conect)
653
{
654
  gmx_conect_t *gc = (gmx_conect_t *)conect;
655
  
656
  sfree(gc->conect);
657
}
658

    
659
bool gmx_conect_exist(gmx_conect conect,int ai,int aj)
660
{
661
  gmx_conect_t *gc = (gmx_conect_t *)conect;
662
  int i;
663
  
664
  /* if (!gc->bSorted) 
665
     sort_conect(gc);*/
666
     
667
  for(i=0; (i<gc->nconect); i++) 
668
    if (((gc->conect[i].ai == ai) &&
669
         (gc->conect[i].aj == aj)) ||
670
        ((gc->conect[i].aj == ai) &&
671
         (gc->conect[i].ai == aj)))
672
      return TRUE;
673
  return FALSE;
674
}
675

    
676
void gmx_conect_add(gmx_conect conect,int ai,int aj)
677
{
678
  gmx_conect_t *gc = (gmx_conect_t *)conect;
679
  int i;
680
  
681
  /* if (!gc->bSorted) 
682
     sort_conect(gc);*/
683
  
684
  if (!gmx_conect_exist(conect,ai,aj)) {   
685
    srenew(gc->conect,++gc->nconect);
686
    gc->conect[gc->nconect-1].ai = ai;
687
    gc->conect[gc->nconect-1].aj = aj;
688
  }
689
}
690

    
691
int read_pdbfile(FILE *in,char *title,int *model_nr,
692
                 t_atoms *atoms,rvec x[],int *ePBC,matrix box,bool bChange,
693
                 gmx_conect conect)
694
{
695
    gmx_conect_t *gc = (gmx_conect_t *)conect;
696
    t_symtab symtab;
697
    bool bCOMPND;
698
    bool bConnWarn = FALSE;
699
    char line[STRLEN+1];
700
    int  line_type;
701
    char *c,*d;
702
    int  natom,chainnum,nres_ter_prev=0;
703
    char chidmax=' ';
704
    bool bStop=FALSE;
705

    
706
    if (ePBC) 
707
    {
708
        /* Only assume pbc when there is a CRYST1 entry */
709
        *ePBC = epbcNONE;
710
    }
711
    if (box != NULL) 
712
    {
713
        clear_mat(box);
714
    }
715
    
716
    open_symtab(&symtab);
717

    
718
    bCOMPND=FALSE;
719
    title[0]='\0';
720
    natom=0;
721
    chainnum=0;
722
    while (!bStop && (fgets2(line,STRLEN,in) != NULL)) 
723
    {
724
        line_type = line2type(line);
725
        
726
        switch(line_type) 
727
        {
728
            case epdbATOM:
729
            case epdbHETATM:
730
                natom = read_atom(&symtab,line,line_type,natom,atoms,x,chainnum,bChange);
731
                break;
732
      
733
            case epdbANISOU:
734
                if (atoms->pdbinfo)
735
                {
736
                    read_anisou(line,natom,atoms);
737
                }
738
                break;
739
                
740
            case epdbCRYST1:
741
                read_cryst1(line,ePBC,box);
742
                break;
743
                
744
            case epdbTITLE:
745
            case epdbHEADER:
746
                if (strlen(line) > 6) 
747
                {
748
                    c=line+6;
749
                    /* skip HEADER or TITLE and spaces */
750
                    while (c && (c[0]!=' ')) c++;
751
                    while (c && (c[0]==' ')) c++;
752
                    /* truncate after title */
753
                    d=strstr(c,"      ");
754
                    if (d) 
755
                    {
756
                        d[0]='\0';
757
                    }
758
                    if (strlen(c)>0)
759
                    {
760
                        strcpy(title,c);
761
                    }
762
                }
763
                break;
764
      
765
            case epdbCOMPND:
766
                if ((!strstr(line,": ")) || (strstr(line+6,"MOLECULE:"))) 
767
                {
768
                    if ( !(c=strstr(line+6,"MOLECULE:")) )
769
                    {
770
                        c=line;
771
                    }
772
                    /* skip 'MOLECULE:' and spaces */
773
                    while (c && (c[0]!=' ')) c++;
774
                    while (c && (c[0]==' ')) c++;
775
                    /* truncate after title */
776
                    d=strstr(c,"   ");
777
                    if (d) 
778
                    {
779
                        while ( (d[-1]==';') && d>c)  d--;
780
                        d[0]='\0';
781
                    }
782
                    if (strlen(c) > 0)
783
                    {
784
                        if (bCOMPND) 
785
                        {
786
                            strcat(title,"; ");
787
                            strcat(title,c);
788
                        } 
789
                        else
790
                        {
791
                            strcpy(title,c);
792
                        }
793
                    }
794
                    bCOMPND=TRUE;
795
                }
796
                break;
797
      
798
            case epdbTER:
799
                chainnum++;
800
                break;
801
                
802
            case epdbMODEL:
803
                if(model_nr)
804
                {
805
                    sscanf(line,"%*s%d",model_nr);
806
                }
807
                break;
808

    
809
            case epdbENDMDL:
810
                bStop=TRUE;
811
                break;
812
            case epdbCONECT:
813
                if (gc) 
814
                {
815
                    gmx_conect_addline(gc,line);
816
                }
817
                else if (!bConnWarn)
818
                {
819
                    fprintf(stderr,"WARNING: all CONECT records are ignored\n");
820
                    bConnWarn = TRUE;
821
                }
822
                break;
823
                
824
            default:
825
                break;
826
        }
827
    }
828

    
829
    free_symtab(&symtab);
830
    return natom;
831
}
832

    
833
void get_pdb_coordnum(FILE *in,int *natoms)
834
{
835
    char line[STRLEN];
836
   
837
    *natoms=0;
838
    while (fgets2(line,STRLEN,in)) 
839
    {
840
        if ( strncmp(line,"ENDMDL",6) == 0 ) 
841
        {
842
            break;
843
        }
844
        if ((strncmp(line,"ATOM  ",6) == 0) || (strncmp(line,"HETATM",6) == 0))
845
        {
846
            (*natoms)++;
847
        }
848
    }
849
}
850

    
851
void read_pdb_conf(const char *infile,char *title, 
852
                   t_atoms *atoms,rvec x[],int *ePBC,matrix box,bool bChange,
853
                   gmx_conect conect)
854
{
855
  FILE *in;
856
  
857
  in = gmx_fio_fopen(infile,"r");
858
  read_pdbfile(in,title,NULL,atoms,x,ePBC,box,bChange,conect);
859
  gmx_fio_fclose(in);
860
}
861

    
862
gmx_conect gmx_conect_generate(t_topology *top)
863
{
864
  int f,i;
865
  gmx_conect gc;
866
  
867
  /* Fill the conect records */
868
  gc  = gmx_conect_init();
869

    
870
  for(f=0; (f<F_NRE); f++) {
871
    if (IS_CHEMBOND(f))
872
      for(i=0; (i<top->idef.il[f].nr); i+=interaction_function[f].nratoms+1) {
873
        gmx_conect_add(gc,top->idef.il[f].iatoms[i+1],
874
                       top->idef.il[f].iatoms[i+2]);
875
    }
876
  }
877
  return gc;
878
}