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: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.54 $       $Date: 2017/05/20 05:37:53 $
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   if (mdio_timestep(gmx->mf, &mdts) < 0) {
00480     if (mdio_errno() == MDIO_EOF || mdio_errno() == MDIO_IOERROR) {
00481       // XXX Lame, why does mdio treat IOERROR like EOF?
00482       return MOLFILE_ERROR;
00483     }
00484     fprintf(stderr, "gromacsplugin) Error reading timestep, %s\n",
00485             mdio_errmsg(mdio_errno()));
00486     return MOLFILE_ERROR;
00487   }
00488   if (mdts.natoms != gmx->natoms) {
00489     fprintf(stderr, "gromacsplugin) Timestep in file contains wrong number of atoms\n");
00490     fprintf(stderr, "gromacsplugin) Found %d, expected %d\n", mdts.natoms, gmx->natoms);
00491     mdio_tsfree(&mdts);
00492     return MOLFILE_ERROR;
00493   }
00494 
00495   if (ts) {
00496     if (mdts.pos) 
00497       memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
00498     else 
00499       printf("gromacsplugin) Warning: skipping empty timestep!\n");
00500 
00501     if (mdts.box) {
00502       ts->A = mdts.box->A;
00503       ts->B = mdts.box->B;
00504       ts->C = mdts.box->C;
00505       ts->alpha = mdts.box->alpha;
00506       ts->beta = mdts.box->beta;
00507       ts->gamma = mdts.box->gamma;
00508     }
00509   }
00510   mdio_tsfree(&mdts);
00511   return MOLFILE_SUCCESS;
00512 }
00513 
00514 static void close_trr_read(void *v) {
00515   gmxdata *gmx = (gmxdata *)v;
00516   mdio_close(gmx->mf);
00517   delete gmx;
00518 }
00519 
00520 // open file for writing
00521 static void *open_trr_write(const char *filename, const char *filetype,
00522     int natoms) {
00523 
00524     md_file *mf;
00525     gmxdata *gmx;
00526     int format;
00527 
00528     if (!strcmp(filetype, "trr"))
00529       format = MDFMT_TRR;
00530     else if (!strcmp(filetype, "xtc"))
00531       format = MDFMT_XTC;
00532     else
00533       return NULL;
00534 
00535     mf = mdio_open(filename, format, MDIO_WRITE);
00536     if (!mf) {
00537         fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
00538                 filename, mdio_errmsg(mdio_errno()));
00539         return NULL;
00540     }
00541     gmx = new gmxdata;
00542     memset(gmx,0,sizeof(gmxdata));
00543     gmx->mf = mf;
00544     gmx->natoms = natoms;
00545     // set some parameters for the output stream:
00546     // start at step 0, convert to big-endian, write single precision.
00547     gmx->step   = 0;
00548     gmx->mf->rev = host_is_little_endian();
00549     gmx->mf->prec = sizeof(float);
00550     return gmx;
00551 }
00552 
00553 // write a trr timestep. the file format has a header with each record
00554 static int write_trr_timestep(void *mydata, const molfile_timestep_t *ts)
00555 {
00556   const float nm=0.1;
00557 
00558   gmxdata *gmx = (gmxdata *)mydata;
00559 
00560   // determine and write header from structure info.
00561   // write trr header. XXX: move this to Gromacs.h ??
00562   if (gmx->mf->fmt == MDFMT_TRR) {
00563     int i;
00564 
00565     if ( put_trx_int(gmx->mf, TRX_MAGIC)            // ID
00566          || put_trx_string(gmx->mf, "GMX_trn_file") // version
00567          || put_trx_int(gmx->mf, 0)                 // ir_size (ignored)
00568          || put_trx_int(gmx->mf, 0)                 // e_size (ignored)
00569          || put_trx_int(gmx->mf, 9*sizeof(float))   // box
00570          || put_trx_int(gmx->mf, 0)                 // vir_size (ignored)
00571          || put_trx_int(gmx->mf, 0)                 // pres_size (ignored)
00572          || put_trx_int(gmx->mf, 0)                 // top_size (ignored)
00573          || put_trx_int(gmx->mf, 0)                 // sym_size (ignored)
00574          || put_trx_int(gmx->mf, 3*sizeof(float)*gmx->natoms) // coordinates
00575          || put_trx_int(gmx->mf, 0)                 // no velocities
00576          || put_trx_int(gmx->mf, 0)                 // no forces
00577          || put_trx_int(gmx->mf, gmx->natoms)       // number of atoms
00578          || put_trx_int(gmx->mf, gmx->step)         // current step number
00579          || put_trx_int(gmx->mf, 0)                 // nre (ignored)
00580          || put_trx_real(gmx->mf, 0.1*gmx->step)    // current time. (dummy value: 0.1)
00581          || put_trx_real(gmx->mf, 0.0))             // current lambda
00582       return MOLFILE_ERROR;
00583 
00584     // set up box according to the VMD unitcell conventions.
00585     // the a-vector is collinear with the x-axis and
00586     // the b-vector is in the xy-plane.
00587     const float sa = sin((double)ts->alpha/180.0*M_PI);
00588     const float ca = cos((double)ts->alpha/180.0*M_PI);
00589     const float cb = cos((double)ts->beta/180.0*M_PI);
00590     const float cg = cos((double)ts->gamma/180.0*M_PI);
00591     const float sg = sin((double)ts->gamma/180.0*M_PI);
00592     float box[9];
00593     box[0] = ts->A;    box[1] = 0.0;      box[2] = 0.0;
00594     box[3] = ts->B*ca; box[4] = ts->B*sa; box[5] = 0.0;
00595     box[6] = ts->C*cb; box[7] = ts->C*(ca - cb*cg)/sg;
00596     box[8] = ts->C*sqrt((double)(1.0 + 2.0*ca*cb*cg
00597                                  - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
00598 
00599     for (i=0; i<9; ++i) {
00600       if (put_trx_real(gmx->mf, box[i]*nm))
00601         return MOLFILE_ERROR;
00602     }
00603 #ifdef TEST_TRR_PLUGIN
00604     fprintf(stderr, "gromacsplugin) box is:\n %f %f %f\n %f %f %f\n %f %f %f\n\n",
00605             box[0], box[1], box[2], box[3], box[4], box[5], box[6], box[7], box[8]);
00606 #endif
00607 
00608     // write coordinates
00609     for (i=0; i<(3*gmx->natoms); ++i) {
00610       if (put_trx_real(gmx->mf, ts->coords[i]*nm))
00611         return MOLFILE_ERROR;
00612     }
00613   } else {
00614     fprintf(stderr, "gromacsplugin) only .trr is supported for writing\n");
00615     return MOLFILE_ERROR;
00616   }
00617 
00618   ++ gmx->step;
00619   return MOLFILE_SUCCESS;
00620   }
00621 
00622 
00623 static void close_trr_write(void *v) {
00624   gmxdata *gmx = (gmxdata *)v;
00625   mdio_close(gmx->mf);
00626   delete gmx;
00627 }
00628 
00629 #define GROMACS_PLUGIN_MAJOR_VERSION 1
00630 #define GROMACS_PLUGIN_MINOR_VERSION 3 
00631 
00632 //
00633 // plugin registration stuff below
00634 //
00635 
00636 static molfile_plugin_t gro_plugin;
00637 static molfile_plugin_t g96_plugin;
00638 static molfile_plugin_t trr_plugin;
00639 static molfile_plugin_t xtc_plugin;
00640 static molfile_plugin_t trj_plugin;
00641 
00642 
00643 VMDPLUGIN_API int VMDPLUGIN_init() {
00644   // GRO plugin init
00645   memset(&gro_plugin, 0, sizeof(molfile_plugin_t));
00646   gro_plugin.abiversion = vmdplugin_ABIVERSION;
00647   gro_plugin.type = MOLFILE_PLUGIN_TYPE;
00648   gro_plugin.name = "gro";
00649   gro_plugin.prettyname = "Gromacs GRO";
00650   gro_plugin.author = "David Norris, Justin Gullingsrud, Magnus Lundborg";
00651   gro_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00652   gro_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00653   gro_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00654   gro_plugin.filename_extension = "gro";
00655   gro_plugin.open_file_read = open_gro_read;
00656   gro_plugin.read_structure = read_gro_structure;
00657   gro_plugin.read_next_timestep = read_gro_timestep;
00658   gro_plugin.close_file_read = close_gro_read;
00659   gro_plugin.open_file_write = open_gro_write;
00660   gro_plugin.write_structure = write_gro_structure;
00661   gro_plugin.write_timestep = write_gro_timestep;
00662   gro_plugin.close_file_write = close_gro_write;
00663   gro_plugin.read_molecule_metadata = read_gro_molecule_metadata;
00664 
00665   // G96 plugin init
00666   memset(&g96_plugin, 0, sizeof(molfile_plugin_t));
00667   g96_plugin.abiversion = vmdplugin_ABIVERSION;
00668   g96_plugin.type = MOLFILE_PLUGIN_TYPE;
00669   g96_plugin.name = "g96";
00670   g96_plugin.prettyname = "Gromacs g96";
00671   g96_plugin.author = "David Norris, Justin Gullingsrud";
00672   g96_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00673   g96_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00674   g96_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00675   g96_plugin.filename_extension = "g96";
00676   g96_plugin.open_file_read = open_g96_read;
00677   g96_plugin.read_structure = read_g96_structure;
00678   g96_plugin.read_next_timestep = read_g96_timestep;
00679   g96_plugin.close_file_read = close_g96_read;
00680 
00681   // TRR plugin
00682   memset(&trr_plugin, 0, sizeof(molfile_plugin_t));
00683   trr_plugin.abiversion = vmdplugin_ABIVERSION;
00684   trr_plugin.type = MOLFILE_PLUGIN_TYPE;
00685   trr_plugin.name = "trr";
00686   trr_plugin.prettyname = "Gromacs TRR Trajectory";
00687   trr_plugin.author = "David Norris, Justin Gullingsrud, Axel Kohlmeyer";
00688   trr_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00689   trr_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00690   trr_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00691   trr_plugin.filename_extension = "trr";
00692   trr_plugin.open_file_read = open_trr_read;
00693   trr_plugin.read_next_timestep = read_trr_timestep;
00694   trr_plugin.close_file_read = close_trr_read;
00695   trr_plugin.open_file_write = open_trr_write;
00696   trr_plugin.write_timestep = write_trr_timestep;
00697   trr_plugin.close_file_write = close_trr_write;
00698 
00699   // XTC plugin 
00700   memset(&xtc_plugin, 0, sizeof(molfile_plugin_t));
00701   xtc_plugin.abiversion = vmdplugin_ABIVERSION;
00702   xtc_plugin.type = MOLFILE_PLUGIN_TYPE;
00703   xtc_plugin.name = "xtc";
00704   xtc_plugin.prettyname = "Gromacs XTC Compressed Trajectory";
00705   xtc_plugin.author = "David Norris, Justin Gullingsrud";
00706   xtc_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00707   xtc_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00708   xtc_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00709   xtc_plugin.filename_extension = "xtc";
00710   xtc_plugin.open_file_read = open_trr_read;
00711   xtc_plugin.read_next_timestep = read_trr_timestep;
00712   xtc_plugin.close_file_read = close_trr_read;
00713 
00714   // TRJ plugin
00715   memset(&trj_plugin, 0, sizeof(molfile_plugin_t));
00716   trj_plugin.abiversion = vmdplugin_ABIVERSION;
00717   trj_plugin.type = MOLFILE_PLUGIN_TYPE;
00718   trj_plugin.name = "trj";
00719   trj_plugin.prettyname = "Gromacs TRJ Trajectory";
00720   trj_plugin.author = "David Norris, Justin Gullingsrud";
00721   trj_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
00722   trj_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
00723   trj_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
00724   trj_plugin.filename_extension = "trj";
00725   trj_plugin.open_file_read = open_trr_read;
00726   trj_plugin.read_next_timestep = read_trr_timestep;
00727   trj_plugin.close_file_read = close_trr_read;
00728 
00729   return 0;
00730 }
00731 
00732 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
00733   (*cb)(v, (vmdplugin_t *)&gro_plugin);
00734   (*cb)(v, (vmdplugin_t *)&g96_plugin);
00735   (*cb)(v, (vmdplugin_t *)&trr_plugin);
00736   (*cb)(v, (vmdplugin_t *)&trj_plugin);
00737   (*cb)(v, (vmdplugin_t *)&xtc_plugin);
00738   return 0;
00739 }
00740 
00741 VMDPLUGIN_API int VMDPLUGIN_fini() {
00742   return 0;
00743 }
00744 
00745 
00746 #ifdef TEST_G96_PLUGIN
00747 
00748 int main(int argc, char *argv[]) {
00749   int natoms;
00750 
00751   molfile_timestep_t timestep;
00752   void *v;
00753   int i;
00754 
00755   if (argc < 2) return 1;
00756   while (--argc) {
00757     ++argv;
00758     v = open_g96_read(*argv, "g96", &natoms);
00759     if (!v) {
00760       fprintf(stderr, "open_g96_read failed for file %s\n", *argv);
00761       return 1;
00762     }
00763     timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00764     i = 0;
00765     while(!read_g96_timestep(v, natoms, &timestep)) {
00766       ++i;
00767     }
00768     fprintf(stderr, "ended read_g96_timestep on step %d\n", i);
00769     free(timestep.coords);
00770     close_g96_read(v);
00771   }
00772   return 0;
00773 }
00774 
00775 #endif
00776 
00777 #ifdef TEST_TRR_PLUGIN
00778 
00779 int main(int argc, char *argv[]) {
00780   int natoms;
00781 
00782   molfile_timestep_t timestep;
00783   void *v, *w;
00784   int i;
00785 
00786   if (argc != 3) return 1;
00787   v = open_trr_read(argv[1], "trr", &natoms);
00788   if (!v) {
00789     fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
00790     return 1;
00791   }
00792   timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
00793   w = open_trr_write(argv[2], "trr", natoms);
00794   if (!w) {
00795     fprintf(stderr, "open_trr_write failed for file %s\n", argv[2]);
00796     return 1;
00797   }
00798 
00799   i = 0;
00800   while(!read_trr_timestep(v, natoms, &timestep)) {
00801     ++i;
00802     if (write_trr_timestep(w, &timestep)) {
00803       fprintf(stderr, "write error\n");
00804       return 1;
00805     }
00806   }
00807 
00808   fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
00809   free(timestep.coords);
00810   close_trr_read(v);
00811   close_trr_write(w);
00812   return 0;
00813 }
00814 
00815 #endif
00816 

Generated on Wed Nov 11 03:06:25 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002