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: gromacsplugin.C,v $
50 : * $Author: johns $ $Locker: $ $State: Exp $
51 : * $Revision: 1.52 $ $Date: 2016/11/28 05:01:54 $
52 : *
53 : ***************************************************************************/
54 :
55 : #include "largefiles.h" /* platform dependent 64-bit file I/O defines */
56 :
57 : #include <math.h>
58 : #include <stdio.h>
59 : #include <stdlib.h>
60 : #include <string.h>
61 : #include <ctype.h>
62 : #include "Gromacs.h"
63 : #include "molfile_plugin.h"
64 :
65 : #if defined(_AIX)
66 : #include <strings.h>
67 : #endif
68 :
69 : namespace PLMD{
70 : namespace molfile{
71 :
72 : #ifndef M_PI
73 : #define M_PI 3.14159265358979323846
74 : #endif
75 :
76 : #if defined(WIN32) || defined(WIN64)
77 : #define strcasecmp stricmp
78 : #endif
79 :
80 : typedef struct {
81 : md_file *mf;
82 : int natoms;
83 : int step;
84 : float timeval;
85 : molfile_atom_t *atomlist;
86 : molfile_metadata_t *meta;
87 : } gmxdata;
88 :
89 0 : static void convert_vmd_box_for_writing(const molfile_timestep_t *ts, float *x, float *y, float *z)
90 : {
91 : // const float sa = sin((double)ts->alpha/180.0*M_PI);
92 0 : const float ca = cos((double)ts->alpha/180.0*M_PI);
93 0 : const float cb = cos((double)ts->beta/180.0*M_PI);
94 0 : const float cg = cos((double)ts->gamma/180.0*M_PI);
95 0 : const float sg = sin((double)ts->gamma/180.0*M_PI);
96 :
97 0 : x[0] = ts->A / ANGS_PER_NM;
98 0 : y[0] = 0.0;
99 0 : z[0] = 0.0;
100 0 : x[1] = ts->B*cg / ANGS_PER_NM; // ts->B*ca when writing trr?!
101 0 : y[1] = ts->B*sg / ANGS_PER_NM; // ts->B*sa when writing trr?!
102 0 : z[1] = 0.0;
103 0 : x[2] = ts->C*cb / ANGS_PER_NM;
104 0 : y[2] = (ts->C / ANGS_PER_NM)*(ca - cb*cg)/sg;
105 0 : z[2] = (ts->C / ANGS_PER_NM)*sqrt((double)(1.0 + 2.0*ca*cb*cg
106 0 : - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
107 0 : }
108 :
109 0 : static void *open_gro_read(const char *filename, const char *,
110 : int *natoms) {
111 :
112 : md_file *mf;
113 : md_header mdh;
114 : gmxdata *gmx;
115 :
116 0 : mf = mdio_open(filename, MDFMT_GRO);
117 0 : if (!mf) {
118 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
119 : filename, mdio_errmsg(mdio_errno()));
120 0 : return NULL;
121 : }
122 :
123 : // read in the header data (careful not to rewind!)
124 0 : if (gro_header(mf, mdh.title, MAX_MDIO_TITLE,
125 : &mdh.timeval, &mdh.natoms, 0) < 0) {
126 0 : fprintf(stderr, "gromacsplugin) Cannot read header fromm '%s', %s\n",
127 : filename, mdio_errmsg(mdio_errno()));
128 : // XXX should free the file handle...
129 0 : return NULL;
130 : }
131 0 : *natoms = mdh.natoms;
132 0 : gmx = new gmxdata;
133 : memset(gmx,0,sizeof(gmxdata));
134 0 : gmx->mf = mf;
135 0 : gmx->natoms = mdh.natoms;
136 0 : gmx->meta = new molfile_metadata_t;
137 : memset(gmx->meta,0,sizeof(molfile_metadata_t));
138 0 : strncpy(gmx->meta->title, mdh.title, 80);
139 0 : gmx->timeval = mdh.timeval;
140 0 : return gmx;
141 : }
142 :
143 0 : static int read_gro_structure(void *mydata, int *optflags,
144 : molfile_atom_t *atoms) {
145 :
146 : md_atom ma;
147 : char buf[MAX_GRO_LINE + 1];
148 : gmxdata *gmx = (gmxdata *)mydata;
149 :
150 0 : *optflags = MOLFILE_NOOPTIONS; // no optional data
151 :
152 : // read in each atom and add it into the molecule
153 0 : for (int i = 0; i < gmx->natoms; i++) {
154 0 : molfile_atom_t *atom = atoms+i;
155 0 : if (gro_rec(gmx->mf, &ma) < 0) {
156 0 : fprintf(stderr, "gromacsplugin) Error reading atom %d from file, %s\n",
157 : i+1, mdio_errmsg(mdio_errno()));
158 0 : return MOLFILE_ERROR;
159 : }
160 0 : strcpy(atom->name, ma.atomname);
161 0 : strcpy(atom->type, ma.atomname);
162 0 : strcpy(atom->resname, ma.resname);
163 0 : atom->resid = atoi(ma.resid);
164 0 : atom->chain[0] = '\0';
165 0 : atom->segid[0] = '\0';
166 : }
167 :
168 0 : if (mdio_readline(gmx->mf, buf, MAX_GRO_LINE + 1, 0) < 0) {
169 0 : fprintf(stderr, "gromacsplugin) Warning, error reading box, %s\n",
170 : mdio_errmsg(mdio_errno()));
171 : }
172 :
173 0 : rewind(gmx->mf->f);
174 : return MOLFILE_SUCCESS;
175 : }
176 :
177 0 : static int read_gro_molecule_metadata(void *v, molfile_metadata_t **metadata) {
178 : gmxdata *gmx = (gmxdata *)v;
179 0 : *metadata = gmx->meta;
180 0 : return MOLFILE_SUCCESS;
181 : }
182 :
183 0 : static int read_gro_timestep(void *v, int natoms, molfile_timestep_t *ts) {
184 : gmxdata *gmx = (gmxdata *)v;
185 : md_ts mdts;
186 : memset(&mdts, 0, sizeof(md_ts));
187 0 : mdts.natoms = natoms;
188 :
189 0 : if (mdio_timestep(gmx->mf, &mdts) < 0)
190 : return MOLFILE_ERROR;
191 0 : if (ts) {
192 0 : memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
193 0 : if (mdts.box) {
194 0 : ts->A = mdts.box->A;
195 0 : ts->B = mdts.box->B;
196 0 : ts->C = mdts.box->C;
197 0 : ts->alpha = mdts.box->alpha;
198 0 : ts->beta = mdts.box->beta;
199 0 : ts->gamma = mdts.box->gamma;
200 : }
201 : }
202 0 : mdio_tsfree(&mdts);
203 0 : return MOLFILE_SUCCESS;
204 : }
205 :
206 0 : static void close_gro_read(void *v) {
207 : gmxdata *gmx = (gmxdata *)v;
208 0 : mdio_close(gmx->mf);
209 0 : delete gmx->meta;
210 0 : delete gmx;
211 0 : }
212 :
213 : // open file for writing
214 0 : static void *open_gro_write(const char *filename, const char *filetype,
215 : int natoms) {
216 :
217 : md_file *mf;
218 : gmxdata *gmx;
219 :
220 0 : mf = mdio_open(filename, MDFMT_GRO, MDIO_WRITE);
221 0 : if (!mf) {
222 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
223 : filename, mdio_errmsg(mdio_errno()));
224 0 : return NULL;
225 : }
226 0 : gmx = new gmxdata;
227 : memset(gmx,0,sizeof(gmxdata));
228 0 : gmx->mf = mf;
229 0 : gmx->natoms = natoms;
230 : gmx->step = 0;
231 0 : gmx->meta = new molfile_metadata_t;
232 : memset(gmx->meta,0,sizeof(molfile_metadata_t));
233 : gmx->meta->title[0] = '\0';
234 :
235 0 : return gmx;
236 : }
237 :
238 0 : static int write_gro_structure(void *v, int optflags,
239 : const molfile_atom_t *atoms) {
240 :
241 : gmxdata *gmx = (gmxdata *)v;
242 0 : int natoms = gmx->natoms;
243 0 : gmx->atomlist = (molfile_atom_t *)malloc(natoms*sizeof(molfile_atom_t));
244 : memcpy(gmx->atomlist, atoms, natoms*sizeof(molfile_atom_t));
245 :
246 0 : return MOLFILE_SUCCESS;
247 : }
248 :
249 0 : static int write_gro_timestep(void *v, const molfile_timestep_t *ts) {
250 : gmxdata *gmx = (gmxdata *)v;
251 : const molfile_atom_t *atom;
252 : const float *pos, *vel;
253 : float x[3], y[3], z[3];
254 : int i;
255 :
256 0 : if (gmx->natoms == 0)
257 : return MOLFILE_SUCCESS;
258 :
259 0 : atom = gmx->atomlist;
260 0 : pos = ts->coords;
261 0 : vel = ts->velocities;
262 :
263 : /* The title cannot be written */
264 : /* fprintf(gmx->mf->f, "%s", gmx->meta->title);*/
265 : /* Write a dummy title instead */
266 0 : fprintf(gmx->mf->f, "generated by VMD");
267 : #if vmdplugin_ABIVERSION > 10
268 0 : fprintf(gmx->mf->f, ", t= %f", ts->physical_time);
269 : #endif
270 0 : fprintf(gmx->mf->f, "\n");
271 :
272 0 : fprintf(gmx->mf->f, "%d\n", gmx->natoms);
273 0 : for (i=0; i<gmx->natoms; i++)
274 : {
275 0 : fprintf(gmx->mf->f, "%5d%-5s%5s%5d%8.3f%8.3f%8.3f",
276 0 : atom->resid, atom->resname, atom->name, i+1,
277 0 : pos[0] / ANGS_PER_NM, pos[1] / ANGS_PER_NM, pos[2] / ANGS_PER_NM);
278 0 : if(vel)
279 : {
280 0 : 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);
281 0 : vel += 3;
282 : }
283 0 : fprintf(gmx->mf->f, "\n");
284 0 : ++atom;
285 0 : pos += 3;
286 : }
287 0 : convert_vmd_box_for_writing(ts, x, y, z);
288 0 : 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]);
289 :
290 : return MOLFILE_SUCCESS;
291 : }
292 :
293 0 : static void close_gro_write(void *v) {
294 : gmxdata *gmx = (gmxdata *)v;
295 0 : mdio_close(gmx->mf);
296 0 : free(gmx->atomlist);
297 0 : delete gmx->meta;
298 0 : delete gmx;
299 0 : }
300 :
301 :
302 0 : static void *open_g96_read(const char *filename, const char *,
303 : int *natoms) {
304 :
305 : md_file *mf;
306 : md_header mdh;
307 : char gbuf[MAX_G96_LINE + 1];
308 :
309 0 : mf = mdio_open(filename, MDFMT_G96);
310 0 : if (!mf) {
311 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
312 : filename, mdio_errmsg(mdio_errno()));
313 0 : return NULL;
314 : }
315 :
316 : // read in the header data
317 0 : if (g96_header(mf, mdh.title, MAX_MDIO_TITLE, &mdh.timeval) < 0) {
318 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
319 : filename, mdio_errmsg(mdio_errno()));
320 0 : return NULL;
321 : }
322 :
323 : // First, look for a timestep block
324 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
325 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
326 : filename, mdio_errmsg(mdio_errno()));
327 0 : return NULL;
328 : }
329 0 : if (!strcasecmp(gbuf, "TIMESTEP")) {
330 : // Read in the value line and the END line, and the next
331 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
332 0 : mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
333 0 : mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
334 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
335 : filename, mdio_errmsg(mdio_errno()));
336 0 : return NULL;
337 : }
338 : }
339 0 : if (strcasecmp(gbuf, "POSITION") && strcasecmp(gbuf, "REFPOSITION")) {
340 0 : fprintf(stderr, "gromacsplugin) No structure information in file %s\n", filename);
341 0 : return NULL;
342 : }
343 0 : *natoms = g96_countatoms(mf);
344 :
345 0 : gmxdata *gmx = new gmxdata;
346 : memset(gmx,0,sizeof(gmxdata));
347 0 : gmx->mf = mf;
348 0 : gmx->natoms = *natoms;
349 0 : return gmx;
350 : }
351 :
352 0 : static int read_g96_structure(void *mydata, int *optflags,
353 : molfile_atom_t *atoms) {
354 :
355 : char gbuf[MAX_G96_LINE + 1];
356 : gmxdata *gmx = (gmxdata *)mydata;
357 : md_atom ma;
358 0 : md_file *mf = gmx->mf;
359 :
360 0 : *optflags = MOLFILE_NOOPTIONS; // no optional data
361 :
362 0 : for (int i = 0; i < gmx->natoms; i++) {
363 0 : molfile_atom_t *atom = atoms+i;
364 0 : if (g96_rec(mf, &ma) < 0) {
365 0 : fprintf(stderr, "gromacsplugin) Error reading atom %d from file, %s\n",
366 : i+1, mdio_errmsg(mdio_errno()));
367 0 : return MOLFILE_ERROR;
368 : }
369 0 : strcpy(atom->name, ma.atomname);
370 0 : strcpy(atom->type, ma.atomname);
371 0 : strcpy(atom->resname, ma.resname);
372 0 : atom->resid = atoi(ma.resid);
373 0 : atom->chain[0] = '\0';
374 0 : atom->segid[0] = '\0';
375 : }
376 :
377 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
378 0 : fprintf(stderr, "gromacsplugin) Warning, error reading END record, %s\n",
379 : mdio_errmsg(mdio_errno()));
380 : }
381 :
382 : // ... another problem: there may or may not be a VELOCITY
383 : // block or a BOX block, so we need to read one line beyond
384 : // the POSITION block to determine this. If neither VEL. nor
385 : // BOX are present we've read a line too far and infringed
386 : // on the next timestep, so we need to keep track of the
387 : // position now for a possible fseek() later to backtrack.
388 0 : long fpos = ftell(mf->f);
389 :
390 : // Now we must read in the velocities and the box, if present
391 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) >= 0) {
392 :
393 : // Is there a velocity block present ?
394 0 : if (!strcasecmp(gbuf, "VELOCITY") || !strcasecmp(gbuf, "VELOCITYRED")) {
395 : // Ignore all the coordinates - VMD doesn't use them
396 : for (;;) {
397 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
398 : return MOLFILE_ERROR;
399 0 : if (!strcasecmp(gbuf, "END")) break;
400 : }
401 :
402 : // Again, record our position because we may need
403 : // to fseek here later if we read too far.
404 0 : fpos = ftell(mf->f);
405 :
406 : // Go ahead and read the next line.
407 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
408 : return MOLFILE_ERROR;
409 : }
410 :
411 : // Is there a box present ?
412 0 : if (!strcasecmp(gbuf, "BOX")) {
413 : // Ignore the box coordinates at this time.
414 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
415 : return MOLFILE_ERROR;
416 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
417 : return MOLFILE_ERROR;
418 0 : if (strcasecmp(gbuf, "END"))
419 : return MOLFILE_ERROR;
420 : }
421 : else {
422 : // We have read too far, so fseek back to the
423 : // last known safe position so we don't return
424 : // with the file pointer set infringing on the
425 : // next timestep data.
426 0 : fseek(mf->f, fpos, SEEK_SET);
427 : }
428 : }
429 : else {
430 : // Go ahead and rewind for good measure
431 0 : fseek(mf->f, fpos, SEEK_SET);
432 : }
433 0 : rewind(mf->f);
434 : return MOLFILE_SUCCESS;
435 : }
436 :
437 0 : static int read_g96_timestep(void *v, int natoms, molfile_timestep_t *ts) {
438 :
439 : gmxdata *gmx = (gmxdata *)v;
440 : md_ts mdts;
441 : memset(&mdts, 0, sizeof(md_ts));
442 0 : mdts.natoms = natoms;
443 :
444 0 : if (mdio_timestep(gmx->mf, &mdts) < 0)
445 : return MOLFILE_ERROR;
446 0 : if (ts) {
447 0 : memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
448 0 : if (mdts.box) {
449 0 : ts->A = mdts.box->A;
450 0 : ts->B = mdts.box->B;
451 0 : ts->C = mdts.box->C;
452 0 : ts->alpha = mdts.box->alpha;
453 0 : ts->beta = mdts.box->beta;
454 0 : ts->gamma = mdts.box->gamma;
455 : }
456 : }
457 0 : mdio_tsfree(&mdts);
458 0 : return MOLFILE_SUCCESS;
459 : }
460 :
461 0 : static void close_g96_read(void *v) {
462 : gmxdata *gmx = (gmxdata *)v;
463 0 : mdio_close(gmx->mf);
464 0 : delete gmx;
465 0 : }
466 :
467 :
468 : //
469 : // TRR and XTC files
470 : //
471 :
472 108 : static void *open_trr_read(const char *filename, const char *filetype,
473 : int *natoms) {
474 :
475 : md_file *mf;
476 : md_header mdh;
477 : gmxdata *gmx;
478 : int format;
479 :
480 108 : if (!strcmp(filetype, "trr"))
481 : format = MDFMT_TRR;
482 103 : else if (!strcmp(filetype, "trj"))
483 : format = MDFMT_TRJ;
484 103 : else if (!strcmp(filetype, "xtc"))
485 : format = MDFMT_XTC;
486 : else
487 : return NULL;
488 :
489 108 : mf = mdio_open(filename, format);
490 108 : if (!mf) {
491 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
492 : filename, mdio_errmsg(mdio_errno()));
493 0 : return NULL;
494 : }
495 108 : if (mdio_header(mf, &mdh) < 0) {
496 0 : mdio_close(mf);
497 0 : fprintf(stderr, "gromacsplugin) Cannot read header fromm '%s', %s\n",
498 : filename, mdio_errmsg(mdio_errno()));
499 0 : return NULL;
500 : }
501 108 : *natoms = mdh.natoms;
502 108 : gmx = new gmxdata;
503 : memset(gmx,0,sizeof(gmxdata));
504 108 : gmx->mf = mf;
505 108 : gmx->natoms = mdh.natoms;
506 108 : return gmx;
507 : }
508 :
509 10070 : static int read_trr_timestep(void *v, int natoms, molfile_timestep_t *ts) {
510 : gmxdata *gmx = (gmxdata *)v;
511 : md_ts mdts;
512 : memset(&mdts, 0, sizeof(md_ts));
513 10070 : mdts.natoms = natoms;
514 :
515 10070 : if (mdio_timestep(gmx->mf, &mdts) < 0) {
516 108 : if (mdio_errno() == MDIO_EOF || mdio_errno() == MDIO_IOERROR) {
517 : // XXX Lame, why does mdio treat IOERROR like EOF?
518 : return MOLFILE_ERROR;
519 : }
520 0 : fprintf(stderr, "gromacsplugin) Error reading timestep, %s\n",
521 : mdio_errmsg(mdio_errno()));
522 0 : return MOLFILE_ERROR;
523 : }
524 9962 : if (mdts.natoms != natoms) {
525 0 : fprintf(stderr, "gromacsplugin) Timestep in file contains wrong number of atoms\n");
526 0 : fprintf(stderr, "gromacsplugin) Found %d, expected %d\n", mdts.natoms, natoms);
527 0 : mdio_tsfree(&mdts);
528 0 : return MOLFILE_ERROR;
529 : }
530 :
531 9962 : if (ts) {
532 9962 : memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
533 9962 : if (mdts.box) {
534 9962 : ts->A = mdts.box->A;
535 9962 : ts->B = mdts.box->B;
536 9962 : ts->C = mdts.box->C;
537 9962 : ts->alpha = mdts.box->alpha;
538 9962 : ts->beta = mdts.box->beta;
539 9962 : ts->gamma = mdts.box->gamma;
540 : }
541 : }
542 9962 : mdio_tsfree(&mdts);
543 9962 : return MOLFILE_SUCCESS;
544 : }
545 :
546 108 : static void close_trr_read(void *v) {
547 : gmxdata *gmx = (gmxdata *)v;
548 108 : mdio_close(gmx->mf);
549 108 : delete gmx;
550 108 : }
551 :
552 : // open file for writing
553 0 : static void *open_trr_write(const char *filename, const char *filetype,
554 : int natoms) {
555 :
556 : md_file *mf;
557 : gmxdata *gmx;
558 : int format;
559 :
560 0 : if (!strcmp(filetype, "trr"))
561 : format = MDFMT_TRR;
562 0 : else if (!strcmp(filetype, "xtc"))
563 : format = MDFMT_XTC;
564 : else
565 : return NULL;
566 :
567 0 : mf = mdio_open(filename, format, MDIO_WRITE);
568 0 : if (!mf) {
569 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
570 : filename, mdio_errmsg(mdio_errno()));
571 0 : return NULL;
572 : }
573 0 : gmx = new gmxdata;
574 : memset(gmx,0,sizeof(gmxdata));
575 0 : gmx->mf = mf;
576 0 : gmx->natoms = natoms;
577 : // set some parameters for the output stream:
578 : // start at step 0, convert to big-endian, write single precision.
579 : gmx->step = 0;
580 0 : gmx->mf->rev = host_is_little_endian();
581 0 : gmx->mf->prec = sizeof(float);
582 0 : return gmx;
583 : }
584 :
585 : // write a trr timestep. the file format has a header with each record
586 0 : static int write_trr_timestep(void *mydata, const molfile_timestep_t *ts)
587 : {
588 : const float nm=0.1;
589 :
590 : gmxdata *gmx = (gmxdata *)mydata;
591 :
592 : // determine and write header from structure info.
593 : // write trr header. XXX: move this to Gromacs.h ??
594 0 : if (gmx->mf->fmt == MDFMT_TRR) {
595 : int i;
596 :
597 0 : if ( put_trx_int(gmx->mf, TRX_MAGIC) // ID
598 0 : || put_trx_string(gmx->mf, "GMX_trn_file") // version
599 0 : || put_trx_int(gmx->mf, 0) // ir_size (ignored)
600 0 : || put_trx_int(gmx->mf, 0) // e_size (ignored)
601 0 : || put_trx_int(gmx->mf, 9*sizeof(float)) // box
602 0 : || put_trx_int(gmx->mf, 0) // vir_size (ignored)
603 0 : || put_trx_int(gmx->mf, 0) // pres_size (ignored)
604 0 : || put_trx_int(gmx->mf, 0) // top_size (ignored)
605 0 : || put_trx_int(gmx->mf, 0) // sym_size (ignored)
606 0 : || put_trx_int(gmx->mf, 3*sizeof(float)*gmx->natoms) // coordinates
607 0 : || put_trx_int(gmx->mf, 0) // no velocities
608 0 : || put_trx_int(gmx->mf, 0) // no forces
609 0 : || put_trx_int(gmx->mf, gmx->natoms) // number of atoms
610 0 : || put_trx_int(gmx->mf, gmx->step) // current step number
611 0 : || put_trx_int(gmx->mf, 0) // nre (ignored)
612 0 : || put_trx_real(gmx->mf, 0.1*gmx->step) // current time. (dummy value: 0.1)
613 0 : || put_trx_real(gmx->mf, 0.0)) // current lambda
614 0 : return MOLFILE_ERROR;
615 :
616 : // set up box according to the VMD unitcell conventions.
617 : // the a-vector is collinear with the x-axis and
618 : // the b-vector is in the xy-plane.
619 0 : const float sa = sin((double)ts->alpha/180.0*M_PI);
620 0 : const float ca = cos((double)ts->alpha/180.0*M_PI);
621 0 : const float cb = cos((double)ts->beta/180.0*M_PI);
622 0 : const float cg = cos((double)ts->gamma/180.0*M_PI);
623 0 : const float sg = sin((double)ts->gamma/180.0*M_PI);
624 : float box[9];
625 0 : box[0] = ts->A; box[1] = 0.0; box[2] = 0.0;
626 0 : box[3] = ts->B*ca; box[4] = ts->B*sa; box[5] = 0.0;
627 0 : box[6] = ts->C*cb; box[7] = ts->C*(ca - cb*cg)/sg;
628 0 : box[8] = ts->C*sqrt((double)(1.0 + 2.0*ca*cb*cg
629 0 : - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
630 :
631 0 : for (i=0; i<9; ++i) {
632 0 : if (put_trx_real(gmx->mf, box[i]*nm))
633 : return MOLFILE_ERROR;
634 : }
635 : #ifdef TEST_TRR_PLUGIN
636 : fprintf(stderr, "gromacsplugin) box is:\n %f %f %f\n %f %f %f\n %f %f %f\n\n",
637 : box[0], box[1], box[2], box[3], box[4], box[5], box[6], box[7], box[8]);
638 : #endif
639 :
640 : // write coordinates
641 0 : for (i=0; i<(3*gmx->natoms); ++i) {
642 0 : if (put_trx_real(gmx->mf, ts->coords[i]*nm))
643 : return MOLFILE_ERROR;
644 : }
645 : } else {
646 0 : fprintf(stderr, "gromacsplugin) only .trr is supported for writing\n");
647 0 : return MOLFILE_ERROR;
648 : }
649 :
650 0 : ++ gmx->step;
651 0 : return MOLFILE_SUCCESS;
652 : }
653 :
654 :
655 0 : static void close_trr_write(void *v) {
656 : gmxdata *gmx = (gmxdata *)v;
657 0 : mdio_close(gmx->mf);
658 0 : delete gmx;
659 0 : }
660 :
661 : #define GROMACS_PLUGIN_MAJOR_VERSION 1
662 : #define GROMACS_PLUGIN_MINOR_VERSION 2
663 :
664 : //
665 : // plugin registration stuff below
666 : //
667 :
668 : static molfile_plugin_t gro_plugin;
669 : static molfile_plugin_t g96_plugin;
670 : static molfile_plugin_t trr_plugin;
671 : static molfile_plugin_t xtc_plugin;
672 : static molfile_plugin_t trj_plugin;
673 :
674 :
675 6946 : VMDPLUGIN_API int VMDPLUGIN_init() {
676 : // GRO plugin init
677 : memset(&gro_plugin, 0, sizeof(molfile_plugin_t));
678 6946 : gro_plugin.abiversion = vmdplugin_ABIVERSION;
679 6946 : gro_plugin.type = MOLFILE_PLUGIN_TYPE;
680 6946 : gro_plugin.name = "gro";
681 6946 : gro_plugin.prettyname = "Gromacs GRO";
682 6946 : gro_plugin.author = "David Norris, Justin Gullingsrud, Magnus Lundborg";
683 6946 : gro_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
684 6946 : gro_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
685 : gro_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
686 6946 : gro_plugin.filename_extension = "gro";
687 6946 : gro_plugin.open_file_read = open_gro_read;
688 6946 : gro_plugin.read_structure = read_gro_structure;
689 6946 : gro_plugin.read_next_timestep = read_gro_timestep;
690 6946 : gro_plugin.close_file_read = close_gro_read;
691 6946 : gro_plugin.open_file_write = open_gro_write;
692 6946 : gro_plugin.write_structure = write_gro_structure;
693 6946 : gro_plugin.write_timestep = write_gro_timestep;
694 6946 : gro_plugin.close_file_write = close_gro_write;
695 6946 : gro_plugin.read_molecule_metadata = read_gro_molecule_metadata;
696 :
697 : // G96 plugin init
698 : memset(&g96_plugin, 0, sizeof(molfile_plugin_t));
699 6946 : g96_plugin.abiversion = vmdplugin_ABIVERSION;
700 6946 : g96_plugin.type = MOLFILE_PLUGIN_TYPE;
701 6946 : g96_plugin.name = "g96";
702 6946 : g96_plugin.prettyname = "Gromacs g96";
703 6946 : g96_plugin.author = "David Norris, Justin Gullingsrud";
704 6946 : g96_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
705 6946 : g96_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
706 : g96_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
707 6946 : g96_plugin.filename_extension = "g96";
708 6946 : g96_plugin.open_file_read = open_g96_read;
709 6946 : g96_plugin.read_structure = read_g96_structure;
710 6946 : g96_plugin.read_next_timestep = read_g96_timestep;
711 6946 : g96_plugin.close_file_read = close_g96_read;
712 :
713 : // TRR plugin
714 : memset(&trr_plugin, 0, sizeof(molfile_plugin_t));
715 6946 : trr_plugin.abiversion = vmdplugin_ABIVERSION;
716 6946 : trr_plugin.type = MOLFILE_PLUGIN_TYPE;
717 6946 : trr_plugin.name = "trr";
718 6946 : trr_plugin.prettyname = "Gromacs TRR Trajectory";
719 6946 : trr_plugin.author = "David Norris, Justin Gullingsrud, Axel Kohlmeyer";
720 6946 : trr_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
721 6946 : trr_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
722 : trr_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
723 6946 : trr_plugin.filename_extension = "trr";
724 6946 : trr_plugin.open_file_read = open_trr_read;
725 6946 : trr_plugin.read_next_timestep = read_trr_timestep;
726 6946 : trr_plugin.close_file_read = close_trr_read;
727 6946 : trr_plugin.open_file_write = open_trr_write;
728 6946 : trr_plugin.write_timestep = write_trr_timestep;
729 6946 : trr_plugin.close_file_write = close_trr_write;
730 :
731 : // XTC plugin
732 : memset(&xtc_plugin, 0, sizeof(molfile_plugin_t));
733 6946 : xtc_plugin.abiversion = vmdplugin_ABIVERSION;
734 6946 : xtc_plugin.type = MOLFILE_PLUGIN_TYPE;
735 6946 : xtc_plugin.name = "xtc";
736 6946 : xtc_plugin.prettyname = "Gromacs XTC Compressed Trajectory";
737 6946 : xtc_plugin.author = "David Norris, Justin Gullingsrud";
738 6946 : xtc_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
739 6946 : xtc_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
740 : xtc_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
741 6946 : xtc_plugin.filename_extension = "xtc";
742 6946 : xtc_plugin.open_file_read = open_trr_read;
743 6946 : xtc_plugin.read_next_timestep = read_trr_timestep;
744 6946 : xtc_plugin.close_file_read = close_trr_read;
745 :
746 : // TRJ plugin
747 : memset(&trj_plugin, 0, sizeof(molfile_plugin_t));
748 6946 : trj_plugin.abiversion = vmdplugin_ABIVERSION;
749 6946 : trj_plugin.type = MOLFILE_PLUGIN_TYPE;
750 6946 : trj_plugin.name = "trj";
751 6946 : trj_plugin.prettyname = "Gromacs TRJ Trajectory";
752 6946 : trj_plugin.author = "David Norris, Justin Gullingsrud";
753 6946 : trj_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
754 6946 : trj_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
755 : trj_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
756 6946 : trj_plugin.filename_extension = "trj";
757 6946 : trj_plugin.open_file_read = open_trr_read;
758 6946 : trj_plugin.read_next_timestep = read_trr_timestep;
759 6946 : trj_plugin.close_file_read = close_trr_read;
760 :
761 6946 : return 0;
762 : }
763 :
764 6946 : VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
765 6946 : (*cb)(v, (vmdplugin_t *)&gro_plugin);
766 6946 : (*cb)(v, (vmdplugin_t *)&g96_plugin);
767 6946 : (*cb)(v, (vmdplugin_t *)&trr_plugin);
768 6946 : (*cb)(v, (vmdplugin_t *)&trj_plugin);
769 6946 : (*cb)(v, (vmdplugin_t *)&xtc_plugin);
770 6946 : return 0;
771 : }
772 :
773 0 : VMDPLUGIN_API int VMDPLUGIN_fini() {
774 0 : return 0;
775 : }
776 :
777 : }
778 : }
779 :
780 :
781 : #ifdef TEST_G96_PLUGIN
782 :
783 : int main(int argc, char *argv[]) {
784 : int natoms;
785 :
786 : molfile_timestep_t timestep;
787 : void *v;
788 : int i;
789 :
790 : if (argc < 2) return 1;
791 : while (--argc) {
792 : ++argv;
793 : v = open_g96_read(*argv, "g96", &natoms);
794 : if (!v) {
795 : fprintf(stderr, "open_g96_read failed for file %s\n", *argv);
796 : return 1;
797 : }
798 : timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
799 : i = 0;
800 : while(!read_g96_timestep(v, natoms, ×tep)) {
801 : ++i;
802 : }
803 : fprintf(stderr, "ended read_g96_timestep on step %d\n", i);
804 : free(timestep.coords);
805 : close_g96_read(v);
806 : }
807 : return 0;
808 : }
809 :
810 : #endif
811 :
812 : #ifdef TEST_TRR_PLUGIN
813 :
814 : int main(int argc, char *argv[]) {
815 : int natoms;
816 :
817 : molfile_timestep_t timestep;
818 : void *v, *w;
819 : int i;
820 :
821 : if (argc != 3) return 1;
822 : v = open_trr_read(argv[1], "trr", &natoms);
823 : if (!v) {
824 : fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
825 : return 1;
826 : }
827 : timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
828 : w = open_trr_write(argv[2], "trr", natoms);
829 : if (!w) {
830 : fprintf(stderr, "open_trr_write failed for file %s\n", argv[2]);
831 : return 1;
832 : }
833 :
834 : i = 0;
835 : while(!read_trr_timestep(v, natoms, ×tep)) {
836 : ++i;
837 : if (write_trr_timestep(w, ×tep)) {
838 : fprintf(stderr, "write error\n");
839 : return 1;
840 : }
841 : }
842 :
843 : fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
844 : free(timestep.coords);
845 : close_trr_read(v);
846 : close_trr_write(w);
847 : return 0;
848 : }
849 :
850 : #endif
851 :
852 : #endif
|
