/* Copyright 1995, Columbia University, all rights reserved. * Permission is granted to utilize and disseminate this code or * document without charge, provided that (1) this copyright notice is * not removed, and (2) all changes made by other than members of the * MacroModel Development Group at Columbia University are, if further * disseminated, (a) noted as such; for example, by means of source-code * comment lines, and (b) communicated back to the author for possible * inclusion in subsequent versions. */ /**************************************************************************** * pss. 2/95. Program to read an mmod file, compressed or otherwise, and * attempt to output as compressed a file as is possible without * reordering atoms or bond tables. Also can uncompress. * * The term "comp" is generally used in variable names and comments to * refer to data stored in a compressed CT. The term "conn" is used * to refer to the data that is stored in a full CT, but not a compressed * CT. * * $RCSfile: mmio_convert.c,v $ * $Revision: 1.10 $ * $Date: 1996/01/24 17:51:58 $ ***************************************************************************/ #if defined( GETOPT_STDLIB ) #include #elif defined( GETOPT_UNISTD ) #include #elif defined( GETOPT_GETOPT ) #include #endif #include #include #include #if defined( DBMALLOC ) #include #endif #include "mmio.h" /* variable #definitions: */ #define FALSE 0 #define TRUE 1 /* ...used in declarations of arrays to be malloc()ed: */ #define AR * /* macro #definitions: */ /* ...automagic exit with appropriate message if error is encountered: */ #define err_check( status ) \ if( status == MMIO_ERR ) { \ fprintf( stderr, \ "!!! mmio_uncompress: exiting, file %s, line %d\n", \ __FILE__, __LINE__ ); \ return -1; \ } /* structure declarations: */ /* ...info readable from full CT but not from compressed CT: */ struct ct_conntag { struct conn_atomtag { /* atom data that don't change with coordinates: */ int itype; /* mmod atom type */ int nbond; /* number of bonded neighbors (n's) */ int bond_atom[ MMIO_MAXBOND ]; /* atom numbers of bonded n's */ int bond_order[ MMIO_MAXBOND ]; /* bond orders to bonded n's */ float charge1; /* 1st charge field in .dat */ float charge2; /* 1st charge field in .dat */ char chain; /* pdb chain designator */ int resnum; /* pdb res number */ char resname1; /* 1-char res name */ char resname4[ MMIO_S_STRLEN + 1 ]; /* 4-char pdb res name */ char pdbname[ MMIO_S_STRLEN + 1 ]; /* 4-char pdb atom name */ } AR atom; int natom; /* no. of atoms in this molecule */ }; /* ...all info readable from compressed CT: */ struct ct_comptag { struct comp_atomtag { /* data per atom: */ float xyz[ 3 ]; /* coords */ int color; /* color */ int mmod_iatom; /* MacroModel atom number -- index origin 1 */ } AR atom; int natom; /* number of atoms stored in this struct */ char title[ MMIO_L_STRLEN + 1 ]; /* title from input file */ }; /* definitions of static variables: */ static int verbose_flag = FALSE; struct ct_conntag global_conn[ 2 ]; struct ct_comptag global_comp[ 2 ]; /* prototypes: */ int main( int narg, char *arg[] ); static void verbose( char *fmt, ... ); static void usage( void ); static int diff_conn( struct ct_conntag *conn1, struct ct_conntag *conn2 ); static int get_ct( int iread, int natom, char *title, struct ct_conntag *newconn, struct ct_comptag *newcomp ); static int put_ct( int iwrite, struct ct_conntag *conn, struct ct_comptag *oldcomp, struct ct_comptag *newcomp ); static int filldiff( int AR diff, struct ct_comptag *oldcomp, struct ct_comptag *newcomp ); /************************************************************************/ int main( int narg, char *arg[] ) { int status, ct_type; struct ct_conntag *oldconn, *newconn, *tmpconn; struct ct_comptag *oldcomp, *newcomp, *tmpcomp; int ict= -1, natom; char title[ MMIO_L_STRLEN ]; char readfname[200], writefname[200]; int ct_type_requested = MMIO_COMPRESSED; int opt; extern int optind; int iread, iwrite; /* initializations: */ global_conn[0].atom = global_conn[1].atom = NULL; global_comp[0].atom = global_comp[1].atom = NULL; global_conn[0].natom = global_conn[1].natom = global_comp[0].natom = global_comp[1].natom = 0; oldconn = global_conn + 0; newconn = global_conn + 1; oldcomp = global_comp + 0; newcomp = global_comp + 1; /* parse cmdline options: */ while( (opt=getopt(narg,arg,"cfv")) != -1 ) { switch( opt ) { case 'c': /* write compressed CT's where possible: */ ct_type_requested = MMIO_COMPRESSED; break; case 'f': /* always write full CT's: */ ct_type_requested = MMIO_FULL; break; case 'v': /* put play-by-play on stderr: */ verbose_flag = TRUE; break; default: fprintf( stderr, "!!! mmio_convert main: found" " illegal flag '%c'\n", opt ); usage(); return -1; } } /* see how many extra args there are, after optional flags: */ switch( narg-optind ) { case 0: /* stdin for input, stdout for output: */ strcpy( readfname, "-" ); strcpy( writefname, "-" ); break; case 1: /* arg for input fname, stdout for output: */ strcpy( readfname, arg[optind] ); strcpy( writefname, "-" ); break; case 2: /* args for input & output fnames: */ strcpy( readfname, arg[optind++] ); strcpy( writefname, arg[optind++] ); break; default: /* huh? */ fprintf( stderr, "!!! mmio_convert main: found %d" " args after flags\n", narg-optind ); usage(); return -1; } /* tell the user what he asked for: */ verbose( "mmio_convert main: output format= %s\n", mmio_return_code(ct_type_requested) ); verbose( "mmio_convert main: readfname= %s\n", strcmp("-",readfname)==0?"":readfname ); verbose( "mmio_convert main: writefname= %s\n", strcmp("-",writefname)==0?"":writefname ); /* tell mmio library to direct its errmsgs (if any) to stderr: */ mmio_errfile( stderr ); /* open file for reading: */ status = mmio_open( &iread, readfname, MMIO_READ ); verbose( "mmio_convert main: mmio_open() returns %s for readfile\n", mmio_return_code(status) ); err_check( status ); /* open file for writing: */ status = mmio_open( &iwrite, writefname, MMIO_WRITE ); verbose( "mmio_convert main: mmio_open() returns %s for writefile\n", mmio_return_code(status) ); err_check( status ); /* get CT's from input file until EOF is detected: */ while( (ct_type=mmio_get_ct(iread,ct_type_requested,&natom,title)) != MMIO_EOF ) { ict += 1; verbose( "mmio_convert main: mmio_get_ct() returns %s " "for ict %d; natom= %d, title=\n'%s'\n", mmio_return_code(ct_type), ict, natom, title ); err_check( ct_type ); if( ct_type == MMIO_FULL ) { /* read full CT into "new" data structures: */ status = get_ct( iread, natom, title, newconn, newcomp ); verbose( "mmio_convert main: get_ct() returns %s\n", mmio_return_code(status) ); err_check( status ); if( ict==0 || ct_type_requested==MMIO_FULL || diff_conn(newconn,oldconn)) { /* this is either the first CT, or the * connectivity differs from the previous * CT, or we have requested full CT's; * write out the full CT we just read: */ status = put_ct( iwrite, newconn, newcomp, NULL ); verbose( "mmio_convert main: put_ct() returns" " %s\n", mmio_return_code(status) ); err_check( status ); /* now make the CT we just read the "old" * CT, and prepare to read the next one * into the "new" storage: */ tmpconn = oldconn; oldconn = newconn; newconn = tmpconn; tmpcomp = oldcomp; oldcomp = newcomp; newcomp = tmpcomp; continue; } /* If the conn part of the new CT is the same * as that of the old one, and we have requested * compressed output, then at this moment the conn * is in oldconn and the comp is in newcomp. */ } else if( ct_type == MMIO_COMPRESSED ) { /* read compressed data into newcomp; the conn data * remains untouched in oldconn: */ status = get_ct( iread, natom, title, NULL, newcomp ); verbose( "mmio_convert main: get_ct() returns %s\n", mmio_return_code(status) ); err_check( status ); } else { /* unexpected ct_type */ fprintf( stderr, "!!! mmio_uncompress: mmio_get_ct()" " returns unexpected status %s\n", mmio_return_code(status) ); err_check( MMIO_ERR ); } /* write out a new CT, taking comp info from newcomp * where it overrides the oldcomp info: */ status = put_ct( iwrite, NULL, oldcomp, newcomp ); verbose( "mmio_convert main: put_ct() returns" " %s\n", mmio_return_code(status) ); err_check( status ); } /* this output covers MMIO_EOF returns from mmio_get_ct(): */ verbose( "mmio_convert main: mmio_get_ct() returns %s " "for ict %d\n", mmio_return_code(ct_type), ict ); ict += 1; status = mmio_cleanup(); verbose( "mmio_convert main: mmio_cleanup() returns" " %s\n", mmio_return_code(status) ); err_check( status ); /* announce the number of CT's we processed: */ verbose( "mmio_convert main: %d ct's found in input file\n", ict ); return ict; } /************************************************************************/ static void verbose( char *fmt, ... ) /* if verbose_flag is TRUE, print msg to stderr; otherwise, do nothing: */ { char buf[ 200 ]; va_list args; if( verbose_flag == FALSE ) { return; } va_start( args, fmt ); ( void )vsprintf( buf, fmt, args ); va_end( args ); fprintf( stderr, "%s", buf ); } /************************************************************************/ static void usage( void ) /* simply print usage message: */ { fprintf( stderr, "Usage: mmio_convert [-cfv] [infile [outfile]]\n" ); } /************************************************************************/ static int diff_conn( struct ct_conntag *conn1, struct ct_conntag *conn2 ) /* return TRUE if conn1 and conn2 are different; FALSE otherwise: */ { int iatom, natom, ibond, nbond; struct conn_atomtag *conn1_atom, *conn2_atom; /* the CT's are different if they have different numbers of atoms: */ if( (natom=conn1->natom) != conn2->natom ) { return TRUE; } /* compare relevant data for each atom: */ for( iatom=0; iatomatom + iatom; conn2_atom = conn2->atom + iatom; /* the CT's are different if any atom has different numbers of * bonds in the two CT's: */ if( (nbond=conn1_atom->nbond) != conn2_atom->nbond ) { return TRUE; } for( ibond=0; ibondbond_atom[ibond] != conn2_atom->bond_atom[ibond] || conn1_atom->bond_order[ibond] != conn2_atom->bond_order[ibond] ) { return TRUE; } } } return FALSE; } /************************************************************************/ static int get_ct( int iread, int natom, char *title, struct ct_conntag *conn, struct ct_comptag *comp ) /* call mmio routines to read conn and comp data into structures conn * and comp. * if conn is NULL, then we're expecting compressed info, and we fill * only conn: */ { int iatom, status; struct conn_atomtag *connatom; struct comp_atomtag *compatom; /* fill statically declared comp variables: */ comp->natom = natom; strcpy( comp->title, title ); /* malloc space in comp for atom data: */ if( comp->atom == NULL ) { comp->atom = ( struct comp_atomtag * )malloc( natom*sizeof(struct comp_atomtag) ); } else { comp->atom = ( struct comp_atomtag * )realloc( comp->atom, natom*sizeof(struct comp_atomtag) ); } if( comp->atom == NULL ) { fprintf( stderr, "!!! get_ct: malloc() fails for comp, natom= %d\n", natom ); return MMIO_ERR; } /* fill comp; if conn is needed as well, malloc storage in it for * atoms, and fill it, too: */ if( conn == NULL ) { /* no conn data needed (compressed CT being passed to us): */ /* ...get data atom by atom: */ for( iatom=0; iatomatom + iatom; status = mmio_get_atom( iread, &compatom->mmod_iatom, NULL, NULL, NULL, NULL, compatom->xyz, NULL, NULL, NULL, &compatom->color, NULL, NULL, NULL, NULL ); verbose( "get_ct: mmio_get_atom() returns %s," " iatom= %d\n", mmio_return_code(status), iatom ); if( status == MMIO_ERR ) break; } } else { /* conn data needed (full CT being passed to us): */ /* ...fill statically declared conn variables: */ conn->natom = natom; /* ...malloc space in conn for atom data: */ if( conn->atom == NULL ) { conn->atom = ( struct conn_atomtag * )malloc( natom*sizeof(struct conn_atomtag) ); } else { conn->atom = ( struct conn_atomtag * )realloc( conn->atom, natom*sizeof(struct conn_atomtag) ); } if( conn->atom == NULL ) { fprintf( stderr, "!!! get_ct: malloc() fails for conn," " natom= %d\n", natom ); return MMIO_ERR; } /* ...get data atom by atom: */ for( iatom=0; iatomatom + iatom; connatom = conn->atom + iatom; status = mmio_get_atom( iread, &compatom->mmod_iatom, &conn->atom[iatom].itype, &connatom->nbond, connatom->bond_atom, connatom->bond_order, compatom->xyz, &connatom->charge1, &connatom->charge2, &connatom->chain, &compatom->color, &connatom->resnum, &connatom->resname1, connatom->resname4, connatom->pdbname ); verbose( "get_ct: mmio_get_atom() returns %s," " iatom= %d\n", mmio_return_code(status), iatom ); if( status == MMIO_ERR ) break; } } return status; } /************************************************************************/ static int put_ct( int iwrite, struct ct_conntag *conn, struct ct_comptag *oldcomp, struct ct_comptag *newcomp ) /* call mmio routines to write out a CT. * If conn is NULL, then we will be writing a compressed CT; in this event, * we write lines only for those atoms in newcomp which have coordinates * or colors differing from those of the corresponding atoms in oldcomp. * If conn is not NULL, then we'll write a full CT, using connectivity data * from conn and coords and colors from oldcomp. */ { int status, natom, iatom, idiff; struct conn_atomtag *conn_atom; struct comp_atomtag *comp_atom; static int AR diff=NULL, ndiff=0; if( conn == NULL ) { /* write compressed CT, only for atoms where newcomp * differs from oldcomp; use newcomp data: */ if( ndiff != newcomp->natom ) { ndiff = newcomp->natom; if( diff == NULL ) { diff = (int *)malloc( ndiff*sizeof(int) ); } else { diff = (int *)realloc( diff, ndiff*sizeof(int)); } if( diff == NULL ) { fprintf( stderr, "!!! put_ct: could not malloc diff\n"); return MMIO_ERR; } } /* natom will be the number of atom records written: */ if( (natom=filldiff(diff,oldcomp,newcomp)) < 0 ) { fprintf( stderr, "!!! put_ct: filldiff() fails\n" ); return MMIO_ERR; } status = mmio_put_ct( iwrite, MMIO_COMPRESSED, natom, newcomp->title); verbose( "put_ct: mmio_put_ct() returns %s\n", mmio_return_code(status) ); if( status == MMIO_ERR ) { fprintf( stderr, "!!! put_ct: mmio_put_ct() fails\n" ); return status; } /* traverse the diff list, writing an atom record wherever * diff[] is true; note that since we are writing a * compressed CT, we can pass junk values to mmio_put_atom * for variables that are associated with connectivity: */ for( idiff=0; idiffatom + idiff; status = mmio_put_atom( iwrite, comp_atom->mmod_iatom, 0, 0, NULL, NULL, comp_atom->xyz, 0, 0, '\0', comp_atom->color, 0, '\0', NULL, NULL ); verbose( "put_ct: mmio_put_atom() returns %s," " iatom= %d\n", mmio_return_code(status), iatom ); if( status == MMIO_ERR ) { fprintf( stderr, "!!! put_ct: mmio_put_atom() fails\n"); return status; } } } } else { /* write full CT, taking connectivity data from conn and * other data from oldcomp: */ if( (natom=conn->natom) != oldcomp->natom ) { fprintf( stderr, "!!! put_ct: conn and comp have" " different numbers of atoms\n" ); return MMIO_ERR; } status = mmio_put_ct( iwrite, MMIO_FULL, natom, oldcomp->title); verbose( "put_ct: mmio_put_ct() returns %s\n", mmio_return_code(status) ); if( status == MMIO_ERR ) { fprintf( stderr, "!!! put_ct: mmio_put_ct() fails\n" ); return status; } for( iatom=0; iatomatom + iatom; comp_atom = oldcomp->atom + iatom; status = mmio_put_atom( iwrite, comp_atom->mmod_iatom, conn_atom->itype, conn_atom->nbond, conn_atom->bond_atom, conn_atom->bond_order, comp_atom->xyz, conn_atom->charge1, conn_atom->charge2, conn_atom->chain, comp_atom->color, conn_atom->resnum, conn_atom->resname1, conn_atom->resname4, conn_atom->pdbname ); verbose( "put_ct: mmio_put_atom() returns %s, iatom= %d\n", mmio_return_code(status), iatom ); if( status == MMIO_ERR ) { fprintf( stderr, "!!! put_ct: mmio_put_atom() fails\n" ); return status; } } } return MMIO_OK; } /*******************************************************************/ static int filldiff( int AR diff, struct ct_comptag *oldcomp, struct ct_comptag *newcomp ) /* array diff[] has the same number of elements as newcomp has atoms. * put TRUE in diff in positions where old and new atoms differ in * coords or color; FALSE where they don't; * return the number of places in which there are differences; * this is the number of atom lines in the compressed CT which * the calling function will write. * oldcomp is expected to contain data on all the atoms in the CT, * and thus must have at least as many atoms as newcomp: */ { int ndiff = 0; int iatom, natom = newcomp->natom; struct comp_atomtag *oldatom = oldcomp->atom; struct comp_atomtag *end_oldatom = oldcomp->atom + oldcomp->natom; struct comp_atomtag *newatom = newcomp->atom; /* traverse the newatom list: */ for( iatom=0; iatommmod_iatom != newatom->mmod_iatom ) { oldatom += 1; if( oldatom >= end_oldatom ) { fprintf( stderr, "!!! filldiff: oldatom end" " encountered\n" ); return -1; } } /* now compare newatom and oldatom data: */ if( memcmp(oldatom->xyz,newatom->xyz,3*sizeof(float))!=0 || oldatom->color != newatom->color ) { diff[ iatom ] = TRUE; ndiff += 1; } else { diff[ iatom ] = FALSE; } } return ndiff; } /*******************************************************************/