LCOV - code coverage report
Current view: top level - molfile - pdbplugin.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage (other modules) Lines: 87 211 41.2 %
Date: 2024-10-11 08:09:49 Functions: 5 14 35.7 %

          Line data    Source code
       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             : University of Illinois Open Source License
       3             : Copyright 2003 Theoretical and Computational Biophysics Group, 
       4             : All rights reserved.
       5             : 
       6             : Developed by:           Theoretical and Computational Biophysics Group
       7             :                         University of Illinois at Urbana-Champaign
       8             :                         http://www.ks.uiuc.edu/
       9             : 
      10             : Permission is hereby granted, free of charge, to any person obtaining a copy of
      11             : this software and associated documentation files (the Software), to deal with 
      12             : the Software without restriction, including without limitation the rights to 
      13             : use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies 
      14             : of the Software, and to permit persons to whom the Software is furnished to 
      15             : do so, subject to the following conditions:
      16             : 
      17             : Redistributions of source code must retain the above copyright notice, 
      18             : this list of conditions and the following disclaimers.
      19             : 
      20             : Redistributions in binary form must reproduce the above copyright notice, 
      21             : this list of conditions and the following disclaimers in the documentation 
      22             : and/or other materials provided with the distribution.
      23             : 
      24             : Neither the names of Theoretical and Computational Biophysics Group, 
      25             : University of Illinois at Urbana-Champaign, nor the names of its contributors 
      26             : may be used to endorse or promote products derived from this Software without 
      27             : specific prior written permission.
      28             : 
      29             : THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
      30             : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
      31             : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL 
      32             : THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 
      33             : OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
      34             : ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 
      35             : OTHER DEALINGS WITH THE SOFTWARE.
      36             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      37             : #if defined(__PLUMED_HAS_MOLFILE_PLUGINS) && ! defined(__PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS)
      38             : /***************************************************************************
      39             :  *cr
      40             :  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
      41             :  *cr                        University of Illinois
      42             :  *cr                         All Rights Reserved
      43             :  *cr
      44             :  ***************************************************************************/
      45             : 
      46             : /***************************************************************************
      47             :  * RCS INFORMATION:
      48             :  *
      49             :  *      $RCSfile: pdbplugin.c,v $
      50             :  *      $Author: johns $       $Locker:  $             $State: Exp $
      51             :  *      $Revision: 1.73 $       $Date: 2016/11/28 05:01:54 $
      52             :  *
      53             :  ***************************************************************************/
      54             : 
      55             : /*
      56             :  * PDB file format specifications:
      57             :  *   http://www.rcsb.org/pdb/static.do?p=file_formats/pdb/index.html
      58             :  */
      59             : 
      60             : #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
      61             : 
      62             : #include "molfile_plugin.h"
      63             : #include "readpdb.h"
      64             : #include "periodic_table.h"
      65             : #include <stdio.h>
      66             : #include <stdlib.h>
      67             : #include <string.h>
      68             : 
      69             : namespace PLMD{
      70             : namespace molfile{
      71             : 
      72             : /*
      73             :  * API functions start here
      74             :  */
      75             : 
      76             : typedef struct {
      77             :   FILE *fd;
      78             :   int first_frame;
      79             :   int natoms;
      80             :   molfile_atom_t *atomlist;
      81             :   molfile_metadata_t *meta;
      82             :   int nconect;
      83             :   int nbonds, maxbnum;
      84             :   int *from, *to, *idxmap;
      85             : } pdbdata;
      86             : 
      87          48 : static void *open_pdb_read(const char *filepath, const char *filetype, 
      88             :     int *natoms) {
      89             :   FILE *fd;
      90             :   pdbdata *pdb;
      91             :   char pdbstr[PDB_BUFFER_LENGTH];
      92             :   int indx, nconect;
      93             : 
      94          48 :   fd = fopen(filepath, "r");
      95          48 :   if (!fd) 
      96             :     return NULL;
      97          48 :   pdb = (pdbdata *)malloc(sizeof(pdbdata));
      98          48 :   pdb->fd = fd;
      99          48 :   pdb->meta = (molfile_metadata_t *) malloc(sizeof(molfile_metadata_t));
     100             :   memset(pdb->meta, 0, sizeof(molfile_metadata_t));
     101             : 
     102             :   pdb->meta->remarklen = 0;
     103             :   pdb->meta->remarks = NULL;
     104             : 
     105          48 :   *natoms=0;
     106             :   nconect=0;
     107             :   do {
     108      100481 :     indx = read_pdb_record(pdb->fd, pdbstr);
     109      100481 :     if (indx == PDB_ATOM) {
     110      100260 :       *natoms += 1;
     111         221 :     } else if (indx == PDB_CONECT) {
     112           0 :       nconect++;
     113         221 :     } else if (indx == PDB_HEADER) {
     114           0 :       get_pdb_header(pdbstr, pdb->meta->accession, pdb->meta->date, NULL);
     115           0 :       if (strlen(pdb->meta->accession) > 0) 
     116           0 :         strcpy(pdb->meta->database, "PDB");
     117         221 :     } else if (indx == PDB_REMARK || indx == PDB_CONECT || indx == PDB_UNKNOWN) {
     118         142 :       int len=strlen(pdbstr);
     119         142 :       int newlen = len + pdb->meta->remarklen;
     120             : 
     121         142 :       char *newstr=(char*)realloc(pdb->meta->remarks, newlen + 1);
     122         142 :       if (newstr != NULL) {
     123         142 :         pdb->meta->remarks = newstr;
     124         142 :         pdb->meta->remarks[pdb->meta->remarklen] = '\0';
     125         142 :         memcpy(pdb->meta->remarks + pdb->meta->remarklen, pdbstr, len);
     126         142 :         pdb->meta->remarks[newlen] = '\0';
     127         142 :         pdb->meta->remarklen = newlen;
     128             :       }
     129             :     }
     130             :  
     131      100481 :   } while (indx != PDB_END && indx != PDB_EOF);
     132             : 
     133             :   /* If no atoms were found, this is probably not a PDB file! */
     134          48 :   if (!*natoms) {
     135           0 :     fprintf(stderr, "PDB file '%s' contains no atoms.\n", filepath);
     136           0 :     if (pdb->meta->remarks != NULL)
     137           0 :       free(pdb->meta->remarks);
     138           0 :     if (pdb->meta != NULL)
     139           0 :       free(pdb->meta);
     140           0 :     free(pdb);
     141           0 :     return NULL;
     142             :   }
     143             : 
     144          48 :   rewind(pdb->fd); /* if ok, rewind file and prepare to parse it for real */
     145          48 :   pdb->natoms = *natoms;
     146          48 :   pdb->nconect = nconect;
     147          48 :   pdb->nbonds = 0;
     148          48 :   pdb->maxbnum = 0;
     149          48 :   pdb->from = NULL;
     150          48 :   pdb->to = NULL;
     151          48 :   pdb->idxmap = NULL;
     152          48 :   pdb->atomlist = NULL;
     153             : 
     154             : #if defined(VMDUSECONECTRECORDS)
     155             :   /* allocate atom index translation table if we have 99,999 atoms or less */
     156             :   /* and we have conect records to process                                 */
     157          48 :   if (pdb->natoms < 100000 && pdb->nconect > 0) {
     158           0 :     pdb->idxmap = (int *) malloc(100000 * sizeof(int));
     159             :     memset(pdb->idxmap, 0, 100000 * sizeof(int));
     160             :   }
     161             : #endif
     162             :  
     163             :   return pdb; 
     164             : }
     165             : 
     166           0 : static int read_pdb_structure(void *mydata, int *optflags, 
     167             :     molfile_atom_t *atoms) { 
     168             :   pdbdata *pdb = (pdbdata *)mydata;
     169             :   molfile_atom_t *atom;
     170             :   char pdbrec[PDB_BUFFER_LENGTH];
     171             :   int i, rectype, atomserial, pteidx;
     172             :   char ridstr[8];
     173             :   char elementsymbol[3];
     174             :   int badptecount = 0;
     175           0 :   long fpos = ftell(pdb->fd);
     176             : 
     177           0 :   *optflags = MOLFILE_INSERTION | MOLFILE_OCCUPANCY | MOLFILE_BFACTOR |
     178             :               MOLFILE_ALTLOC | MOLFILE_ATOMICNUMBER | MOLFILE_BONDSSPECIAL;
     179             : 
     180             :   i = 0;
     181             :   do {
     182           0 :     rectype = read_pdb_record(pdb->fd, pdbrec);
     183           0 :     switch (rectype) {
     184           0 :     case PDB_ATOM:
     185           0 :       atom = atoms+i;
     186           0 :       get_pdb_fields(pdbrec, strlen(pdbrec), &atomserial, 
     187           0 :           atom->name, atom->resname, atom->chain, atom->segid, 
     188           0 :           ridstr, atom->insertion, atom->altloc, elementsymbol,
     189             :           NULL, NULL, NULL, &atom->occupancy, &atom->bfactor);
     190             : 
     191           0 :       if (pdb->idxmap != NULL && atomserial < 100000) {
     192           0 :         pdb->idxmap[atomserial] = i; /* record new serial number translation */ 
     193             :       }
     194             :  
     195           0 :       atom->resid = atoi(ridstr);
     196             : 
     197             :       /* determine atomic number from the element symbol */
     198           0 :       pteidx = get_pte_idx_from_string(elementsymbol);
     199           0 :       atom->atomicnumber = pteidx;
     200           0 :       if (pteidx != 0) {
     201           0 :         atom->mass = get_pte_mass(pteidx);
     202           0 :         atom->radius = get_pte_vdw_radius(pteidx);
     203             :       } else {
     204           0 :         badptecount++; /* unrecognized element */
     205             :       }
     206             :  
     207           0 :       strcpy(atom->type, atom->name);
     208           0 :       i++;
     209           0 :       break;
     210             : 
     211           0 :     case PDB_CONECT:
     212             :       /* only read CONECT records for structures where we know they can */
     213             :       /* be valid for all of the atoms in the structure                 */
     214           0 :       if (pdb->idxmap != NULL) {
     215           0 :         get_pdb_conect(pdbrec, pdb->natoms, pdb->idxmap, 
     216             :                        &pdb->maxbnum, &pdb->nbonds, &pdb->from, &pdb->to);
     217             :       }
     218             :       break;
     219             : 
     220             :     default:
     221             :       /* other record types are ignored in the structure callback */
     222             :       /* and are dealt with in the timestep callback or elsewhere */
     223             :       break;
     224             :     }
     225           0 :   } while (rectype != PDB_END && rectype != PDB_EOF);
     226             : 
     227           0 :   fseek(pdb->fd, fpos, SEEK_SET);
     228             : 
     229             :   /* if all atoms are recognized, set the mass and radius flags too,  */
     230             :   /* otherwise let VMD guess these for itself using it's own methods  */
     231           0 :   if (badptecount == 0) {
     232           0 :     *optflags |= MOLFILE_MASS | MOLFILE_RADIUS;
     233             :   }
     234             : 
     235           0 :   return MOLFILE_SUCCESS;
     236             : }
     237             : 
     238           0 : static int read_bonds(void *v, int *nbonds, int **fromptr, int **toptr, 
     239             :                       float ** bondorder,int **bondtype, 
     240             :                       int *nbondtypes, char ***bondtypename) {
     241             :   pdbdata *pdb = (pdbdata *)v;
     242             :   
     243           0 :   *nbonds = 0;
     244           0 :   *fromptr = NULL;
     245           0 :   *toptr = NULL;
     246           0 :   *bondorder = NULL; /* PDB files don't have bond order information */
     247           0 :   *bondtype = NULL;
     248           0 :   *nbondtypes = 0;
     249           0 :   *bondtypename = NULL;
     250             : 
     251             : /* The newest plugin API allows us to return CONECT records as 
     252             :  * additional bonds above and beyond what the distance search returns.
     253             :  * Without that feature, we otherwise have to check completeness and
     254             :  * ignore them if they don't look to be fully specified for this molecule */
     255             : #if !defined(MOLFILE_BONDSSPECIAL)
     256             :   if (pdb->natoms >= 100000) {
     257             :     printf("pdbplugin) Warning: more than 99,999 atoms, ignored CONECT records\n");
     258             :     return MOLFILE_SUCCESS;
     259             :   } else if (((float) pdb->nconect / (float) pdb->natoms) <= 0.85) {
     260             :     printf("pdbplugin) Warning: Probable incomplete bond structure specified,\n");
     261             :     printf("pdbplugin)          ignoring CONECT records\n");
     262             :     return MOLFILE_SUCCESS;
     263             :   } else if (pdb->nconect == 0) {
     264             :     return MOLFILE_SUCCESS;
     265             :   }
     266             : #endif
     267             : 
     268           0 :   *nbonds = pdb->nbonds;
     269           0 :   *fromptr = pdb->from;
     270           0 :   *toptr = pdb->to;
     271             : 
     272           0 :   return MOLFILE_SUCCESS;
     273             : }
     274             : 
     275             : 
     276             : /* 
     277             :  * 
     278             :  */
     279         107 : static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) {
     280             :   pdbdata *pdb = (pdbdata *)v;
     281             :   char pdbstr[PDB_BUFFER_LENGTH];
     282             :   int indx, i;
     283             :   float *x, *y, *z;
     284             :   float occup, bfac;
     285         107 :   if (pdb->natoms == 0) 
     286             :     return MOLFILE_ERROR; /* EOF */
     287         107 :   if (ts) {
     288         107 :     x = ts->coords;
     289         107 :     y = x+1;
     290         107 :     z = x+2;
     291             :   } else {
     292             :     x = y = z = 0;
     293             :   } 
     294             :   i = 0;
     295             :   do {
     296      104656 :     indx = read_pdb_record(pdb->fd, pdbstr);
     297      104656 :     if((indx == PDB_END || indx == PDB_EOF) && (i < pdb->natoms)) {
     298             :       return MOLFILE_ERROR;
     299      104608 :     } else if(indx == PDB_ATOM) {
     300      104358 :       if(i++ >= pdb->natoms) {
     301             :         break;      
     302             :       }
     303             :       /* just get the coordinates, and store them */
     304      104358 :       if (ts) {
     305      104358 :         get_pdb_coordinates(pdbstr, x, y, z, &occup, &bfac);
     306      104358 :         x += 3;
     307      104358 :         y += 3;
     308      104358 :         z += 3;
     309             :       } 
     310         250 :     } else if (indx == PDB_CRYST1) {
     311          33 :       if (ts) {
     312          33 :         get_pdb_cryst1(pdbstr, &ts->alpha, &ts->beta, &ts->gamma,
     313             :                                &ts->A, &ts->B, &ts->C);
     314             :       }
     315             :     }
     316      104608 :   } while(!(indx == PDB_END || indx == PDB_EOF));
     317             : 
     318             :   return MOLFILE_SUCCESS;
     319             : }
     320             : 
     321          48 : static void close_pdb_read(void *v) { 
     322             :   pdbdata *pdb = (pdbdata *)v;
     323          48 :   if (pdb->fd != NULL)
     324          48 :     fclose(pdb->fd);
     325          48 :   if (pdb->idxmap != NULL)
     326           0 :     free(pdb->idxmap);
     327          48 :   if (pdb->meta->remarks != NULL)
     328          38 :     free(pdb->meta->remarks);
     329          48 :   if (pdb->meta != NULL) 
     330          48 :     free(pdb->meta);
     331          48 :   free(pdb);
     332          48 : }
     333             : 
     334           0 : static void *open_file_write(const char *path, const char *filetype, 
     335             :     int natoms) {
     336             : 
     337             :   FILE *fd;
     338             :   pdbdata *pdb;
     339           0 :   fd = fopen(path, "w");
     340           0 :   if (!fd) {
     341           0 :     fprintf(stderr, "Unable to open file %s for writing\n", path);
     342           0 :     return NULL;
     343             :   }
     344           0 :   pdb = (pdbdata *)malloc(sizeof(pdbdata));
     345           0 :   pdb->fd = fd;
     346           0 :   pdb->natoms = natoms; 
     347           0 :   pdb->atomlist = NULL;
     348           0 :   pdb->first_frame = 1;
     349           0 :   return pdb;
     350             : }
     351             :  
     352           0 : static int write_structure(void *v, int optflags, 
     353             :     const molfile_atom_t *atoms) {
     354             : 
     355             :   int i;
     356             :   pdbdata *pdb = (pdbdata *)v;
     357           0 :   int natoms = pdb->natoms;
     358           0 :   pdb->atomlist = (molfile_atom_t *)malloc(natoms*sizeof(molfile_atom_t));
     359             :   memcpy(pdb->atomlist, atoms, natoms*sizeof(molfile_atom_t));
     360             : 
     361             :   /* If occ, bfactor, and insertion aren't given, we assign defaultvalues. */
     362           0 :   if (!(optflags & MOLFILE_OCCUPANCY)) {
     363           0 :     for (i=0; i<natoms; i++) pdb->atomlist[i].occupancy = 0.0f;
     364             :   }
     365           0 :   if (!(optflags & MOLFILE_BFACTOR)) {
     366           0 :     for (i=0; i<natoms; i++) pdb->atomlist[i].bfactor= 0.0f;
     367             :   }
     368           0 :   if (!(optflags & MOLFILE_INSERTION)) {
     369           0 :     for (i=0; i<natoms; i++) {
     370           0 :       pdb->atomlist[i].insertion[0] =' ';
     371           0 :       pdb->atomlist[i].insertion[1] ='\0';
     372             :     }
     373             :   }
     374           0 :   if (!(optflags & MOLFILE_ALTLOC)) {
     375           0 :     for (i=0; i<natoms; i++) {
     376           0 :       pdb->atomlist[i].altloc[0]=' ';
     377           0 :       pdb->atomlist[i].altloc[1]='\0';
     378             :     }
     379             :   }
     380           0 :   if (!(optflags & MOLFILE_ATOMICNUMBER)) {
     381           0 :     for (i=0; i<natoms; i++) pdb->atomlist[i].atomicnumber = 0;
     382             :   }
     383             : 
     384             :   /* TODO: put bonds into CONECT records? */
     385           0 :   return MOLFILE_SUCCESS;
     386             : }
     387             : 
     388             : /* SEQRES records look like this:
     389             : 
     390             : COLUMNS        DATA TYPE       FIELD         DEFINITION
     391             : ---------------------------------------------------------------------------------
     392             :  1 -  6        Record name     "SEQRES"
     393             : 
     394             :  9 - 10        Integer         serNum        Serial number of the SEQRES record
     395             :                                              for the current chain.  Starts at 1
     396             :                                              and increments by one each line.
     397             :                                              Reset to 1 for each chain.
     398             : 
     399             : 12             Character       chainID       Chain identifier.  This may be any
     400             :                                              single legal character, including a
     401             :                                              blank which is used if there is
     402             :                                              only one chain.
     403             : 
     404             : 14 - 17        Integer         numRes        Number of residues in the chain.
     405             :                                              This value is repeated on every
     406             :                                              record.
     407             : 
     408             : 20 - 22        Residue name    resName       Residue name.
     409             : 
     410             : 24 - 26        Residue name    resName       Residue name.
     411             : 
     412             : ... and so forth out to 68-70, for a total of 13 in each line (except possibly
     413             : the last.
     414             : 
     415             : source:
     416             : http://www.rcsb.org/pdb/file_formats/pdb/pdbguide2.2/part_35.html
     417             : */
     418             : 
     419             : /*
     420             :  * However, we don't use them right now because of several issues that
     421             :  * can't presently be resolved satisfactorily in VMD:
     422             : 
     423             : According to the RCSB, SEQRES records have to contain all residues, not
     424             : just those in the structure, which means VMD will usually produce incorrect
     425             : output and there's nothing we can do about it.  The RCSB actually specifies
     426             : that all residues in the chain have to present in the SEQRES records, even
     427             : if they're not in the structure.
     428             :   
     429             : We can never know which residues to output.  Our current system of outputting   
     430             : everything is just terrible when you have 20,000 waters in your system; we
     431             : have to fix this immediately.  We could almost get away with making a hash
     432             : table of the names of protein and nucleic acid residues and only write chains
     433             : containing those residues.  However, there's this little snippet from the
     434             : specification:
     435             :   
     436             : * Heterogens which are integrated into the backbone of the chain are listed
     437             :   as being part of the chain and are included in the SEQRES records for
     438             :   that chain.
     439             :   
     440             : That means that we can never know what might appear in the sequence unless we
     441             : also read HET records and keep track of them in VMD as well.  We shouldn't 
     442             : get people depending on such fallible SEQRES records.
     443             :   
     444             : And of course, there's the fact that no other program that we know of besides   
     445             : CE needs these SEQRES records.
     446             : 
     447             :  * Uncomment the write_seqres line in write_timestep to turn them back on.
     448             :  */
     449             : 
     450             : 
     451             : #if 0
     452             : static void write_seqres(FILE * fd, int natoms, const molfile_atom_t *atomlist) {
     453             :   int i=0;
     454             :   while (i < natoms) {
     455             :     int k, serNum;
     456             :     int j = i;
     457             :     int ires, nres = 1;
     458             :     int resid = atomlist[i].resid;
     459             :     /* Count up the number of residues in the chain */
     460             :     const char *chain = atomlist[i].chain;
     461             :     while (j < natoms && !strcmp(chain, atomlist[j].chain)) {
     462             :       if (resid != atomlist[j].resid) {
     463             :         nres++;
     464             :         resid = atomlist[j].resid;
     465             :       }
     466             :       j++;
     467             :     }
     468             :     /* There are nres residues in the chain, from atoms i to j. */
     469             :     serNum = 1;
     470             :     ires = 1;
     471             :     resid = atomlist[i].resid;
     472             :     fprintf(fd, "SEQRES  %2d %c %4d  ",  serNum, chain[0], nres);
     473             :     serNum = 2;
     474             :     fprintf(fd, "%3s ", atomlist[i].resname);
     475             :     for (k=i; k<j; k++) {
     476             :       if (resid != atomlist[k].resid) {
     477             :         resid = atomlist[k].resid;
     478             :         if (!(ires % 13)) {
     479             :           fprintf(fd, "\nSEQRES  %2d %c %4d  ",  serNum, chain[0], nres);
     480             :           serNum++;
     481             :         }
     482             :         fprintf(fd, "%3s ", atomlist[k].resname);
     483             :         ires++;
     484             :       }
     485             :     }
     486             :     i = j;
     487             :     fprintf(fd, "\n");
     488             :   }
     489             : }
     490             : #endif
     491             : 
     492             : /*
     493             : CRYST1 records look like this:
     494             : The CRYST1 record presents the unit cell parameters, space group, and Z value. If the structure was not determined by crystallographic means, CRYST1 simply defines a unit cube. 
     495             : 
     496             : 
     497             : Record Format 
     498             : 
     499             : COLUMNS       DATA TYPE      FIELD         DEFINITION
     500             : -------------------------------------------------------------
     501             :  1 -  6       Record name    "CRYST1"
     502             : 
     503             :  7 - 15       Real(9.3)      a             a (Angstroms).
     504             : 
     505             : 16 - 24       Real(9.3)      b             b (Angstroms).
     506             : 
     507             : 25 - 33       Real(9.3)      c             c (Angstroms).
     508             : 
     509             : 34 - 40       Real(7.2)      alpha         alpha (degrees).
     510             : 
     511             : 41 - 47       Real(7.2)      beta          beta (degrees).
     512             : 
     513             : 48 - 54       Real(7.2)      gamma         gamma (degrees).
     514             : 
     515             : 56 - 66       LString        sGroup        Space group.
     516             : 
     517             : 67 - 70       Integer        z             Z value.
     518             : 
     519             : * If the coordinate entry describes a structure determined by a technique
     520             : other than crystallography, CRYST1 contains a = b = c = 1.0, alpha =
     521             : beta = gamma = 90 degrees, space group = P 1, and Z = 1.
     522             : 
     523             : We will use "P 1" and "1" for space group and z value, as recommended, but
     524             : we'll populate the other fields with the unit cell information we do have.
     525             : 
     526             : */
     527             :   
     528           0 : static void write_cryst1(FILE *fd, const molfile_timestep_t *ts) {
     529           0 :   fprintf(fd, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n", 
     530           0 :     ts->A, ts->B, ts->C, ts->alpha, ts->beta, ts->gamma);
     531           0 : }
     532             : 
     533             : 
     534           0 : static int write_timestep(void *v, const molfile_timestep_t *ts) {
     535             :   pdbdata *pdb = (pdbdata *)v; 
     536             :   const molfile_atom_t *atom;
     537             :   const float *pos;
     538             :   int i;
     539             :   char elementsymbol[3];
     540             : 
     541           0 :   if (pdb->natoms == 0)
     542             :     return MOLFILE_SUCCESS;
     543             : 
     544           0 :   if (pdb->first_frame) {
     545             :     /* Turn off SEQRES writing for now; see comments above.
     546             :     write_seqres(pdb->fd, pdb->natoms, pdb->atomlist);
     547             :     */
     548           0 :     write_cryst1(pdb->fd, ts);
     549           0 :     pdb->first_frame = 0;
     550             :   }
     551           0 :   atom = pdb->atomlist;
     552           0 :   pos = ts->coords;
     553           0 :   for (i=0; i<pdb->natoms; i++) {
     554             :     /*
     555             :      * The 8.3 format for position, occupancy, and bfactor permits values 
     556             :      * only in the range of -999.9994 to 9999.9994 (so that they round
     557             :      * to the range [-999.999, 9999.999]).  If values fall outside of that
     558             :      * range, fail and emit an error message rather than generate a
     559             :      * misformatted PDB file.
     560             :      */
     561             : #define PDBBAD(x) ((x) < -999.9994f || (x) > 9999.9994f)
     562           0 :     if (PDBBAD(pos[0]) || PDBBAD(pos[1]) || PDBBAD(pos[2]) ||
     563           0 :                 PDBBAD(atom->occupancy) || PDBBAD(atom->bfactor)) {
     564           0 :             fprintf(stderr, "PDB WRITE ERROR: Position, occupancy, or b-factor (beta) for atom %d\n", i);
     565           0 :       fprintf(stderr, "                 cannot be written in PDB format.\n");
     566           0 :       fprintf(stderr, "                 File will be truncated.\n");
     567           0 :       return MOLFILE_ERROR;
     568             :     }
     569             : 
     570             :     /* check the atomicnumber and format the atomic element symbol string */
     571           0 :     strcpy(elementsymbol, (atom->atomicnumber < 1) ? "  " : get_pte_label(atom->atomicnumber));
     572           0 :     elementsymbol[0] = toupper(elementsymbol[0]);
     573           0 :     elementsymbol[1] = toupper(elementsymbol[1]);
     574             :  
     575           0 :     if (!write_raw_pdb_record(pdb->fd,  
     576           0 :         "ATOM  ", i+1, atom->name, atom->resname, atom->resid, 
     577           0 :         atom->insertion, atom->altloc, elementsymbol,
     578             :         pos[0], pos[1], pos[2], 
     579           0 :         atom->occupancy, atom->bfactor, atom->chain, atom->segid)) {
     580           0 :       fprintf(stderr, 
     581             :           "PDB: Error encountered writing atom %d; file may be incomplete.\n",
     582             :           i+1);
     583           0 :       return MOLFILE_ERROR;
     584             :     }
     585           0 :     ++atom;
     586           0 :     pos += 3;
     587             :   }
     588           0 :   fprintf(pdb->fd, "END\n");
     589             : 
     590             :   return MOLFILE_SUCCESS;
     591             : }
     592             :  
     593           0 : static void close_file_write(void *v) {
     594             :   pdbdata *pdb = (pdbdata *)v; 
     595           0 :   fclose(pdb->fd);
     596           0 :   free(pdb->atomlist);
     597           0 :   free(pdb);
     598           0 : }
     599             : 
     600           0 : static int read_molecule_metadata(void *v, molfile_metadata_t **metadata) {
     601             :   pdbdata *pdb = (pdbdata *)v; 
     602           0 :   *metadata = pdb->meta;
     603           0 :   return MOLFILE_SUCCESS;
     604             : }
     605             : 
     606             : /*
     607             :  * Initialization stuff down here
     608             :  */
     609             : 
     610             : static molfile_plugin_t plugin;
     611             :  
     612        6946 : VMDPLUGIN_API int VMDPLUGIN_init() {
     613             :   memset(&plugin, 0, sizeof(molfile_plugin_t));
     614        6946 :   plugin.abiversion = vmdplugin_ABIVERSION;
     615        6946 :   plugin.type = MOLFILE_PLUGIN_TYPE;
     616        6946 :   plugin.name = "pdb";
     617        6946 :   plugin.prettyname = "PDB";
     618        6946 :   plugin.author = "Justin Gullingsrud, John Stone";
     619        6946 :   plugin.majorv = 1;
     620        6946 :   plugin.minorv = 16;
     621        6946 :   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
     622        6946 :   plugin.filename_extension = "pdb,ent";
     623        6946 :   plugin.open_file_read = open_pdb_read;
     624        6946 :   plugin.read_structure = read_pdb_structure;
     625        6946 :   plugin.read_bonds = read_bonds;
     626        6946 :   plugin.read_next_timestep = read_next_timestep;
     627        6946 :   plugin.close_file_read = close_pdb_read;
     628        6946 :   plugin.open_file_write = open_file_write;
     629        6946 :   plugin.write_structure = write_structure;
     630        6946 :   plugin.write_timestep = write_timestep;
     631        6946 :   plugin.close_file_write = close_file_write;
     632        6946 :   plugin.read_molecule_metadata = read_molecule_metadata;
     633        6946 :   return VMDPLUGIN_SUCCESS;
     634             : }
     635             : 
     636        6946 : VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
     637        6946 :   (*cb)(v, (vmdplugin_t *)&plugin);
     638        6946 :   return VMDPLUGIN_SUCCESS;
     639             : }
     640             : 
     641           0 : VMDPLUGIN_API int VMDPLUGIN_fini() {
     642           0 :   return VMDPLUGIN_SUCCESS;
     643             : }
     644             : 
     645             : }
     646             : }
     647             : 
     648             : #endif

Generated by: LCOV version 1.15