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.54 $ $Date: 2017/05/20 05:37:53 $
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,
277 0 : (i+1) % 100000, // GRO format only supports indices up to 99999
278 : // but since GROMACS ignores indices, modular
279 : // arithmetic prevents formatting problems for
280 : // very large structures
281 0 : pos[0] / ANGS_PER_NM, pos[1] / ANGS_PER_NM, pos[2] / ANGS_PER_NM);
282 0 : if(vel)
283 : {
284 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);
285 0 : vel += 3;
286 : }
287 0 : fprintf(gmx->mf->f, "\n");
288 0 : ++atom;
289 0 : pos += 3;
290 : }
291 0 : convert_vmd_box_for_writing(ts, x, y, z);
292 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]);
293 :
294 : return MOLFILE_SUCCESS;
295 : }
296 :
297 0 : static void close_gro_write(void *v) {
298 : gmxdata *gmx = (gmxdata *)v;
299 0 : mdio_close(gmx->mf);
300 0 : free(gmx->atomlist);
301 0 : delete gmx->meta;
302 0 : delete gmx;
303 0 : }
304 :
305 :
306 0 : static void *open_g96_read(const char *filename, const char *,
307 : int *natoms) {
308 :
309 : md_file *mf;
310 : md_header mdh;
311 : char gbuf[MAX_G96_LINE + 1];
312 :
313 0 : mf = mdio_open(filename, MDFMT_G96);
314 0 : if (!mf) {
315 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
316 : filename, mdio_errmsg(mdio_errno()));
317 0 : return NULL;
318 : }
319 :
320 : // read in the header data
321 0 : if (g96_header(mf, mdh.title, MAX_MDIO_TITLE, &mdh.timeval) < 0) {
322 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
323 : filename, mdio_errmsg(mdio_errno()));
324 0 : return NULL;
325 : }
326 :
327 : // First, look for a timestep block
328 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
329 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
330 : filename, mdio_errmsg(mdio_errno()));
331 0 : return NULL;
332 : }
333 0 : if (!strcasecmp(gbuf, "TIMESTEP")) {
334 : // Read in the value line and the END line, and the next
335 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
336 0 : mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0 ||
337 0 : mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
338 0 : fprintf(stderr, "gromacsplugin) Cannot read header from '%s', %s\n",
339 : filename, mdio_errmsg(mdio_errno()));
340 0 : return NULL;
341 : }
342 : }
343 0 : if (strcasecmp(gbuf, "POSITION") && strcasecmp(gbuf, "REFPOSITION")) {
344 0 : fprintf(stderr, "gromacsplugin) No structure information in file %s\n", filename);
345 0 : return NULL;
346 : }
347 0 : *natoms = g96_countatoms(mf);
348 :
349 0 : gmxdata *gmx = new gmxdata;
350 : memset(gmx,0,sizeof(gmxdata));
351 0 : gmx->mf = mf;
352 0 : gmx->natoms = *natoms;
353 0 : return gmx;
354 : }
355 :
356 0 : static int read_g96_structure(void *mydata, int *optflags,
357 : molfile_atom_t *atoms) {
358 :
359 : char gbuf[MAX_G96_LINE + 1];
360 : gmxdata *gmx = (gmxdata *)mydata;
361 : md_atom ma;
362 0 : md_file *mf = gmx->mf;
363 :
364 0 : *optflags = MOLFILE_NOOPTIONS; // no optional data
365 :
366 0 : for (int i = 0; i < gmx->natoms; i++) {
367 0 : molfile_atom_t *atom = atoms+i;
368 0 : if (g96_rec(mf, &ma) < 0) {
369 0 : fprintf(stderr, "gromacsplugin) Error reading atom %d from file, %s\n",
370 : i+1, mdio_errmsg(mdio_errno()));
371 0 : return MOLFILE_ERROR;
372 : }
373 0 : strcpy(atom->name, ma.atomname);
374 0 : strcpy(atom->type, ma.atomname);
375 0 : strcpy(atom->resname, ma.resname);
376 0 : atom->resid = atoi(ma.resid);
377 0 : atom->chain[0] = '\0';
378 0 : atom->segid[0] = '\0';
379 : }
380 :
381 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0) {
382 0 : fprintf(stderr, "gromacsplugin) Warning, error reading END record, %s\n",
383 : mdio_errmsg(mdio_errno()));
384 : }
385 :
386 : // ... another problem: there may or may not be a VELOCITY
387 : // block or a BOX block, so we need to read one line beyond
388 : // the POSITION block to determine this. If neither VEL. nor
389 : // BOX are present we've read a line too far and infringed
390 : // on the next timestep, so we need to keep track of the
391 : // position now for a possible fseek() later to backtrack.
392 0 : long fpos = ftell(mf->f);
393 :
394 : // Now we must read in the velocities and the box, if present
395 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) >= 0) {
396 :
397 : // Is there a velocity block present ?
398 0 : if (!strcasecmp(gbuf, "VELOCITY") || !strcasecmp(gbuf, "VELOCITYRED")) {
399 : // Ignore all the coordinates - VMD doesn't use them
400 : for (;;) {
401 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
402 : return MOLFILE_ERROR;
403 0 : if (!strcasecmp(gbuf, "END")) break;
404 : }
405 :
406 : // Again, record our position because we may need
407 : // to fseek here later if we read too far.
408 0 : fpos = ftell(mf->f);
409 :
410 : // Go ahead and read the next line.
411 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
412 : return MOLFILE_ERROR;
413 : }
414 :
415 : // Is there a box present ?
416 0 : if (!strcasecmp(gbuf, "BOX")) {
417 : // Ignore the box coordinates at this time.
418 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
419 : return MOLFILE_ERROR;
420 0 : if (mdio_readline(mf, gbuf, MAX_G96_LINE + 1) < 0)
421 : return MOLFILE_ERROR;
422 0 : if (strcasecmp(gbuf, "END"))
423 : return MOLFILE_ERROR;
424 : }
425 : else {
426 : // We have read too far, so fseek back to the
427 : // last known safe position so we don't return
428 : // with the file pointer set infringing on the
429 : // next timestep data.
430 0 : fseek(mf->f, fpos, SEEK_SET);
431 : }
432 : }
433 : else {
434 : // Go ahead and rewind for good measure
435 0 : fseek(mf->f, fpos, SEEK_SET);
436 : }
437 0 : rewind(mf->f);
438 : return MOLFILE_SUCCESS;
439 : }
440 :
441 0 : static int read_g96_timestep(void *v, int natoms, molfile_timestep_t *ts) {
442 :
443 : gmxdata *gmx = (gmxdata *)v;
444 : md_ts mdts;
445 : memset(&mdts, 0, sizeof(md_ts));
446 0 : mdts.natoms = natoms;
447 :
448 0 : if (mdio_timestep(gmx->mf, &mdts) < 0)
449 : return MOLFILE_ERROR;
450 0 : if (ts) {
451 0 : memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
452 0 : if (mdts.box) {
453 0 : ts->A = mdts.box->A;
454 0 : ts->B = mdts.box->B;
455 0 : ts->C = mdts.box->C;
456 0 : ts->alpha = mdts.box->alpha;
457 0 : ts->beta = mdts.box->beta;
458 0 : ts->gamma = mdts.box->gamma;
459 : }
460 : }
461 0 : mdio_tsfree(&mdts);
462 0 : return MOLFILE_SUCCESS;
463 : }
464 :
465 0 : static void close_g96_read(void *v) {
466 : gmxdata *gmx = (gmxdata *)v;
467 0 : mdio_close(gmx->mf);
468 0 : delete gmx;
469 0 : }
470 :
471 :
472 : //
473 : // TRR and XTC files
474 : //
475 :
476 209 : static void *open_trr_read(const char *filename, const char *filetype,
477 : int *natoms) {
478 :
479 : md_file *mf;
480 : md_header mdh;
481 : gmxdata *gmx;
482 : int format;
483 :
484 209 : if (!strcmp(filetype, "trr"))
485 : format = MDFMT_TRR;
486 204 : else if (!strcmp(filetype, "trj"))
487 : format = MDFMT_TRJ;
488 204 : else if (!strcmp(filetype, "xtc"))
489 : format = MDFMT_XTC;
490 : else
491 : return NULL;
492 :
493 209 : mf = mdio_open(filename, format);
494 209 : if (!mf) {
495 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
496 : filename, mdio_errmsg(mdio_errno()));
497 0 : return NULL;
498 : }
499 209 : if (mdio_header(mf, &mdh) < 0) {
500 0 : mdio_close(mf);
501 0 : fprintf(stderr, "gromacsplugin) Cannot read header fromm '%s', %s\n",
502 : filename, mdio_errmsg(mdio_errno()));
503 0 : return NULL;
504 : }
505 209 : *natoms = mdh.natoms;
506 209 : gmx = new gmxdata;
507 : memset(gmx,0,sizeof(gmxdata));
508 209 : gmx->mf = mf;
509 209 : gmx->natoms = mdh.natoms;
510 209 : return gmx;
511 : }
512 :
513 18758 : static int read_trr_timestep(void *v, int natoms, molfile_timestep_t *ts) {
514 : gmxdata *gmx = (gmxdata *)v;
515 : md_ts mdts;
516 : memset(&mdts, 0, sizeof(md_ts));
517 18758 : mdts.natoms = natoms;
518 :
519 18758 : if (mdio_timestep(gmx->mf, &mdts) < 0) {
520 209 : if (mdio_errno() == MDIO_EOF || mdio_errno() == MDIO_IOERROR) {
521 : // XXX Lame, why does mdio treat IOERROR like EOF?
522 : return MOLFILE_ERROR;
523 : }
524 0 : fprintf(stderr, "gromacsplugin) Error reading timestep, %s\n",
525 : mdio_errmsg(mdio_errno()));
526 0 : return MOLFILE_ERROR;
527 : }
528 18549 : if (mdts.natoms != gmx->natoms) {
529 0 : fprintf(stderr, "gromacsplugin) Timestep in file contains wrong number of atoms\n");
530 0 : fprintf(stderr, "gromacsplugin) Found %d, expected %d\n", mdts.natoms, gmx->natoms);
531 0 : mdio_tsfree(&mdts);
532 0 : return MOLFILE_ERROR;
533 : }
534 :
535 18549 : if (ts) {
536 18549 : if (mdts.pos)
537 18549 : memcpy(ts->coords, mdts.pos, 3 * sizeof(float) * gmx->natoms);
538 : else
539 : printf("gromacsplugin) Warning: skipping empty timestep!\n");
540 :
541 18549 : if (mdts.box) {
542 18549 : ts->A = mdts.box->A;
543 18549 : ts->B = mdts.box->B;
544 18549 : ts->C = mdts.box->C;
545 18549 : ts->alpha = mdts.box->alpha;
546 18549 : ts->beta = mdts.box->beta;
547 18549 : ts->gamma = mdts.box->gamma;
548 : }
549 : }
550 18549 : mdio_tsfree(&mdts);
551 18549 : return MOLFILE_SUCCESS;
552 : }
553 :
554 209 : static void close_trr_read(void *v) {
555 : gmxdata *gmx = (gmxdata *)v;
556 209 : mdio_close(gmx->mf);
557 209 : delete gmx;
558 209 : }
559 :
560 : // open file for writing
561 0 : static void *open_trr_write(const char *filename, const char *filetype,
562 : int natoms) {
563 :
564 : md_file *mf;
565 : gmxdata *gmx;
566 : int format;
567 :
568 0 : if (!strcmp(filetype, "trr"))
569 : format = MDFMT_TRR;
570 0 : else if (!strcmp(filetype, "xtc"))
571 : format = MDFMT_XTC;
572 : else
573 : return NULL;
574 :
575 0 : mf = mdio_open(filename, format, MDIO_WRITE);
576 0 : if (!mf) {
577 0 : fprintf(stderr, "gromacsplugin) Cannot open file '%s', %s\n",
578 : filename, mdio_errmsg(mdio_errno()));
579 0 : return NULL;
580 : }
581 0 : gmx = new gmxdata;
582 : memset(gmx,0,sizeof(gmxdata));
583 0 : gmx->mf = mf;
584 0 : gmx->natoms = natoms;
585 : // set some parameters for the output stream:
586 : // start at step 0, convert to big-endian, write single precision.
587 : gmx->step = 0;
588 0 : gmx->mf->rev = host_is_little_endian();
589 0 : gmx->mf->prec = sizeof(float);
590 0 : return gmx;
591 : }
592 :
593 : // write a trr timestep. the file format has a header with each record
594 0 : static int write_trr_timestep(void *mydata, const molfile_timestep_t *ts)
595 : {
596 : const float nm=0.1;
597 :
598 : gmxdata *gmx = (gmxdata *)mydata;
599 :
600 : // determine and write header from structure info.
601 : // write trr header. XXX: move this to Gromacs.h ??
602 0 : if (gmx->mf->fmt == MDFMT_TRR) {
603 : int i;
604 :
605 0 : if ( put_trx_int(gmx->mf, TRX_MAGIC) // ID
606 0 : || put_trx_string(gmx->mf, "GMX_trn_file") // version
607 0 : || put_trx_int(gmx->mf, 0) // ir_size (ignored)
608 0 : || put_trx_int(gmx->mf, 0) // e_size (ignored)
609 0 : || put_trx_int(gmx->mf, 9*sizeof(float)) // box
610 0 : || put_trx_int(gmx->mf, 0) // vir_size (ignored)
611 0 : || put_trx_int(gmx->mf, 0) // pres_size (ignored)
612 0 : || put_trx_int(gmx->mf, 0) // top_size (ignored)
613 0 : || put_trx_int(gmx->mf, 0) // sym_size (ignored)
614 0 : || put_trx_int(gmx->mf, 3*sizeof(float)*gmx->natoms) // coordinates
615 0 : || put_trx_int(gmx->mf, 0) // no velocities
616 0 : || put_trx_int(gmx->mf, 0) // no forces
617 0 : || put_trx_int(gmx->mf, gmx->natoms) // number of atoms
618 0 : || put_trx_int(gmx->mf, gmx->step) // current step number
619 0 : || put_trx_int(gmx->mf, 0) // nre (ignored)
620 0 : || put_trx_real(gmx->mf, 0.1*gmx->step) // current time. (dummy value: 0.1)
621 0 : || put_trx_real(gmx->mf, 0.0)) // current lambda
622 0 : return MOLFILE_ERROR;
623 :
624 : // set up box according to the VMD unitcell conventions.
625 : // the a-vector is collinear with the x-axis and
626 : // the b-vector is in the xy-plane.
627 0 : const float sa = sin((double)ts->alpha/180.0*M_PI);
628 0 : const float ca = cos((double)ts->alpha/180.0*M_PI);
629 0 : const float cb = cos((double)ts->beta/180.0*M_PI);
630 0 : const float cg = cos((double)ts->gamma/180.0*M_PI);
631 0 : const float sg = sin((double)ts->gamma/180.0*M_PI);
632 : float box[9];
633 0 : box[0] = ts->A; box[1] = 0.0; box[2] = 0.0;
634 0 : box[3] = ts->B*ca; box[4] = ts->B*sa; box[5] = 0.0;
635 0 : box[6] = ts->C*cb; box[7] = ts->C*(ca - cb*cg)/sg;
636 0 : box[8] = ts->C*sqrt((double)(1.0 + 2.0*ca*cb*cg
637 0 : - ca*ca - cb*cb - cg*cg)/(1.0 - cg*cg));
638 :
639 0 : for (i=0; i<9; ++i) {
640 0 : if (put_trx_real(gmx->mf, box[i]*nm))
641 : return MOLFILE_ERROR;
642 : }
643 : #ifdef TEST_TRR_PLUGIN
644 : fprintf(stderr, "gromacsplugin) box is:\n %f %f %f\n %f %f %f\n %f %f %f\n\n",
645 : box[0], box[1], box[2], box[3], box[4], box[5], box[6], box[7], box[8]);
646 : #endif
647 :
648 : // write coordinates
649 0 : for (i=0; i<(3*gmx->natoms); ++i) {
650 0 : if (put_trx_real(gmx->mf, ts->coords[i]*nm))
651 : return MOLFILE_ERROR;
652 : }
653 : } else {
654 0 : fprintf(stderr, "gromacsplugin) only .trr is supported for writing\n");
655 0 : return MOLFILE_ERROR;
656 : }
657 :
658 0 : ++ gmx->step;
659 0 : return MOLFILE_SUCCESS;
660 : }
661 :
662 :
663 0 : static void close_trr_write(void *v) {
664 : gmxdata *gmx = (gmxdata *)v;
665 0 : mdio_close(gmx->mf);
666 0 : delete gmx;
667 0 : }
668 :
669 : #define GROMACS_PLUGIN_MAJOR_VERSION 1
670 : #define GROMACS_PLUGIN_MINOR_VERSION 3
671 :
672 : //
673 : // plugin registration stuff below
674 : //
675 :
676 : static molfile_plugin_t gro_plugin;
677 : static molfile_plugin_t g96_plugin;
678 : static molfile_plugin_t trr_plugin;
679 : static molfile_plugin_t xtc_plugin;
680 : static molfile_plugin_t trj_plugin;
681 :
682 :
683 10632 : VMDPLUGIN_API int VMDPLUGIN_init() {
684 : // GRO plugin init
685 : memset(&gro_plugin, 0, sizeof(molfile_plugin_t));
686 10632 : gro_plugin.abiversion = vmdplugin_ABIVERSION;
687 10632 : gro_plugin.type = MOLFILE_PLUGIN_TYPE;
688 10632 : gro_plugin.name = "gro";
689 10632 : gro_plugin.prettyname = "Gromacs GRO";
690 10632 : gro_plugin.author = "David Norris, Justin Gullingsrud, Magnus Lundborg";
691 10632 : gro_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
692 10632 : gro_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
693 : gro_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
694 10632 : gro_plugin.filename_extension = "gro";
695 10632 : gro_plugin.open_file_read = open_gro_read;
696 10632 : gro_plugin.read_structure = read_gro_structure;
697 10632 : gro_plugin.read_next_timestep = read_gro_timestep;
698 10632 : gro_plugin.close_file_read = close_gro_read;
699 10632 : gro_plugin.open_file_write = open_gro_write;
700 10632 : gro_plugin.write_structure = write_gro_structure;
701 10632 : gro_plugin.write_timestep = write_gro_timestep;
702 10632 : gro_plugin.close_file_write = close_gro_write;
703 10632 : gro_plugin.read_molecule_metadata = read_gro_molecule_metadata;
704 :
705 : // G96 plugin init
706 : memset(&g96_plugin, 0, sizeof(molfile_plugin_t));
707 10632 : g96_plugin.abiversion = vmdplugin_ABIVERSION;
708 10632 : g96_plugin.type = MOLFILE_PLUGIN_TYPE;
709 10632 : g96_plugin.name = "g96";
710 10632 : g96_plugin.prettyname = "Gromacs g96";
711 10632 : g96_plugin.author = "David Norris, Justin Gullingsrud";
712 10632 : g96_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
713 10632 : g96_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
714 : g96_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
715 10632 : g96_plugin.filename_extension = "g96";
716 10632 : g96_plugin.open_file_read = open_g96_read;
717 10632 : g96_plugin.read_structure = read_g96_structure;
718 10632 : g96_plugin.read_next_timestep = read_g96_timestep;
719 10632 : g96_plugin.close_file_read = close_g96_read;
720 :
721 : // TRR plugin
722 : memset(&trr_plugin, 0, sizeof(molfile_plugin_t));
723 10632 : trr_plugin.abiversion = vmdplugin_ABIVERSION;
724 10632 : trr_plugin.type = MOLFILE_PLUGIN_TYPE;
725 10632 : trr_plugin.name = "trr";
726 10632 : trr_plugin.prettyname = "Gromacs TRR Trajectory";
727 10632 : trr_plugin.author = "David Norris, Justin Gullingsrud, Axel Kohlmeyer";
728 10632 : trr_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
729 10632 : trr_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
730 : trr_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
731 10632 : trr_plugin.filename_extension = "trr";
732 10632 : trr_plugin.open_file_read = open_trr_read;
733 10632 : trr_plugin.read_next_timestep = read_trr_timestep;
734 10632 : trr_plugin.close_file_read = close_trr_read;
735 10632 : trr_plugin.open_file_write = open_trr_write;
736 10632 : trr_plugin.write_timestep = write_trr_timestep;
737 10632 : trr_plugin.close_file_write = close_trr_write;
738 :
739 : // XTC plugin
740 : memset(&xtc_plugin, 0, sizeof(molfile_plugin_t));
741 10632 : xtc_plugin.abiversion = vmdplugin_ABIVERSION;
742 10632 : xtc_plugin.type = MOLFILE_PLUGIN_TYPE;
743 10632 : xtc_plugin.name = "xtc";
744 10632 : xtc_plugin.prettyname = "Gromacs XTC Compressed Trajectory";
745 10632 : xtc_plugin.author = "David Norris, Justin Gullingsrud";
746 10632 : xtc_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
747 10632 : xtc_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
748 : xtc_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
749 10632 : xtc_plugin.filename_extension = "xtc";
750 10632 : xtc_plugin.open_file_read = open_trr_read;
751 10632 : xtc_plugin.read_next_timestep = read_trr_timestep;
752 10632 : xtc_plugin.close_file_read = close_trr_read;
753 :
754 : // TRJ plugin
755 : memset(&trj_plugin, 0, sizeof(molfile_plugin_t));
756 10632 : trj_plugin.abiversion = vmdplugin_ABIVERSION;
757 10632 : trj_plugin.type = MOLFILE_PLUGIN_TYPE;
758 10632 : trj_plugin.name = "trj";
759 10632 : trj_plugin.prettyname = "Gromacs TRJ Trajectory";
760 10632 : trj_plugin.author = "David Norris, Justin Gullingsrud";
761 10632 : trj_plugin.majorv = GROMACS_PLUGIN_MAJOR_VERSION;
762 10632 : trj_plugin.minorv = GROMACS_PLUGIN_MINOR_VERSION;
763 : trj_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
764 10632 : trj_plugin.filename_extension = "trj";
765 10632 : trj_plugin.open_file_read = open_trr_read;
766 10632 : trj_plugin.read_next_timestep = read_trr_timestep;
767 10632 : trj_plugin.close_file_read = close_trr_read;
768 :
769 10632 : return 0;
770 : }
771 :
772 10632 : VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
773 10632 : (*cb)(v, (vmdplugin_t *)&gro_plugin);
774 10632 : (*cb)(v, (vmdplugin_t *)&g96_plugin);
775 10632 : (*cb)(v, (vmdplugin_t *)&trr_plugin);
776 10632 : (*cb)(v, (vmdplugin_t *)&trj_plugin);
777 10632 : (*cb)(v, (vmdplugin_t *)&xtc_plugin);
778 10632 : return 0;
779 : }
780 :
781 0 : VMDPLUGIN_API int VMDPLUGIN_fini() {
782 0 : return 0;
783 : }
784 :
785 : }
786 : }
787 :
788 :
789 : #ifdef TEST_G96_PLUGIN
790 :
791 : int main(int argc, char *argv[]) {
792 : int natoms;
793 :
794 : molfile_timestep_t timestep;
795 : void *v;
796 : int i;
797 :
798 : if (argc < 2) return 1;
799 : while (--argc) {
800 : ++argv;
801 : v = open_g96_read(*argv, "g96", &natoms);
802 : if (!v) {
803 : fprintf(stderr, "open_g96_read failed for file %s\n", *argv);
804 : return 1;
805 : }
806 : timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
807 : i = 0;
808 : while(!read_g96_timestep(v, natoms, ×tep)) {
809 : ++i;
810 : }
811 : fprintf(stderr, "ended read_g96_timestep on step %d\n", i);
812 : free(timestep.coords);
813 : close_g96_read(v);
814 : }
815 : return 0;
816 : }
817 :
818 : #endif
819 :
820 : #ifdef TEST_TRR_PLUGIN
821 :
822 : int main(int argc, char *argv[]) {
823 : int natoms;
824 :
825 : molfile_timestep_t timestep;
826 : void *v, *w;
827 : int i;
828 :
829 : if (argc != 3) return 1;
830 : v = open_trr_read(argv[1], "trr", &natoms);
831 : if (!v) {
832 : fprintf(stderr, "open_trr_read failed for file %s\n", argv[1]);
833 : return 1;
834 : }
835 : timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
836 : w = open_trr_write(argv[2], "trr", natoms);
837 : if (!w) {
838 : fprintf(stderr, "open_trr_write failed for file %s\n", argv[2]);
839 : return 1;
840 : }
841 :
842 : i = 0;
843 : while(!read_trr_timestep(v, natoms, ×tep)) {
844 : ++i;
845 : if (write_trr_timestep(w, ×tep)) {
846 : fprintf(stderr, "write error\n");
847 : return 1;
848 : }
849 : }
850 :
851 : fprintf(stderr, "ended read_trr_timestep on step %d\n", i);
852 : free(timestep.coords);
853 : close_trr_read(v);
854 : close_trr_write(w);
855 : return 0;
856 : }
857 :
858 : #endif
859 :
860 : #endif
|
