00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "largefiles.h"
00046 #include "fastio.h"
00047
00048 #include <stdio.h>
00049 #include <sys/stat.h>
00050 #include <sys/types.h>
00051 #include <stdlib.h>
00052 #include <stddef.h>
00053 #include <string.h>
00054 #include <math.h>
00055 #include <time.h>
00056 #include "endianswap.h"
00057 #include "molfile_plugin.h"
00058
00059 #ifndef M_PI_2
00060 #define M_PI_2 1.57079632679489661922
00061 #endif
00062
00063 #define RECSCALE32BIT 1
00064 #define RECSCALE64BIT 2
00065 #define RECSCALEMAX 2
00066
00067 typedef struct {
00068 fio_fd fd;
00069 int natoms;
00070 int nsets;
00071 int setsread;
00072 int istart;
00073 int nsavc;
00074 double delta;
00075 int nfixed;
00076 float *x, *y, *z;
00077 int *freeind;
00078 float *fixedcoords;
00079 int reverse;
00080 int charmm;
00081 int first;
00082 int with_unitcell;
00083 } dcdhandle;
00084
00085
00086 #define DCD_SUCCESS 0
00087 #define DCD_EOF -1
00088 #define DCD_DNE -2
00089 #define DCD_OPENFAILED -3
00090 #define DCD_BADREAD -4
00091 #define DCD_BADEOF -5
00092 #define DCD_BADFORMAT -6
00093 #define DCD_FILEEXISTS -7
00094 #define DCD_BADMALLOC -8
00095 #define DCD_BADWRITE -9
00096
00097
00098 #define DCD_IS_XPLOR 0x00
00099 #define DCD_IS_CHARMM 0x01
00100 #define DCD_HAS_4DIMS 0x02
00101 #define DCD_HAS_EXTRA_BLOCK 0x04
00102 #define DCD_HAS_64BIT_REC 0x08
00103
00104
00105 #define NFILE_POS 8L
00106 #define NSTEP_POS 20L
00107
00108
00109 #define READ(fd, buf, size) fio_fread(((void *) buf), (size), 1, (fd))
00110
00111
00112 #define WRITE(fd, buf, size) fio_fwrite(((void *) buf), (size), 1, (fd))
00113
00114
00115 #define CHECK_FREAD(X, msg) if (X==-1) { return(DCD_BADREAD); }
00116 #define CHECK_FEOF(X, msg) if (X==0) { return(DCD_BADEOF); }
00117
00118
00119
00120 static void print_dcderror(const char *func, int errcode) {
00121 const char *errstr;
00122
00123 switch (errcode) {
00124 case DCD_EOF: errstr = "end of file"; break;
00125 case DCD_DNE: errstr = "file not found"; break;
00126 case DCD_OPENFAILED: errstr = "file open failed"; break;
00127 case DCD_BADREAD: errstr = "error during read"; break;
00128 case DCD_BADEOF: errstr = "premature end of file"; break;
00129 case DCD_BADFORMAT: errstr = "corruption or unrecognized file structure"; break;
00130 case DCD_FILEEXISTS: errstr = "output file already exists"; break;
00131 case DCD_BADMALLOC: errstr = "memory allocation failed"; break;
00132 case DCD_BADWRITE: errstr = "error during write"; break;
00133 case DCD_SUCCESS:
00134 default:
00135 errstr = "no error";
00136 break;
00137 }
00138 printf("dcdplugin) %s: %s\n", func, errstr);
00139 }
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART,
00157 int *NSAVC, double *DELTA, int *NAMNF,
00158 int **FREEINDEXES, float **fixedcoords, int *reverseEndian,
00159 int *charmm)
00160 {
00161 unsigned int input_integer[2];
00162 int i, ret_val, rec_scale;
00163 char hdrbuf[84];
00164 int NTITLE;
00165 int dcdcordmagic;
00166 int hugefile = 0;
00167 char *corp = (char *) &dcdcordmagic;
00168
00169
00170 corp[0] = 'C';
00171 corp[1] = 'O';
00172 corp[2] = 'R';
00173 corp[3] = 'D';
00174
00175
00176
00177
00178
00179 ret_val = READ(fd, input_integer, 2*sizeof(unsigned int));
00180 CHECK_FREAD(ret_val, "reading first int from dcd file");
00181 CHECK_FEOF(ret_val, "reading first int from dcd file");
00182
00183
00184 if ((input_integer[0]+input_integer[1]) == 84) {
00185 *reverseEndian=0;
00186 rec_scale=RECSCALE64BIT;
00187 printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of native endianness\n");
00188 } else if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
00189 *reverseEndian=0;
00190 rec_scale=RECSCALE32BIT;
00191 printf("dcdplugin) detected standard 32-bit DCD file of native endianness\n");
00192 } else {
00193
00194 swap4_aligned(input_integer, 2);
00195 if ((input_integer[0]+input_integer[1]) == 84) {
00196 *reverseEndian=1;
00197 rec_scale=RECSCALE64BIT;
00198 printf("dcdplugin) detected CHARMM -i8 64-bit DCD file of opposite endianness\n");
00199 } else {
00200 swap4_aligned(&input_integer[1], 1);
00201 if (input_integer[0] == 84 && input_integer[1] == dcdcordmagic) {
00202 *reverseEndian=1;
00203 rec_scale=RECSCALE32BIT;
00204 printf("dcdplugin) detected standard 32-bit DCD file of opposite endianness\n");
00205 } else {
00206
00207 printf("dcdplugin) unrecognized DCD header:\n");
00208 printf("dcdplugin) [0]: %10d [1]: %10d\n", input_integer[0], input_integer[1]);
00209 printf("dcdplugin) [0]: 0x%08x [1]: 0x%08x\n", input_integer[0], input_integer[1]);
00210 return DCD_BADFORMAT;
00211
00212 }
00213 }
00214 }
00215
00216
00217 if (rec_scale == RECSCALE64BIT) {
00218 ret_val = READ(fd, input_integer, sizeof(unsigned int));
00219 if (input_integer[0] != dcdcordmagic) {
00220 printf("dcdplugin) failed to find CORD magic in CHARMM -i8 64-bit DCD file\n");
00221 return DCD_BADFORMAT;
00222 }
00223 }
00224
00225
00226 ret_val = READ(fd, hdrbuf, 80);
00227 CHECK_FREAD(ret_val, "buffering header");
00228 CHECK_FEOF(ret_val, "buffering header");
00229
00230
00231
00232
00233
00234 if (*((int *) (hdrbuf + 76)) != 0) {
00235 (*charmm) = DCD_IS_CHARMM;
00236 if (*((int *) (hdrbuf + 40)) != 0)
00237 (*charmm) |= DCD_HAS_EXTRA_BLOCK;
00238
00239 if (*((int *) (hdrbuf + 44)) == 1)
00240 (*charmm) |= DCD_HAS_4DIMS;
00241
00242 if (rec_scale == RECSCALE64BIT)
00243 (*charmm) |= DCD_HAS_64BIT_REC;
00244
00245 } else {
00246 (*charmm) = DCD_IS_XPLOR;
00247 }
00248
00249 if (*charmm & DCD_IS_CHARMM) {
00250
00251 printf("dcdplugin) CHARMM format DCD file (also NAMD 2.1 and later)\n");
00252 } else {
00253
00254 printf("dcdplugin) X-PLOR format DCD file (also NAMD 2.0 and earlier)\n");
00255 }
00256
00257
00258 (*NSET) = *((int *) (hdrbuf));
00259 if (*reverseEndian) swap4_unaligned(NSET, 1);
00260
00261
00262 (*ISTART) = *((int *) (hdrbuf + 4));
00263 if (*reverseEndian) swap4_unaligned(ISTART, 1);
00264
00265
00266 (*NSAVC) = *((int *) (hdrbuf + 8));
00267 if (*reverseEndian) swap4_unaligned(NSAVC, 1);
00268
00269
00270 (*NAMNF) = *((int *) (hdrbuf + 32));
00271 if (*reverseEndian) swap4_unaligned(NAMNF, 1);
00272
00273
00274
00275 if ((*charmm) & DCD_IS_CHARMM) {
00276 float ftmp;
00277 ftmp = *((float *)(hdrbuf+36));
00278 if (*reverseEndian)
00279 swap4_aligned(&ftmp, 1);
00280
00281 *DELTA = (double)ftmp;
00282 } else {
00283 (*DELTA) = *((double *)(hdrbuf + 36));
00284 if (*reverseEndian) swap8_unaligned(DELTA, 1);
00285 }
00286
00287
00288 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00289 CHECK_FREAD(ret_val, "reading second 84 from dcd file");
00290 CHECK_FEOF(ret_val, "reading second 84 from dcd file");
00291 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00292
00293 if (rec_scale == RECSCALE64BIT) {
00294 if ((input_integer[0]+input_integer[1]) != 84) {
00295 return DCD_BADFORMAT;
00296 }
00297 } else {
00298 if (input_integer[0] != 84) {
00299 return DCD_BADFORMAT;
00300 }
00301 }
00302
00303
00304 input_integer[1] = 0;
00305 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00306 CHECK_FREAD(ret_val, "reading size of title block");
00307 CHECK_FEOF(ret_val, "reading size of title block");
00308 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00309
00310 if ((((input_integer[0]+input_integer[1])-4) % 80) == 0) {
00311
00312 ret_val = READ(fd, &NTITLE, sizeof(int));
00313 CHECK_FREAD(ret_val, "reading NTITLE");
00314 CHECK_FEOF(ret_val, "reading NTITLE");
00315 if (*reverseEndian) swap4_aligned(&NTITLE, 1);
00316
00317 if (NTITLE < 0) {
00318 printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n",
00319 NTITLE, NTITLE);
00320 return DCD_BADFORMAT;
00321 }
00322
00323 if (NTITLE > 1000) {
00324 printf("dcdplugin) WARNING: Bogus NTITLE value: %d (hex: %08x)\n",
00325 NTITLE, NTITLE);
00326 if (NTITLE == 1095062083) {
00327 printf("dcdplugin) WARNING: Broken Vega ZZ 2.4.0 DCD file detected\n");
00328 printf("dcdplugin) Assuming 2 title lines, good luck...\n");
00329 NTITLE = 2;
00330 } else {
00331 printf("dcdplugin) Assuming zero title lines, good luck...\n");
00332 NTITLE = 0;
00333 }
00334 }
00335
00336 for (i=0; i<NTITLE; i++) {
00337 fio_fseek(fd, 80, FIO_SEEK_CUR);
00338 CHECK_FEOF(ret_val, "reading TITLE");
00339 }
00340
00341
00342 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00343 CHECK_FREAD(ret_val, "reading size of title block");
00344 CHECK_FEOF(ret_val, "reading size of title block");
00345 } else {
00346 return DCD_BADFORMAT;
00347 }
00348
00349
00350 input_integer[1] = 0;
00351 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00352
00353 CHECK_FREAD(ret_val, "reading a '4'");
00354 CHECK_FEOF(ret_val, "reading a '4'");
00355 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00356
00357 if ((input_integer[0]+input_integer[1]) != 4) {
00358 return DCD_BADFORMAT;
00359 }
00360
00361
00362 ret_val = READ(fd, N, sizeof(int));
00363 CHECK_FREAD(ret_val, "reading number of atoms");
00364 CHECK_FEOF(ret_val, "reading number of atoms");
00365 if (*reverseEndian) swap4_aligned(N, 1);
00366
00367 if (*N > (1L<<30)) {
00368 hugefile=1;
00369 printf("dcdplugin) ***\n");
00370 printf("dcdplugin) *** Trajectory contains over 2^30 atoms.\n");
00371 printf("dcdplugin) *** Huge file integer wraparound handling enabled.\n");
00372 printf("dcdplugin) ***\n");
00373 }
00374
00375
00376 input_integer[1] = 0;
00377 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00378 CHECK_FREAD(ret_val, "reading a '4'");
00379 CHECK_FEOF(ret_val, "reading a '4'");
00380 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00381
00382 if ((input_integer[0]+input_integer[1]) != 4) {
00383 return DCD_BADFORMAT;
00384 }
00385
00386 *FREEINDEXES = NULL;
00387 *fixedcoords = NULL;
00388 if (*NAMNF != 0) {
00389 (*FREEINDEXES) = (int *) calloc((((ptrdiff_t)(*N))-(*NAMNF)), sizeof(int));
00390 if (*FREEINDEXES == NULL)
00391 return DCD_BADMALLOC;
00392
00393 *fixedcoords = (float *) calloc(((ptrdiff_t)(*N))*4L - (*NAMNF), sizeof(float));
00394 if (*fixedcoords == NULL)
00395 return DCD_BADMALLOC;
00396
00397
00398 input_integer[1]=0;
00399 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00400 CHECK_FREAD(ret_val, "reading size of index array");
00401 CHECK_FEOF(ret_val, "reading size of index array");
00402 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00403
00404
00405
00406 if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
00407 return DCD_BADFORMAT;
00408 }
00409
00410 ret_val = READ(fd, (*FREEINDEXES), ((ptrdiff_t) ((*N)-(*NAMNF)))*sizeof(int));
00411 CHECK_FREAD(ret_val, "reading size of index array");
00412 CHECK_FEOF(ret_val, "reading size of index array");
00413
00414 if (*reverseEndian)
00415 swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF)));
00416
00417 input_integer[1]=0;
00418 ret_val = READ(fd, input_integer, rec_scale*sizeof(int));
00419 CHECK_FREAD(ret_val, "reading size of index array");
00420 CHECK_FEOF(ret_val, "reading size of index array");
00421 if (*reverseEndian) swap4_aligned(input_integer, rec_scale);
00422
00423
00424
00425 if (!hugefile && ((input_integer[0]+input_integer[1]) != ((*N)-(*NAMNF))*4L)) {
00426 return DCD_BADFORMAT;
00427 }
00428 }
00429
00430 return DCD_SUCCESS;
00431 }
00432
00433
00434 static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian,
00435 float *unitcell) {
00436 int i, input_integer[2], rec_scale;
00437
00438 if (charmm & DCD_HAS_64BIT_REC) {
00439 rec_scale = RECSCALE64BIT;
00440 } else {
00441 rec_scale = RECSCALE32BIT;
00442 }
00443
00444 if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
00445
00446 input_integer[1] = 0;
00447 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale)
00448 return DCD_BADREAD;
00449 if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00450 if ((input_integer[0]+input_integer[1]) == 48) {
00451 double tmp[6];
00452 if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD;
00453 if (reverseEndian)
00454 swap8_aligned(tmp, 6);
00455 for (i=0; i<6; i++) unitcell[i] = (float)tmp[i];
00456 } else {
00457
00458 if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
00459 }
00460 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00461 }
00462
00463 return DCD_SUCCESS;
00464 }
00465
00466
00467 static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes,
00468 int reverseEndian, const float *fixedcoords,
00469 float *freeatoms, float *pos, int charmm) {
00470 int i, input_integer[2], rec_scale;
00471
00472 if(charmm & DCD_HAS_64BIT_REC) {
00473 rec_scale=RECSCALE64BIT;
00474 } else {
00475 rec_scale=RECSCALE32BIT;
00476 }
00477
00478
00479 input_integer[1]=0;
00480 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00481 if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00482 if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
00483
00484
00485 if (fio_fread(freeatoms, 4L*num_free, 1, fd) != 1) return DCD_BADREAD;
00486 if (reverseEndian)
00487 swap4_aligned(freeatoms, num_free);
00488
00489
00490 memcpy(pos, fixedcoords, 4L*N);
00491 for (i=0; i<num_free; i++)
00492 pos[indexes[i]-1] = freeatoms[i];
00493
00494
00495 input_integer[1]=0;
00496 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00497 if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00498 if ((input_integer[0]+input_integer[1]) != 4L*num_free) return DCD_BADFORMAT;
00499
00500 return DCD_SUCCESS;
00501 }
00502
00503
00504 static int read_charmm_4dim(fio_fd fd, int charmm, int reverseEndian) {
00505 int input_integer[2], rec_scale;
00506
00507 if (charmm & DCD_HAS_64BIT_REC) {
00508 rec_scale=RECSCALE64BIT;
00509 } else {
00510 rec_scale=RECSCALE32BIT;
00511 }
00512
00513
00514
00515 if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
00516 input_integer[1]=0;
00517 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00518 if (reverseEndian) swap4_aligned(input_integer, rec_scale);
00519 if (fio_fseek(fd, (input_integer[0]+input_integer[1]), FIO_SEEK_CUR)) return DCD_BADREAD;
00520 if (fio_fread(input_integer, sizeof(int), rec_scale, fd) != rec_scale) return DCD_BADREAD;
00521 }
00522
00523 return DCD_SUCCESS;
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541 static int read_dcdstep(fio_fd fd, int N, float *X, float *Y, float *Z,
00542 float *unitcell, int num_fixed,
00543 int first, int *indexes, float *fixedcoords,
00544 int reverseEndian, int charmm) {
00545 int ret_val;
00546 ptrdiff_t rec_scale;
00547 int hugefile = (N > (1L<<30)) ? 1 : 0;
00548 int check_reclen = 1;
00549
00550
00551 if (hugefile || (getenv("VMDDCDNOCHECKRECLEN") != NULL))
00552 check_reclen = 0;
00553
00554 if (charmm & DCD_HAS_64BIT_REC) {
00555 rec_scale=RECSCALE64BIT;
00556 } else {
00557 rec_scale=RECSCALE32BIT;
00558 }
00559
00560 if ((num_fixed==0) || first) {
00561
00562
00563 int tmpbuf[6L*RECSCALEMAX];
00564
00565 fio_iovec iov[7];
00566 int i;
00567
00568
00569
00570
00571
00572
00573
00574
00575 ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
00576 if (ret_val) return ret_val;
00577
00578
00579 iov[0].iov_base = (fio_caddr_t) &tmpbuf[0];
00580 iov[0].iov_len = rec_scale*sizeof(int);
00581
00582 iov[1].iov_base = (fio_caddr_t) X;
00583 iov[1].iov_len = sizeof(float)*N;
00584
00585 iov[2].iov_base = (fio_caddr_t) &tmpbuf[1*rec_scale];
00586 iov[2].iov_len = rec_scale*sizeof(int) * 2L;
00587
00588 iov[3].iov_base = (fio_caddr_t) Y;
00589 iov[3].iov_len = sizeof(float)*N;
00590
00591 iov[4].iov_base = (fio_caddr_t) &tmpbuf[3L*rec_scale];
00592 iov[4].iov_len = rec_scale*sizeof(int) * 2L;
00593
00594 iov[5].iov_base = (fio_caddr_t) Z;
00595 iov[5].iov_len = sizeof(float)*N;
00596
00597 iov[6].iov_base = (fio_caddr_t) &tmpbuf[5L*rec_scale];
00598 iov[6].iov_len = rec_scale*sizeof(int);
00599
00600 #if 1
00601
00602
00603
00604
00605
00606
00607
00608 {
00609 int readcnt = 0;
00610 readcnt = fio_fread(iov[0].iov_base, iov[0].iov_len, 1, fd);
00611 readcnt += fio_fread(iov[1].iov_base, iov[1].iov_len, 1, fd);
00612 readcnt += fio_fread(iov[2].iov_base, iov[2].iov_len, 1, fd);
00613 readcnt += fio_fread(iov[3].iov_base, iov[3].iov_len, 1, fd);
00614 readcnt += fio_fread(iov[4].iov_base, iov[4].iov_len, 1, fd);
00615 readcnt += fio_fread(iov[5].iov_base, iov[5].iov_len, 1, fd);
00616 readcnt += fio_fread(iov[6].iov_base, iov[6].iov_len, 1, fd);
00617
00618
00619 if (readcnt != 7)
00620 return DCD_BADREAD;
00621 }
00622 #else
00623
00624 if (fio_readv(fd, &iov[0], 7) != ((fio_size_t) (rec_scale*6L*sizeof(int) + 3L*N*sizeof(float))))
00625 return DCD_BADREAD;
00626 #endif
00627
00628
00629 if (reverseEndian) {
00630 swap4_aligned(&tmpbuf[0], rec_scale*6L);
00631 swap4_aligned(X, N);
00632 swap4_aligned(Y, N);
00633 swap4_aligned(Z, N);
00634 }
00635
00636
00637
00638 if (check_reclen) {
00639
00640 if (rec_scale == 1) {
00641 for (i=0; i<6; i++) {
00642 if (tmpbuf[i] != sizeof(float)*N) return DCD_BADFORMAT;
00643 }
00644 } else {
00645 for (i=0; i<6; i++) {
00646 if ((tmpbuf[2L*i]+tmpbuf[2L*i+1L]) != sizeof(float)*N) return DCD_BADFORMAT;
00647 }
00648 }
00649 }
00650
00651
00652
00653 if (num_fixed && first) {
00654 memcpy(fixedcoords, X, N*sizeof(float));
00655 memcpy(fixedcoords+N, Y, N*sizeof(float));
00656 memcpy(fixedcoords+2L*N, Z, N*sizeof(float));
00657 }
00658
00659
00660
00661
00662
00663 ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00664 if (ret_val) return ret_val;
00665 } else {
00666
00667
00668 ret_val = read_charmm_extrablock(fd, charmm, reverseEndian, unitcell);
00669 if (ret_val) return ret_val;
00670 ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00671 fixedcoords, fixedcoords+3L*N, X, charmm);
00672 if (ret_val) return ret_val;
00673 ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00674 fixedcoords+N, fixedcoords+3L*N, Y, charmm);
00675 if (ret_val) return ret_val;
00676 ret_val = read_fixed_atoms(fd, N, N-num_fixed, indexes, reverseEndian,
00677 fixedcoords+2*N, fixedcoords+3L*N, Z, charmm);
00678 if (ret_val) return ret_val;
00679 ret_val = read_charmm_4dim(fd, charmm, reverseEndian);
00680 if (ret_val) return ret_val;
00681 }
00682
00683 return DCD_SUCCESS;
00684 }
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698 static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm) {
00699 ptrdiff_t seekoffset = 0;
00700 ptrdiff_t rec_scale;
00701
00702 if (charmm & DCD_HAS_64BIT_REC) {
00703 rec_scale=RECSCALE64BIT;
00704 } else {
00705 rec_scale=RECSCALE32BIT;
00706 }
00707
00708
00709 if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) {
00710 seekoffset += 4L*rec_scale + 48L + 4L*rec_scale;
00711 }
00712
00713
00714 seekoffset += 3L * (2L*rec_scale + natoms - nfixed) * 4L;
00715
00716
00717 if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_4DIMS)) {
00718 seekoffset += (2L*rec_scale + natoms - nfixed) * 4L;
00719 }
00720
00721 if (fio_fseek(fd, seekoffset, FIO_SEEK_CUR)) return DCD_BADEOF;
00722
00723 return DCD_SUCCESS;
00724 }
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737 static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N,
00738 const float *X, const float *Y, const float *Z,
00739 const double *unitcell, int charmm) {
00740 int out_integer;
00741
00742 if (charmm) {
00743
00744 if (unitcell != NULL) {
00745 out_integer = 48;
00746 fio_write_int32(fd, out_integer);
00747 WRITE(fd, unitcell, out_integer);
00748 fio_write_int32(fd, out_integer);
00749 }
00750 }
00751
00752
00753 out_integer = N*4;
00754 fio_write_int32(fd, out_integer);
00755 if (fio_fwrite((void *) X, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00756 fio_write_int32(fd, out_integer);
00757 fio_write_int32(fd, out_integer);
00758 if (fio_fwrite((void *) Y, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00759 fio_write_int32(fd, out_integer);
00760 fio_write_int32(fd, out_integer);
00761 if (fio_fwrite((void *) Z, out_integer, 1, fd) != 1) return DCD_BADWRITE;
00762 fio_write_int32(fd, out_integer);
00763
00764
00765 fio_fseek(fd, NFILE_POS, FIO_SEEK_SET);
00766 fio_write_int32(fd, curframe);
00767 fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET);
00768 fio_write_int32(fd, curstep);
00769 fio_fseek(fd, 0, FIO_SEEK_END);
00770
00771 return DCD_SUCCESS;
00772 }
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784 static int write_dcdheader(fio_fd fd, const char *remarks, int N,
00785 int ISTART, int NSAVC, double DELTA, int with_unitcell,
00786 int charmm) {
00787 int out_integer;
00788 float out_float;
00789 char title_string[200];
00790 time_t cur_time;
00791 struct tm *tmbuf;
00792 char time_str[81];
00793
00794 out_integer = 84;
00795 WRITE(fd, (char *) & out_integer, sizeof(int));
00796 strcpy(title_string, "CORD");
00797 WRITE(fd, title_string, 4);
00798 fio_write_int32(fd, 0);
00799 fio_write_int32(fd, ISTART);
00800 fio_write_int32(fd, NSAVC);
00801 fio_write_int32(fd, 0);
00802 fio_write_int32(fd, 0);
00803 fio_write_int32(fd, 0);
00804 fio_write_int32(fd, 0);
00805 fio_write_int32(fd, 0);
00806 fio_write_int32(fd, 0);
00807 if (charmm) {
00808 out_float = DELTA;
00809 WRITE(fd, (char *) &out_float, sizeof(float));
00810 if (with_unitcell) {
00811 fio_write_int32(fd, 1);
00812 } else {
00813 fio_write_int32(fd, 0);
00814 }
00815 } else {
00816 WRITE(fd, (char *) &DELTA, sizeof(double));
00817 }
00818 fio_write_int32(fd, 0);
00819 fio_write_int32(fd, 0);
00820 fio_write_int32(fd, 0);
00821 fio_write_int32(fd, 0);
00822 fio_write_int32(fd, 0);
00823 fio_write_int32(fd, 0);
00824 fio_write_int32(fd, 0);
00825 fio_write_int32(fd, 0);
00826 if (charmm) {
00827 fio_write_int32(fd, 24);
00828 } else {
00829 fio_write_int32(fd, 0);
00830 }
00831 fio_write_int32(fd, 84);
00832 fio_write_int32(fd, 164);
00833 fio_write_int32(fd, 2);
00834
00835 strncpy(title_string, remarks, 80);
00836 title_string[79] = '\0';
00837 WRITE(fd, title_string, 80);
00838
00839 cur_time=time(NULL);
00840 tmbuf=localtime(&cur_time);
00841 strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf);
00842 WRITE(fd, time_str, 80);
00843
00844 fio_write_int32(fd, 164);
00845 fio_write_int32(fd, 4);
00846 fio_write_int32(fd, N);
00847 fio_write_int32(fd, 4);
00848
00849 return DCD_SUCCESS;
00850 }
00851
00852
00853
00854
00855
00856
00857
00858
00859 static void close_dcd_read(int *indexes, float *fixedcoords) {
00860 free(indexes);
00861 free(fixedcoords);
00862 }
00863
00864
00865 static void *open_dcd_read(const char *path, const char *filetype,
00866 int *natoms) {
00867 dcdhandle *dcd;
00868 fio_fd fd;
00869 int rc;
00870 struct stat stbuf;
00871
00872 if (!path) return NULL;
00873
00874 #if !(defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32))
00875
00876 memset(&stbuf, 0, sizeof(struct stat));
00877 if (stat(path, &stbuf)) {
00878 printf("dcdplugin) Could not access file '%s'.\n", path);
00879 return NULL;
00880 }
00881 #endif
00882
00883 if (fio_open(path, FIO_READ, &fd) < 0) {
00884 printf("dcdplugin) Could not open file '%s' for reading.\n", path);
00885 return NULL;
00886 }
00887
00888 dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
00889 memset(dcd, 0, sizeof(dcdhandle));
00890 dcd->fd = fd;
00891
00892 if ((rc = read_dcdheader(dcd->fd, &dcd->natoms, &dcd->nsets, &dcd->istart,
00893 &dcd->nsavc, &dcd->delta, &dcd->nfixed, &dcd->freeind,
00894 &dcd->fixedcoords, &dcd->reverse, &dcd->charmm))) {
00895 print_dcderror("read_dcdheader", rc);
00896 fio_fclose(dcd->fd);
00897 free(dcd);
00898 return NULL;
00899 }
00900
00901
00902
00903
00904
00905
00906 {
00907 fio_size_t ndims, firstframesize, framesize, extrablocksize;
00908 fio_size_t trjsize, filesize, curpos;
00909 int newnsets;
00910
00911 extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0;
00912 ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3;
00913 firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize;
00914 framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float)
00915 + extrablocksize;
00916
00917
00918
00919
00920
00921
00922 curpos = fio_ftell(dcd->fd);
00923
00924 #if defined(_MSC_VER) && defined(FASTIO_NATIVEWIN32)
00925
00926
00927
00928 fio_fseek(dcd->fd, 0, FIO_SEEK_END);
00929 filesize = fio_ftell(dcd->fd);
00930 fio_fseek(dcd->fd, curpos, FIO_SEEK_SET);
00931 #else
00932 filesize = stbuf.st_size;
00933 #endif
00934 trjsize = filesize - curpos - firstframesize;
00935 if (trjsize < 0) {
00936 printf("dcdplugin) file '%s' appears to contain no timesteps.\n", path);
00937 fio_fclose(dcd->fd);
00938 free(dcd);
00939 return NULL;
00940 }
00941
00942 newnsets = trjsize / framesize + 1;
00943
00944 if (dcd->nsets > 0 && newnsets != dcd->nsets) {
00945 printf("dcdplugin) Warning: DCD header claims %d frames, but \n"
00946 "dcdplugin) file size (%ld) indicates there are actually \n"
00947 "%d frames of size (%ld)\n",
00948 dcd->nsets, trjsize, newnsets, framesize);
00949 }
00950
00951 dcd->nsets = newnsets;
00952 dcd->setsread = 0;
00953 }
00954
00955 dcd->first = 1;
00956 dcd->x = (float *)malloc(dcd->natoms * sizeof(float));
00957 dcd->y = (float *)malloc(dcd->natoms * sizeof(float));
00958 dcd->z = (float *)malloc(dcd->natoms * sizeof(float));
00959 if (!dcd->x || !dcd->y || !dcd->z) {
00960 printf("dcdplugin) Unable to allocate space for %d atoms.\n", dcd->natoms);
00961 if (dcd->x)
00962 free(dcd->x);
00963 if (dcd->y)
00964 free(dcd->y);
00965 if (dcd->z)
00966 free(dcd->z);
00967 fio_fclose(dcd->fd);
00968 free(dcd);
00969 return NULL;
00970 }
00971 *natoms = dcd->natoms;
00972 return dcd;
00973 }
00974
00975
00976 static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00977 dcdhandle *dcd;
00978 int i, j, rc;
00979 float unitcell[6];
00980 unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
00981 unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
00982 dcd = (dcdhandle *)v;
00983
00984
00985 if (dcd->setsread == dcd->nsets) return MOLFILE_EOF;
00986 dcd->setsread++;
00987 if (!ts) {
00988 if (dcd->first && dcd->nfixed) {
00989
00990 rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z,
00991 unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords,
00992 dcd->reverse, dcd->charmm);
00993 dcd->first = 0;
00994 return rc;
00995 }
00996 dcd->first = 0;
00997
00998 return skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm);
00999 }
01000 rc = read_dcdstep(dcd->fd, dcd->natoms, dcd->x, dcd->y, dcd->z, unitcell,
01001 dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords,
01002 dcd->reverse, dcd->charmm);
01003 dcd->first = 0;
01004 if (rc < 0) {
01005 print_dcderror("read_dcdstep", rc);
01006 return MOLFILE_ERROR;
01007 }
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017 {
01018 int natoms = dcd->natoms;
01019 float *nts = ts->coords;
01020 const float *bufx = dcd->x;
01021 const float *bufy = dcd->y;
01022 const float *bufz = dcd->z;
01023
01024 for (i=0, j=0; i<natoms; i++, j+=3) {
01025 nts[j ] = bufx[i];
01026 nts[j + 1] = bufy[i];
01027 nts[j + 2] = bufz[i];
01028 }
01029 }
01030
01031 ts->A = unitcell[0];
01032 ts->B = unitcell[2];
01033 ts->C = unitcell[5];
01034
01035 if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 &&
01036 unitcell[3] >= -1.0 && unitcell[3] <= 1.0 &&
01037 unitcell[4] >= -1.0 && unitcell[4] <= 1.0) {
01038
01039
01040
01041
01042 ts->alpha = 90.0 - asin(unitcell[4]) * 90.0 / M_PI_2;
01043 ts->beta = 90.0 - asin(unitcell[3]) * 90.0 / M_PI_2;
01044 ts->gamma = 90.0 - asin(unitcell[1]) * 90.0 / M_PI_2;
01045 } else {
01046
01047
01048 ts->alpha = unitcell[4];
01049 ts->beta = unitcell[3];
01050 ts->gamma = unitcell[1];
01051 }
01052
01053 return MOLFILE_SUCCESS;
01054 }
01055
01056
01057 static void close_file_read(void *v) {
01058 dcdhandle *dcd = (dcdhandle *)v;
01059 close_dcd_read(dcd->freeind, dcd->fixedcoords);
01060 fio_fclose(dcd->fd);
01061 free(dcd->x);
01062 free(dcd->y);
01063 free(dcd->z);
01064 free(dcd);
01065 }
01066
01067
01068 static void *open_dcd_write(const char *path, const char *filetype,
01069 int natoms) {
01070 dcdhandle *dcd;
01071 fio_fd fd;
01072 int rc;
01073 int istart, nsavc;
01074 double delta;
01075 int with_unitcell;
01076 int charmm;
01077
01078 if (fio_open(path, FIO_WRITE, &fd) < 0) {
01079 printf("dcdplugin) Could not open file '%s' for writing\n", path);
01080 return NULL;
01081 }
01082
01083 dcd = (dcdhandle *)malloc(sizeof(dcdhandle));
01084 memset(dcd, 0, sizeof(dcdhandle));
01085 dcd->fd = fd;
01086
01087 istart = 0;
01088 nsavc = 1;
01089 delta = 1.0;
01090
01091 if (getenv("VMDDCDWRITEXPLORFORMAT") != NULL) {
01092 with_unitcell = 0;
01093 charmm = DCD_IS_XPLOR;
01094 printf("dcdplugin) WARNING: Writing DCD file in X-PLOR format, \n");
01095 printf("dcdplugin) WARNING: unit cell information will be lost!\n");
01096 } else {
01097 with_unitcell = 1;
01098 charmm = DCD_IS_CHARMM;
01099 if (with_unitcell)
01100 charmm |= DCD_HAS_EXTRA_BLOCK;
01101 }
01102
01103 rc = write_dcdheader(dcd->fd, "Created by DCD plugin", natoms,
01104 istart, nsavc, delta, with_unitcell, charmm);
01105
01106 if (rc < 0) {
01107 print_dcderror("write_dcdheader", rc);
01108 fio_fclose(dcd->fd);
01109 free(dcd);
01110 return NULL;
01111 }
01112
01113 dcd->natoms = natoms;
01114 dcd->nsets = 0;
01115 dcd->istart = istart;
01116 dcd->nsavc = nsavc;
01117 dcd->with_unitcell = with_unitcell;
01118 dcd->charmm = charmm;
01119 dcd->x = (float *)malloc(natoms * sizeof(float));
01120 dcd->y = (float *)malloc(natoms * sizeof(float));
01121 dcd->z = (float *)malloc(natoms * sizeof(float));
01122 return dcd;
01123 }
01124
01125
01126 static int write_timestep(void *v, const molfile_timestep_t *ts) {
01127 dcdhandle *dcd = (dcdhandle *)v;
01128 int i, rc, curstep;
01129 float *pos = ts->coords;
01130 double unitcell[6];
01131 unitcell[0] = unitcell[2] = unitcell[5] = 0.0f;
01132 unitcell[1] = unitcell[3] = unitcell[4] = 90.0f;
01133
01134
01135 for (i=0; i<dcd->natoms; i++) {
01136 dcd->x[i] = *(pos++);
01137 dcd->y[i] = *(pos++);
01138 dcd->z[i] = *(pos++);
01139 }
01140 dcd->nsets++;
01141 curstep = dcd->istart + dcd->nsets * dcd->nsavc;
01142
01143 unitcell[0] = ts->A;
01144 unitcell[2] = ts->B;
01145 unitcell[5] = ts->C;
01146 unitcell[1] = sin((M_PI_2 / 90.0) * (90.0 - ts->gamma));
01147 unitcell[3] = sin((M_PI_2 / 90.0) * (90.0 - ts->beta));
01148 unitcell[4] = sin((M_PI_2 / 90.0) * (90.0 - ts->alpha));
01149
01150 rc = write_dcdstep(dcd->fd, dcd->nsets, curstep, dcd->natoms,
01151 dcd->x, dcd->y, dcd->z,
01152 dcd->with_unitcell ? unitcell : NULL,
01153 dcd->charmm);
01154 if (rc < 0) {
01155 print_dcderror("write_dcdstep", rc);
01156 return MOLFILE_ERROR;
01157 }
01158
01159 return MOLFILE_SUCCESS;
01160 }
01161
01162 static void close_file_write(void *v) {
01163 dcdhandle *dcd = (dcdhandle *)v;
01164 fio_fclose(dcd->fd);
01165 free(dcd->x);
01166 free(dcd->y);
01167 free(dcd->z);
01168 free(dcd);
01169 }
01170
01171
01172
01173
01174
01175 static molfile_plugin_t plugin;
01176
01177 VMDPLUGIN_API int VMDPLUGIN_init() {
01178 memset(&plugin, 0, sizeof(molfile_plugin_t));
01179 plugin.abiversion = vmdplugin_ABIVERSION;
01180 plugin.type = MOLFILE_PLUGIN_TYPE;
01181 plugin.name = "dcd";
01182 plugin.prettyname = "CHARMM,NAMD,XPLOR DCD Trajectory";
01183 plugin.author = "Axel Kohlmeyer, Justin Gullingsrud, John Stone";
01184 plugin.majorv = 1;
01185 plugin.minorv = 18;
01186 plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
01187 plugin.filename_extension = "dcd";
01188 plugin.open_file_read = open_dcd_read;
01189 plugin.read_next_timestep = read_next_timestep;
01190 plugin.close_file_read = close_file_read;
01191 plugin.open_file_write = open_dcd_write;
01192 plugin.write_timestep = write_timestep;
01193 plugin.close_file_write = close_file_write;
01194 return VMDPLUGIN_SUCCESS;
01195 }
01196
01197 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01198 (*cb)(v, (vmdplugin_t *)&plugin);
01199 return VMDPLUGIN_SUCCESS;
01200 }
01201
01202 VMDPLUGIN_API int VMDPLUGIN_fini() {
01203 return VMDPLUGIN_SUCCESS;
01204 }
01205
01206
01207 #ifdef TEST_DCDPLUGIN
01208
01209 #include <sys/time.h>
01210
01211
01212 double time_of_day(void) {
01213 #if defined(_MSC_VER)
01214 double t;
01215
01216 t = GetTickCount();
01217 t = t / 1000.0;
01218
01219 return t;
01220 #else
01221 struct timeval tm;
01222 struct timezone tz;
01223
01224 gettimeofday(&tm, &tz);
01225 return((double)(tm.tv_sec) + (double)(tm.tv_usec)/1000000.0);
01226 #endif
01227 }
01228
01229 int main(int argc, char *argv[]) {
01230 molfile_timestep_t timestep;
01231 void *v;
01232 dcdhandle *dcd;
01233 int i, natoms;
01234 float sizeMB =0.0, totalMB = 0.0;
01235 double starttime, endtime, totaltime = 0.0;
01236
01237 while (--argc) {
01238 ++argv;
01239 natoms = 0;
01240 v = open_dcd_read(*argv, "dcd", &natoms);
01241 if (!v) {
01242 fprintf(stderr, "main) open_dcd_read failed for file %s\n", *argv);
01243 return 1;
01244 }
01245 dcd = (dcdhandle *)v;
01246 sizeMB = ((natoms * 3.0) * dcd->nsets * 4.0) / (1024.0 * 1024.0);
01247 totalMB += sizeMB;
01248 printf("main) file: %s\n", *argv);
01249 printf(" %d atoms, %d frames, size: %6.1fMB\n", natoms, dcd->nsets, sizeMB);
01250
01251 starttime = time_of_day();
01252 timestep.coords = (float *)malloc(3*sizeof(float)*natoms);
01253 for (i=0; i<dcd->nsets; i++) {
01254 int rc = read_next_timestep(v, natoms, ×tep);
01255 if (rc) {
01256 fprintf(stderr, "error in read_next_timestep on frame %d\n", i);
01257 return 1;
01258 }
01259 }
01260 endtime = time_of_day();
01261 close_file_read(v);
01262 totaltime += endtime - starttime;
01263 printf(" Time: %5.1f seconds\n", endtime - starttime);
01264 printf(" Speed: %5.1f MB/sec, %5.1f timesteps/sec\n", sizeMB / (endtime - starttime), (dcd->nsets / (endtime - starttime)));
01265 }
01266 printf("Overall Size: %6.1f MB\n", totalMB);
01267 printf("Overall Time: %6.1f seconds\n", totaltime);
01268 printf("Overall Speed: %5.1f MB/sec\n", totalMB / totaltime);
01269 return 0;
01270 }
01271
01272 #endif
01273