Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

gromacsplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: gromacsplugin.C,v $
00013  *      $Author: dgomes $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.55 $       $Date: 2024/05/16 14:03:31 $
00015  *
00016  ***************************************************************************/
00017 
00018 #include "largefiles.h"   /* platform dependent 64-bit file I/O defines */
00019 
00020 #include <math.h>
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024 #include <ctype.h>
00025 #include "Gromacs.h"
00026 #include "molfile_plugin.h"
00027 
00028 #if defined(_AIX)
00029 #include <strings.h>
00030 #endif
00031 
00032 #ifndef M_PI
00033 #define M_PI           3.14159265358979323846
00034 #endif
00035 
00036 #if defined(WIN32) || defined(WIN64)
00037 #define strcasecmp stricmp
00038 #endif
00039 
00040 typedef struct {
00041   md_file *mf;
00042   int natoms;
00043   int step;
00044   float timeval;
00045   molfile_atom_t *atomlist;
00046   molfile_metadata_t *meta;
00047 } gmxdata;
00048 
00049 static void convert_vmd_box_for_writing(const molfile_timestep_t *ts, float *x, float *y, float *z)
00050 {
00051 //     const float sa = sin((double)ts->alpha/180.0*M_PI);
00052     const float ca = cos((double)ts->alpha/180.0*M_PI);
00053     const float cb = cos((double)ts->beta/180.0*M_PI);
00054     const float cg = cos((double)ts->gamma/180.0*M_PI);
00055     const float sg = sin((double)ts->gamma/180.0*M_PI);
00056 
00057     x[0] = ts->A / ANGS_PER_NM;
00058     y[0] = 0.0;
00059     z[0] = 0.0;
00060     x[1] = ts->B*cg / ANGS_PER_NM; // ts->B*ca when writing trr?!
00061     y[1] = ts->B*sg / ANGS_PER_NM; // ts->B*sa when writing trr?!
00062     z[1] = 0.0;
00063     x[2] = ts->C*cb / ANGS_PER_NM;
00064     y[2] = (ts->C / ANGS_PER_NM)*(ca - cb*cg)/sg;
00065     z[2] = (ts->C / ANGS_PER_NM)*sqrt((double)(1.0 + 2.0*ca*cb*cg
00066                                - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
00067 }
00068 
00069 static void *open_gro_read(const char *filename, const char *,
00070     int *natoms) {
00071 
00072     md_file *mf;
00073     md_header mdh;
00074     gmxdata *gmx;
00075 
00076     mf = mdio_open(filename, MDFMT_GRO);
00077     if (!mf) {
00078         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00079                 filename, mdio_errmsg(mdio_errno()));
00080         return NULL;
00081     }
00082 
00083     // read in the header data (careful not to rewind!)
00084     if (gro_header(mf, mdh.title, MAX_MDIO_TITLE,
00085     &mdh.timeval, &mdh.natoms, 0) < 0) {
00086         fprintf(stderr, "gromacsplugin) Cannot read header fromm '%s', %s\n",
00087                 filename, mdio_errmsg(mdio_errno()));
00088             // XXX should free the file handle...
00089         return NULL;
00090     }
00091     *natoms = mdh.natoms;
00092     gmx = new gmxdata;
00093     memset(gmx,0,sizeof(gmxdata));
00094     gmx->mf = mf;
00095     gmx->natoms = mdh.natoms;
00096     gmx->meta = new molfile_metadata_t;
00097     memset(gmx->meta,0,sizeof(molfile_metadata_t));
00098     strncpy(gmx->meta->title, mdh.title, 80);
00099     gmx->timeval = mdh.timeval;
00100     return gmx;
00101 }
00102 
00103 static int read_gro_structure(void *mydata, int *optflags,
00104     molfile_atom_t *atoms) {
00105 
00106   md_atom ma;
00107   char buf[MAX_GRO_LINE + 1];
00108   gmxdata *gmx = (gmxdata *)mydata;
00109 
00110   *optflags = MOLFILE_NOOPTIONS; // no optional data
00111 
00112   // read in each atom and add it into the molecule
00113   for (int i = 0; i < gmx->natoms; i++) {
00114     molfile_atom_t *atom = atoms+i;
00115     if (gro_rec(gmx->mf, &ma) < 0) {
00116       fprintf(stderr, "gromacsplugin) Error reading atom %d from file, %s\n",
00117               i+1, mdio_errmsg(mdio_errno()));
00118       return MOLFILE_ERROR;
00119     }
00120     strcpy(atom->name, ma.atomname);
00121     strcpy(atom->type, ma.atomname);
00122     strcpy(atom->resname, ma.resname);
00123     atom->resid = atoi(ma.resid);
00124     atom->chain[0] = '\0';
00125     atom->segid[0] = '\0';
00126   }
00127 
00128   if (mdio_readline(gmx->mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
00129     fprintf(stderr, "gromacsplugin) Warning, error reading box, %s\n",
00130             mdio_errmsg(mdio_errno()));
00131   }
00132 
00133   rewind(gmx->mf->f);
00134   return MOLFILE_SUCCESS;
00135 }
00136 
00137 static int read_gro_molecule_metadata(void *v, molfile_metadata_t **metadata) {
00138   gmxdata *gmx = (gmxdata *)v;
00139   *metadata = gmx->meta;
00140   return MOLFILE_SUCCESS;
00141 }
00142 
00143 static int read_gro_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00144   gmxdata *gmx = (gmxdata *)v;
00145   md_ts mdts;
00146   memset(&mdts, 0, sizeof(md_ts));
00147   mdts.natoms = natoms;
00148 
00149   if (mdio_timestep(gmx->mf, &mdts) < 0)
00150     return MOLFILE_ERROR;
00151   if (ts) {
00152     memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
00153     if (mdts.box) {
00154       ts->A = mdts.box->A;
00155       ts->B = mdts.box->B;
00156       ts->C = mdts.box->C;
00157       ts->alpha = mdts.box->alpha;
00158       ts->beta = mdts.box->beta;
00159       ts->gamma = mdts.box->gamma;
00160     }
00161   }
00162   mdio_tsfree(&mdts);
00163   return MOLFILE_SUCCESS;
00164 }
00165 
00166 static void close_gro_read(void *v) {
00167   gmxdata *gmx = (gmxdata *)v;
00168   mdio_close(gmx->mf);
00169   delete gmx->meta;
00170   delete gmx;
00171 }
00172 
00173 // open file for writing
00174 static void *open_gro_write(const char *filename, const char *filetype,
00175     int natoms) {
00176 
00177     md_file *mf;
00178     gmxdata *gmx;
00179 
00180     mf = mdio_open(filename, MDFMT_GRO, MDIO_WRITE);
00181     if (!mf) {
00182         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00183                 filename, mdio_errmsg(mdio_errno()));
00184         return NULL;
00185     }
00186     gmx = new gmxdata;
00187     memset(gmx,0,sizeof(gmxdata));
00188     gmx->mf = mf;
00189     gmx->natoms = natoms;
00190     gmx->step   = 0;
00191     gmx->meta = new molfile_metadata_t;
00192     memset(gmx->meta,0,sizeof(molfile_metadata_t));
00193     gmx->meta->title[0] = '\0';
00194 
00195     return gmx;
00196 }
00197 
00198 static int write_gro_structure(void *v, int optflags,
00199     const molfile_atom_t *atoms) {
00200 
00201   gmxdata *gmx = (gmxdata *)v;
00202   int natoms = gmx->natoms;
00203   gmx->atomlist = (molfile_atom_t *)malloc(natoms*sizeof(molfile_atom_t));
00204   memcpy(gmx->atomlist, atoms, natoms*sizeof(molfile_atom_t));
00205 
00206   return MOLFILE_SUCCESS;
00207 }
00208 
00209 static int write_gro_timestep(void *v, const molfile_timestep_t *ts) {
00210   gmxdata *gmx = (gmxdata *)v;
00211   const molfile_atom_t *atom;
00212   const float *pos, *vel;
00213   float x[3], y[3], z[3];
00214   int i;
00215 
00216   if (gmx->natoms == 0)
00217     return MOLFILE_SUCCESS;
00218 
00219   atom = gmx->atomlist;
00220   pos = ts->coords;
00221   vel = ts->velocities;
00222 
00223   /* The title cannot be written */
00224 /*  fprintf(gmx->mf->f, "%s", gmx->meta->title);*/
00225   /* Write a dummy title instead */
00226   fprintf(gmx->mf->f, "generated by VMD");
00227 #if vmdplugin_ABIVERSION > 10
00228   fprintf(gmx->mf->f, ", t= %f", ts->physical_time);
00229 #endif
00230   fprintf(gmx->mf->f, "\n");
00231 
00232   fprintf(gmx->mf->f, "%d\n", gmx->natoms);
00233   for (i=0; i<gmx->natoms; i++)
00234   {
00235      fprintf(gmx->mf->f, "%5d%-5s%5s%5d%8.3f%8.3f%8.3f",
00236              atom->resid, atom->resname, atom->name,
00237              (i+1) % 100000, // GRO format only supports indices up to 99999
00238                              // but since GROMACS ignores indices, modular
00239                              // arithmetic prevents formatting problems for 
00240                              // very large structures
00241              pos[0] / ANGS_PER_NM, pos[1] / ANGS_PER_NM, pos[2] / ANGS_PER_NM);
00242      if(vel)
00243      {
00244          fprintf(gmx->mf->f, "%8.4f%8.4f%8.4f", vel[0] / ANGS_PER_NM, vel[1] / ANGS_PER_NM, vel[2] / ANGS_PER_NM);
00245          vel += 3;
00246      }
00247      fprintf(gmx->mf->f, "\n");
00248      ++atom;
00249      pos += 3;
00250   }
00251   convert_vmd_box_for_writing(ts, x, y, z);
00252   fprintf(gmx->mf->f, "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n", x[0], y[1], z[2], y[0], z[0], x[1], z[1], x[2], y[2]);
00253 
00254   return MOLFILE_SUCCESS;
00255 }
00256 
00257 static void close_gro_write(void *v) {
00258   gmxdata *gmx = (gmxdata *)v;
00259   mdio_close(gmx->mf);
00260   free(gmx->atomlist);
00261   delete gmx->meta;
00262   delete gmx;
00263 }
00264 
00265 
00266 static void *open_g96_read(const char *filename, const char *,
00267     int *natoms) {
00268 
00269     md_file *mf;
00270     md_header mdh;
00271     char gbuf[MAX_G96_LINE + 1];
00272 
00273     mf = mdio_open(filename, MDFMT_G96);
00274     if (!mf) {
00275         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00276                 filename, mdio_errmsg(mdio_errno()));
00277         return NULL;
00278     }
00279 
00280         // read in the header data
00281         if (g96_header(mf, mdh.title, MAX_MDIO_TITLE, &mdh.timeval) < 0) {
00282             fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
00283                     filename, mdio_errmsg(mdio_errno()));
00284             return NULL;
00285         }
00286 
00287         // First, look for a timestep block
00288         if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
00289             fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
00290                     filename, mdio_errmsg(mdio_errno()));
00291             return NULL;
00292         }
00293         if (!strcasecmp(gbuf, "TIMESTEP")) {
00294             // Read in the value line and the END line, and the next
00295             if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
00296                 mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
00297                 mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
00298               fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
00299                       filename, mdio_errmsg(mdio_errno()));
00300               return NULL;
00301             }
00302         }
00303         if (strcasecmp(gbuf, "POSITION") && strcasecmp(gbuf, "REFPOSITION")) {
00304           fprintf(stderr, "gromacsplugin) No structure information in file %s\n", filename);
00305           return NULL;
00306         }
00307         *natoms = g96_countatoms(mf);
00308 
00309         gmxdata *gmx = new gmxdata;
00310         memset(gmx,0,sizeof(gmxdata));
00311         gmx->mf = mf;
00312         gmx->natoms = *natoms; 
00313         return gmx;
00314 }
00315 
00316 static int read_g96_structure(void *mydata, int *optflags,
00317     molfile_atom_t *atoms) {
00318 
00319     char gbuf[MAX_G96_LINE + 1];
00320     gmxdata *gmx = (gmxdata *)mydata;
00321     md_atom ma;
00322     md_file *mf = gmx->mf;
00323 
00324     *optflags = MOLFILE_NOOPTIONS; // no optional data
00325 
00326         for (int i = 0; i < gmx->natoms; i++) {
00327             molfile_atom_t *atom = atoms+i;
00328             if (g96_rec(mf, &ma) < 0) {
00329                 fprintf(stderr, "gromacsplugin) Error reading atom %d from file, %s\n",
00330                   i+1, mdio_errmsg(mdio_errno()));
00331                 return MOLFILE_ERROR;
00332             }
00333             strcpy(atom->name, ma.atomname);
00334             strcpy(atom->type, ma.atomname);
00335             strcpy(atom->resname, ma.resname);
00336             atom->resid = atoi(ma.resid);
00337             atom->chain[0] = '\0';
00338             atom->segid[0] = '\0';
00339         }
00340 
00341         if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
00342             fprintf(stderr, "gromacsplugin) Warning, error reading END record, %s\n",
00343                 mdio_errmsg(mdio_errno()));
00344         }
00345 
00346             // ... another problem: there may or may not be a VELOCITY
00347             // block or a BOX block, so we need to read one line beyond
00348             // the POSITION block to determine this. If neither VEL. nor
00349             // BOX are present we've read a line too far and infringed
00350             // on the next timestep, so we need to keep track of the
00351             // position now for a possible fseek() later to backtrack.
00352             long fpos = ftell(mf->f);
00353 
00354             // Now we must read in the velocities and the box, if present
00355             if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) >= 0) {
00356 
00357                 // Is there a velocity block present ?
00358                 if (!strcasecmp(gbuf, "VELOCITY") || !strcasecmp(gbuf, "VELOCITYRED")) {
00359                         // Ignore all the coordinates - VMD doesn't use them
00360                         for (;;) {
00361                                 if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
00362                                         return MOLFILE_ERROR;
00363                                 if (!strcasecmp(gbuf, "END")) break;
00364                         }
00365 
00366                         // Again, record our position because we may need
00367                         // to fseek here later if we read too far.
00368                         fpos = ftell(mf->f);
00369 
00370                         // Go ahead and read the next line.
00371                         if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
00372                     return MOLFILE_ERROR;
00373                 }
00374 
00375                 // Is there a box present ?
00376                 if (!strcasecmp(gbuf, "BOX")) {
00377                         // Ignore the box coordinates at this time.
00378                         if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
00379                     return MOLFILE_ERROR;
00380                         if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
00381                     return MOLFILE_ERROR;
00382                         if (strcasecmp(gbuf, "END"))
00383                     return MOLFILE_ERROR;
00384                 }
00385                 else {
00386                         // We have read too far, so fseek back to the
00387                         // last known safe position so we don't return
00388                         // with the file pointer set infringing on the
00389                         // next timestep data.
00390                         fseek(mf->f, fpos, SEEK_SET);
00391                 }
00392         }
00393         else {
00394             // Go ahead and rewind for good measure
00395             fseek(mf->f, fpos, SEEK_SET);
00396         }
00397         rewind(mf->f);
00398         return MOLFILE_SUCCESS;
00399 }
00400 
00401 static int read_g96_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00402 
00403   gmxdata *gmx = (gmxdata *)v;
00404   md_ts mdts;
00405   memset(&mdts, 0, sizeof(md_ts));
00406   mdts.natoms = natoms;
00407 
00408   if (mdio_timestep(gmx->mf, &mdts) < 0)
00409     return MOLFILE_ERROR;
00410   if (ts) {
00411     memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
00412     if (mdts.box) {
00413       ts->A = mdts.box->A;
00414       ts->B = mdts.box->B;
00415       ts->C = mdts.box->C;
00416       ts->alpha = mdts.box->alpha;
00417       ts->beta = mdts.box->beta;
00418       ts->gamma = mdts.box->gamma;
00419     }
00420   }
00421   mdio_tsfree(&mdts);
00422   return MOLFILE_SUCCESS;
00423 }
00424 
00425 static void close_g96_read(void *v) {
00426   gmxdata *gmx = (gmxdata *)v;
00427   mdio_close(gmx->mf);
00428   delete gmx;
00429 }
00430 
00431 
00432 //
00433 // TRR and XTC files
00434 //
00435 
00436 static void *open_trr_read(const char *filename, const char *filetype,
00437     int *natoms) {
00438 
00439     md_file *mf;
00440     md_header mdh;
00441     gmxdata *gmx;
00442     int format;
00443 
00444     if (!strcmp(filetype, "trr"))
00445       format = MDFMT_TRR;
00446     else if (!strcmp(filetype, "trj"))
00447       format = MDFMT_TRJ;
00448     else if (!strcmp(filetype, "xtc"))
00449       format = MDFMT_XTC;
00450     else
00451       return NULL;
00452 
00453     mf = mdio_open(filename, format);
00454     if (!mf) {
00455         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00456                 filename, mdio_errmsg(mdio_errno()));
00457         return NULL;
00458     }
00459     if (mdio_header(mf, &mdh) < 0) {
00460         mdio_close(mf);
00461         fprintf(stderr, "gromacsplugin) Cannot read header fromm '%s', %s\n",
00462                 filename, mdio_errmsg(mdio_errno()));
00463         return NULL;
00464     }
00465     *natoms = mdh.natoms;
00466     gmx = new gmxdata;
00467     memset(gmx,0,sizeof(gmxdata));
00468     gmx->mf = mf;
00469     gmx->natoms = mdh.natoms;
00470     return gmx;
00471 }
00472 
00473 static int read_trr_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00474   gmxdata *gmx = (gmxdata *)v;
00475   md_ts mdts;
00476   memset(&mdts, 0, sizeof(md_ts));
00477   mdts.natoms = natoms;
00478 
00479   //Take advantage of boolean short-circuiting.
00480   //If ts is null, we want to skip the timestep and call mdio_timeskip.
00481   //Otherwse, read in the full timestep with mdio_timestep
00482   //Error handling is shared
00483   if ((ts == NULL && (mdio_timeskip(gmx->mf, &mdts) < 0)) || (ts != NULL && (mdio_timestep(gmx->mf, &mdts) < 0))) {
00484     if (mdio_errno() == MDIO_EOF) {
00485         mdio_tsfree(&mdts);
00486         return MOLFILE_EOF;
00487     }
00488     if (mdio_errno() == MDIO_IOERROR) {
00489         mdio_tsfree(&mdts);
00490       // XXX Lame, why does mdio treat IOERROR like EOF?
00491       return MOLFILE_ERROR;
00492     }
00493     fprintf(stderr, "gromacsplugin) Error reading timestep, %s\n",
00494             mdio_errmsg(mdio_errno()));
00495     return MOLFILE_ERROR;
00496   }
00497   if (mdts.natoms != gmx->natoms) {
00498     fprintf(stderr, "gromacsplugin) Timestep in file contains wrong number of atoms\n");
00499     fprintf(stderr, "gromacsplugin) Found %d, expected %d\n", mdts.natoms, gmx->natoms);
00500     mdio_tsfree(&mdts);
00501     return MOLFILE_ERROR;
00502   }
00503 
00504   if (ts) {
00505     if (mdts.pos) 
00506       memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
00507     else 
00508       printf("gromacsplugin) Warning: skipping empty timestep!\n");
00509 
00510     if (mdts.box) {
00511       ts->A = mdts.box->A;
00512       ts->B = mdts.box->B;
00513       ts->C = mdts.box->C;
00514       ts->alpha = mdts.box->alpha;
00515       ts->beta = mdts.box->beta;
00516       ts->gamma = mdts.box->gamma;
00517     }
00518   }
00519   mdio_tsfree(&mdts);
00520   return MOLFILE_SUCCESS;
00521 }
00522 
00523 static void close_trr_read(void *v) {
00524   gmxdata *gmx = (gmxdata *)v;
00525   mdio_close(gmx->mf);
00526   delete gmx;
00527 }
00528 
00529 // open file for writing
00530 static void *open_trr_write(const char *filename, const char *filetype,
00531     int natoms) {
00532 
00533     md_file *mf;
00534     gmxdata *gmx;
00535     int format;
00536 
00537     if (!strcmp(filetype, "trr"))
00538       format = MDFMT_TRR;
00539     else if (!strcmp(filetype, "xtc"))
00540       format = MDFMT_XTC;
00541     else
00542       return NULL;
00543 
00544     mf = mdio_open(filename, format, MDIO_WRITE);
00545     if (!mf) {
00546         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00547                 filename, mdio_errmsg(mdio_errno()));
00548         return NULL;
00549     }
00550     gmx = new gmxdata;
00551     memset(gmx,0,sizeof(gmxdata));
00552     gmx->mf = mf;
00553     gmx->natoms = natoms;
00554     // set some parameters for the output stream:
00555     // start at step 0, convert to big-endian, write single precision.
00556     gmx->step   = 0;
00557     gmx->mf->rev = host_is_little_endian();
00558     gmx->mf->prec = sizeof(float);
00559     return gmx;
00560 }
00561 
00562 // write a trr timestep. the file format has a header with each record
00563 static int write_trr_timestep(void *mydata, const molfile_timestep_t *ts)
00564 {
00565   const float nm=0.1;
00566 
00567   gmxdata *gmx = (gmxdata *)mydata;
00568 
00569   // determine and write header from structure info.
00570   // write trr header. XXX: move this to Gromacs.h ??
00571   if (gmx->mf->fmt == MDFMT_TRR) {
00572     int i;
00573 
00574     if ( put_trx_int(gmx->mf, TRX_MAGIC)            // ID
00575          || put_trx_string(gmx->mf, "GMX_trn_file") // version
00576          || put_trx_int(gmx->mf, 0)                 // ir_size (ignored)
00577          || put_trx_int(gmx->mf, 0)                 // e_size (ignored)
00578          || put_trx_int(gmx->mf, 9*sizeof(float))   // box
00579          || put_trx_int(gmx->mf, 0)                 // vir_size (ignored)
00580          || put_trx_int(gmx->mf, 0)                 // pres_size (ignored)
00581          || put_trx_int(gmx->mf, 0)                 // top_size (ignored)
00582          || put_trx_int(gmx->mf, 0)                 // sym_size (ignored)
00583          || put_trx_int(gmx->mf, 3*sizeof(float)*gmx->natoms) // coordinates
00584          || put_trx_int(gmx->mf, 0)                 // no velocities
00585          || put_trx_int(gmx->mf, 0)                 // no forces
00586          || put_trx_int(gmx->mf, gmx->natoms)       // number of atoms
00587          || put_trx_int(gmx->mf, gmx->step)         // current step number
00588          || put_trx_int(gmx->mf, 0)                 // nre (ignored)
00589          || put_trx_real(gmx->mf, 0.1*gmx->step)    // current time. (dummy value: 0.1)
00590          || put_trx_real(gmx->mf, 0.0))             // current lambda
00591       return MOLFILE_ERROR;
00592 
00593     // set up box according to the VMD unitcell conventions.
00594     // the a-vector is collinear with the x-axis and
00595     // the b-vector is in the xy-plane.
00596     const float sa = sin((double)ts->alpha/180.0*M_PI);
00597     const float ca = cos((double)ts->alpha/180.0*M_PI);
00598     const float cb = cos((double)ts->beta/180.0*M_PI);
00599     const float cg = cos((double)ts->gamma/180.0*M_PI);
00600     const float sg = sin((double)ts->gamma/180.0*M_PI);
00601     float box[9];
00602     box[0] = ts->A;    box[1] = 0.0;      box[2] = 0.0;
00603     box[3] = ts->B*ca; box[4] = ts->B*sa; box[5] = 0.0;
00604     box[6] = ts->C*cb; box[7] = ts->C*(ca - cb*cg)/sg;
00605     box[8] = ts->C*sqrt((double)(1.0 + 2.0*ca*cb*cg
00606                                  - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
00607 
00608     for (i=0; i<9; ++i) {
00609       if (put_trx_real(gmx->mf, box[i]*nm))
00610         return MOLFILE_ERROR;
00611     }
00612 #ifdef TEST_TRR_PLUGIN
00613     fprintf(stderr, "gromacsplugin) box is:\n %f %f %f\n %f %f %f\n %f %f %f\n\n",
00614             box[0], box[1], box[2], box[3], box[4], box[5], box[6], box[7], box[8]);
00615 #endif
00616 
00617     // write coordinates
00618     for (i=0; i<(3*gmx->natoms); ++i) {
00619       if (put_trx_real(gmx->mf, ts->coords[i]*nm))
00620         return MOLFILE_ERROR;
00621     }
00622   } else {
00623     fprintf(stderr, "gromacsplugin) only .trr is supported for writing\n");
00624     return MOLFILE_ERROR;
00625   }
00626 
00627   ++ gmx->step;
00628   return MOLFILE_SUCCESS;
00629   }
00630 
00631 
00632 static void close_trr_write(void *v) {
00633   gmxdata *gmx = (gmxdata *)v;
00634   mdio_close(gmx->mf);
00635   delete gmx;
00636 }
00637 
00638 #define GROMACS_PLUGIN_MAJOR_VERSION 1
00639 #define GROMACS_PLUGIN_MINOR_VERSION 4 
00640 
00641 //
00642 // plugin registration stuff below
00643 //
00644 
00645 static molfile_plugin_t gro_plugin;
00646 static molfile_plugin_t g96_plugin;
00647 static molfile_plugin_t trr_plugin;
00648 static molfile_plugin_t xtc_plugin;
00649 static molfile_plugin_t trj_plugin;
00650 
00651 
00652 VMDPLUGIN_API int VMDPLUGIN_init() {
00653   // GRO plugin init
00654   memset(&gro_plugin, 0, sizeof(molfile_plugin_t));
00655   gro_plugin.abiversion = vmdplugin_ABIVERSION;
00656   gro_plugin.type = MOLFILE_PLUGIN_TYPE;
00657   gro_plugin.name = "gro";
00658   gro_plugin.prettyname = "Gromacs GRO";
00659   gro_plugin.author = "David Norris, Justin Gullingsrud, Magnus Lundborg";
00660   gro_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00661   gro_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00662   gro_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00663   gro_plugin.filename_extension = "gro";
00664   gro_plugin.open_file_read = open_gro_read;
00665   gro_plugin.read_structure = read_gro_structure;
00666   gro_plugin.read_next_timestep = read_gro_timestep;
00667   gro_plugin.close_file_read = close_gro_read;
00668   gro_plugin.open_file_write = open_gro_write;
00669   gro_plugin.write_structure = write_gro_structure;
00670   gro_plugin.write_timestep = write_gro_timestep;
00671   gro_plugin.close_file_write = close_gro_write;
00672   gro_plugin.read_molecule_metadata = read_gro_molecule_metadata;
00673 
00674   // G96 plugin init
00675   memset(&g96_plugin, 0, sizeof(molfile_plugin_t));
00676   g96_plugin.abiversion = vmdplugin_ABIVERSION;
00677   g96_plugin.type = MOLFILE_PLUGIN_TYPE;
00678   g96_plugin.name = "g96";
00679   g96_plugin.prettyname = "Gromacs g96";
00680   g96_plugin.author = "David Norris, Justin Gullingsrud";
00681   g96_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00682   g96_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00683   g96_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00684   g96_plugin.filename_extension = "g96";
00685   g96_plugin.open_file_read = open_g96_read;
00686   g96_plugin.read_structure = read_g96_structure;
00687   g96_plugin.read_next_timestep = read_g96_timestep;
00688   g96_plugin.close_file_read = close_g96_read;
00689 
00690   // TRR plugin
00691   memset(&trr_plugin, 0, sizeof(molfile_plugin_t));
00692   trr_plugin.abiversion = vmdplugin_ABIVERSION;
00693   trr_plugin.type = MOLFILE_PLUGIN_TYPE;
00694   trr_plugin.name = "trr";
00695   trr_plugin.prettyname = "Gromacs TRR Trajectory";
00696   trr_plugin.author = "David Norris, Justin Gullingsrud, Axel Kohlmeyer, Josh Vermaas";
00697   trr_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00698   trr_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00699   trr_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00700   trr_plugin.filename_extension = "trr";
00701   trr_plugin.open_file_read = open_trr_read;
00702   trr_plugin.read_next_timestep = read_trr_timestep;
00703   trr_plugin.close_file_read = close_trr_read;
00704   trr_plugin.open_file_write = open_trr_write;
00705   trr_plugin.write_timestep = write_trr_timestep;
00706   trr_plugin.close_file_write = close_trr_write;
00707 
00708   // XTC plugin 
00709   memset(&xtc_plugin, 0, sizeof(molfile_plugin_t));
00710   xtc_plugin.abiversion = vmdplugin_ABIVERSION;
00711   xtc_plugin.type = MOLFILE_PLUGIN_TYPE;
00712   xtc_plugin.name = "xtc";
00713   xtc_plugin.prettyname = "Gromacs XTC Compressed Trajectory";
00714   xtc_plugin.author = "David Norris, Justin Gullingsrud, Josh Vermaas";
00715   xtc_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00716   xtc_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00717   xtc_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00718   xtc_plugin.filename_extension = "xtc";
00719   xtc_plugin.open_file_read = open_trr_read;
00720   xtc_plugin.read_next_timestep = read_trr_timestep;
00721   xtc_plugin.close_file_read = close_trr_read;
00722 
00723   // TRJ plugin
00724   memset(&trj_plugin, 0, sizeof(molfile_plugin_t));
00725   trj_plugin.abiversion = vmdplugin_ABIVERSION;
00726   trj_plugin.type = MOLFILE_PLUGIN_TYPE;
00727   trj_plugin.name = "trj";
00728   trj_plugin.prettyname = "Gromacs TRJ Trajectory";
00729   trj_plugin.author = "David Norris, Justin Gullingsrud";
00730   trj_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00731   trj_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00732   trj_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00733   trj_plugin.filename_extension = "trj";
00734   trj_plugin.open_file_read = open_trr_read;
00735   trj_plugin.read_next_timestep = read_trr_timestep;
00736   trj_plugin.close_file_read = close_trr_read;
00737 
00738   return 0;
00739 }
00740 
00741 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00742   (*cb)(v, (vmdplugin_t *)&gro_plugin);
00743   (*cb)(v, (vmdplugin_t *)&g96_plugin);
00744   (*cb)(v, (vmdplugin_t *)&trr_plugin);
00745   (*cb)(v, (vmdplugin_t *)&trj_plugin);
00746   (*cb)(v, (vmdplugin_t *)&xtc_plugin);
00747   return 0;
00748 }
00749 
00750 VMDPLUGIN_API int VMDPLUGIN_fini() {
00751   return 0;
00752 }
00753 
00754 
00755 #ifdef TEST_G96_PLUGIN
00756 
00757 int main(int argc, char *argv[]) {
00758   int natoms;
00759 
00760   molfile_timestep_t timestep;
00761   void *v;
00762   int i;
00763 
00764   if (argc < 2) return 1;
00765   while (--argc) {
00766     ++argv;
00767     v = open_g96_read(*argv, "g96", &natoms);
00768     if (!v) {
00769       fprintf(stderr, "open_g96_read failed for file %s\n", *argv);
00770       return 1;
00771     }
00772     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00773     i = 0;
00774     while(!read_g96_timestep(v, natoms, &timestep)) {
00775       ++i;
00776     }
00777     fprintf(stderr, "ended read_g96_timestep on step %d\n", i);
00778     free(timestep.coords);
00779     close_g96_read(v);
00780   }
00781   return 0;
00782 }
00783 
00784 #endif
00785 
00786 #ifdef TEST_TRR_PLUGIN
00787 
00788 int main(int argc, char *argv[]) {
00789   int natoms;
00790 
00791   molfile_timestep_t timestep;
00792   void *v, *w;
00793   int i;
00794 
00795   if (argc != 3) return 1;
00796   v = open_trr_read(argv[1], "trr", &natoms);
00797   if (!v) {
00798     fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
00799     return 1;
00800   }
00801   timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00802   w = open_trr_write(argv[2], "trr", natoms);
00803   if (!w) {
00804     fprintf(stderr, "open_trr_write failed for file %s\n", argv[2]);
00805     return 1;
00806   }
00807 
00808   i = 0;
00809   while(!read_trr_timestep(v, natoms, &timestep)) {
00810     ++i;
00811     if (write_trr_timestep(w, &timestep)) {
00812       fprintf(stderr, "write error\n");
00813       return 1;
00814     }
00815   }
00816 
00817   fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
00818   free(timestep.coords);
00819   close_trr_read(v);
00820   close_trr_write(w);
00821   return 0;
00822 }
00823 
00824 #endif
00825 
00826 #ifdef TEST_TRR_SKIP
00827 
00828 int main(int argc, char *argv[]) {
00829   int natoms;
00830 
00831   molfile_timestep_t timestep;
00832   void *v, *w;
00833   int i;
00834 
00835   if (argc != 2) return 1;
00836   v = open_trr_read(argv[1], "trr", &natoms);
00837   if (!v) {
00838     fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
00839     return 1;
00840   }
00841   timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00842 
00843   i = 0;
00844   printf("natoms: %d\n", natoms);
00845   printf("Frame %d\n", i++);
00846   read_trr_timestep(v, natoms, &timestep);
00847   while(! read_trr_timestep(v, natoms, NULL)){
00848     printf("Frame %d\n", i++);
00849   }
00850 
00851   fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
00852   free(timestep.coords);
00853   close_trr_read(v);
00854   // close_trr_write(w);
00855   return 0;
00856 }
00857 
00858 #endif
00859 
00860 #ifdef TEST_XTCSKIP
00861 
00862 int main(int argc, char *argv[]) {
00863   int natoms;
00864 
00865   molfile_timestep_t timestep;
00866   void *v, *w;
00867   int i;
00868 
00869   if (argc != 2) return 1;
00870   v = open_trr_read(argv[1], "xtc", &natoms);
00871   if (!v) {
00872     fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
00873     return 1;
00874   }
00875   timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00876 
00877   i = 0;
00878   printf("natoms: %d\n", natoms);
00879   printf("Frame %d\n", i++);
00880   read_trr_timestep(v, natoms, &timestep);
00881   while(! read_trr_timestep(v, natoms, NULL)){
00882     printf("Frame %d\n", i++);
00883   }
00884 
00885   fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
00886   free(timestep.coords);
00887   close_trr_read(v);
00888   return 0;
00889 }
00890 
00891 #endif
00892 
00893 

Generated on Fri May 24 03:07:35 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002