1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U)
17 open-source XDR routines from
18 Frans van Hoesel (https://www.rug.nl/staff/f.h.j.van.hoesel/)
19 are included in this file
20 Axel Kohlmeyer (Temple U)
21 port to platforms without XDR support
22 added support for unwrapped trajectories
23 support for groups
24 ------------------------------------------------------------------------- */
25
26 #include "dump_xtc.h"
27 #include <cmath>
28
29 #include <cstring>
30 #include <climits>
31 #include "domain.h"
32 #include "atom.h"
33 #include "update.h"
34 #include "group.h"
35 #include "output.h"
36 #include "force.h"
37 #include "comm.h"
38 #include "memory.h"
39 #include "error.h"
40
41 using namespace LAMMPS_NS;
42
43 #define EPS 1e-5
44 #define XTC_MAGIC 1995
45
46 #define MYMIN(a,b) ((a) < (b) ? (a) : (b))
47 #define MYMAX(a,b) ((a) > (b) ? (a) : (b))
48
49 int xdropen(XDR *, const char *, const char *);
50 int xdrclose(XDR *);
51 void xdrfreebuf();
52 int xdr3dfcoord(XDR *, float *, int *, float *);
53
54 /* ---------------------------------------------------------------------- */
55
DumpXTC(LAMMPS * lmp,int narg,char ** arg)56 DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg),
57 coords(nullptr)
58 {
59 if (narg != 5) error->all(FLERR,"Illegal dump xtc command");
60 if (binary || compressed || multifile || multiproc)
61 error->all(FLERR,"Invalid dump xtc filename");
62
63 size_one = 3;
64 sort_flag = 1;
65 sortcol = 0;
66 format_default = nullptr;
67 flush_flag = 0;
68 unwrap_flag = 0;
69 precision = 1000.0;
70
71 // allocate global array for atom coords
72
73 bigint n = group->count(igroup);
74 if (n > static_cast<bigint>(MAXSMALLINT/3/sizeof(float)))
75 error->all(FLERR,"Too many atoms for dump xtc");
76 natoms = static_cast<int> (n);
77
78 memory->create(coords,3*natoms,"dump:coords");
79
80 // sfactor = conversion of coords to XTC units
81 // tfactor = conversion of simulation time to XTC units
82 // GROMACS standard is nanometers and picoseconds
83
84 sfactor = 0.1 / force->angstrom;
85 tfactor = 0.001 / force->femtosecond;
86
87 // in reduced units we do not scale anything
88 if (strcmp(update->unit_style,"lj") == 0) {
89 sfactor = tfactor = 1.0;
90 if (comm->me == 0)
91 error->warning(FLERR,"No automatic unit conversion to XTC file "
92 "format conventions possible for units lj");
93 }
94
95 DumpXTC::openfile();
96 nevery_save = 0;
97 ntotal = 0;
98 }
99
100 /* ---------------------------------------------------------------------- */
101
~DumpXTC()102 DumpXTC::~DumpXTC()
103 {
104 memory->destroy(coords);
105
106 if (me == 0) {
107 xdrclose(&xd);
108 xdrfreebuf();
109 }
110 }
111
112 /* ---------------------------------------------------------------------- */
113
init_style()114 void DumpXTC::init_style()
115 {
116 if (sort_flag == 0 || sortcol != 0)
117 error->all(FLERR,"Dump xtc requires sorting by atom ID");
118
119 // check that flush_flag is not set since dump::write() will use it
120
121 if (flush_flag) error->all(FLERR,"Cannot set dump_modify flush for dump xtc");
122
123 // check that dump frequency has not changed and is not a variable
124
125 int idump;
126 for (idump = 0; idump < output->ndump; idump++)
127 if (strcmp(id,output->dump[idump]->id) == 0) break;
128 if (output->every_dump[idump] == 0)
129 error->all(FLERR,"Cannot use variable every setting for dump xtc");
130
131 if (nevery_save == 0) nevery_save = output->every_dump[idump];
132 else if (nevery_save != output->every_dump[idump])
133 error->all(FLERR,"Cannot change dump_modify every for dump xtc");
134 }
135
136 /* ---------------------------------------------------------------------- */
137
openfile()138 void DumpXTC::openfile()
139 {
140 // XTC maintains it's own XDR file ptr
141 // set fp to a null pointer so parent dump class will not use it
142
143 fp = nullptr;
144 if (me == 0)
145 if (xdropen(&xd,filename,"w") == 0) error->one(FLERR,"Cannot open dump file");
146 }
147
148 /* ---------------------------------------------------------------------- */
149
write_header(bigint nbig)150 void DumpXTC::write_header(bigint nbig)
151 {
152 if (nbig > MAXSMALLINT) error->all(FLERR,"Too many atoms for dump xtc");
153 int n = nbig;
154 if (update->ntimestep > MAXSMALLINT)
155 error->one(FLERR,"Too big a timestep for dump xtc");
156 int ntimestep = update->ntimestep;
157
158 // all procs realloc coords if total count grew
159
160 if (n != natoms) {
161 natoms = n;
162 memory->destroy(coords);
163 memory->create(coords,3*natoms,"dump:coords");
164 }
165
166 // only proc 0 writes header
167
168 if (me != 0) return;
169
170 int tmp = XTC_MAGIC;
171 xdr_int(&xd,&tmp);
172 xdr_int(&xd,&n);
173 xdr_int(&xd,&ntimestep);
174 float time_value = ntimestep * tfactor * update->dt;
175 xdr_float(&xd,&time_value);
176
177 // cell basis vectors
178 if (domain->triclinic) {
179 float zero = 0.0;
180 float xdim = sfactor * (domain->boxhi[0] - domain->boxlo[0]);
181 float ydim = sfactor * (domain->boxhi[1] - domain->boxlo[1]);
182 float zdim = sfactor * (domain->boxhi[2] - domain->boxlo[2]);
183 float xy = sfactor * domain->xy;
184 float xz = sfactor * domain->xz;
185 float yz = sfactor * domain->yz;
186
187 xdr_float(&xd,&xdim); xdr_float(&xd,&zero); xdr_float(&xd,&zero);
188 xdr_float(&xd,&xy ); xdr_float(&xd,&ydim); xdr_float(&xd,&zero);
189 xdr_float(&xd,&xz ); xdr_float(&xd,&yz ); xdr_float(&xd,&zdim);
190 } else {
191 float zero = 0.0;
192 float xdim = sfactor * (domain->boxhi[0] - domain->boxlo[0]);
193 float ydim = sfactor * (domain->boxhi[1] - domain->boxlo[1]);
194 float zdim = sfactor * (domain->boxhi[2] - domain->boxlo[2]);
195
196 xdr_float(&xd,&xdim); xdr_float(&xd,&zero); xdr_float(&xd,&zero);
197 xdr_float(&xd,&zero); xdr_float(&xd,&ydim); xdr_float(&xd,&zero);
198 xdr_float(&xd,&zero); xdr_float(&xd,&zero); xdr_float(&xd,&zdim);
199 }
200 }
201
202 /* ---------------------------------------------------------------------- */
203
pack(tagint * ids)204 void DumpXTC::pack(tagint *ids)
205 {
206 int m,n;
207
208 tagint *tag = atom->tag;
209 double **x = atom->x;
210 imageint *image = atom->image;
211 int *mask = atom->mask;
212 int nlocal = atom->nlocal;
213
214 m = n = 0;
215 if (unwrap_flag == 1) {
216 double xprd = domain->xprd;
217 double yprd = domain->yprd;
218 double zprd = domain->zprd;
219 double xy = domain->xy;
220 double xz = domain->xz;
221 double yz = domain->yz;
222
223 for (int i = 0; i < nlocal; i++)
224 if (mask[i] & groupbit) {
225 int ix = (image[i] & IMGMASK) - IMGMAX;
226 int iy = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
227 int iz = (image[i] >> IMG2BITS) - IMGMAX;
228
229 if (domain->triclinic) {
230 buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz);
231 buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz);
232 buf[m++] = sfactor * (x[i][2] + iz * zprd);
233 } else {
234 buf[m++] = sfactor * (x[i][0] + ix * xprd);
235 buf[m++] = sfactor * (x[i][1] + iy * yprd);
236 buf[m++] = sfactor * (x[i][2] + iz * zprd);
237 }
238 ids[n++] = tag[i];
239 }
240
241 } else {
242 for (int i = 0; i < nlocal; i++)
243 if (mask[i] & groupbit) {
244 buf[m++] = sfactor*x[i][0];
245 buf[m++] = sfactor*x[i][1];
246 buf[m++] = sfactor*x[i][2];
247 ids[n++] = tag[i];
248 }
249 }
250 }
251
252 /* ---------------------------------------------------------------------- */
253
write_data(int n,double * mybuf)254 void DumpXTC::write_data(int n, double *mybuf)
255 {
256 // copy buf atom coords into global array
257
258 int m = 0;
259 int k = 3*ntotal;
260 for (int i = 0; i < n; i++) {
261 coords[k++] = mybuf[m++];
262 coords[k++] = mybuf[m++];
263 coords[k++] = mybuf[m++];
264 ntotal++;
265 }
266
267 // if last chunk of atoms in this snapshot, write global arrays to file
268
269 if (ntotal == natoms) {
270 write_frame();
271 ntotal = 0;
272 }
273 }
274
275 /* ---------------------------------------------------------------------- */
276
modify_param(int narg,char ** arg)277 int DumpXTC::modify_param(int narg, char **arg)
278 {
279 if (strcmp(arg[0],"unwrap") == 0) {
280 if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
281 if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
282 else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
283 else error->all(FLERR,"Illegal dump_modify command");
284 return 2;
285 } else if (strcmp(arg[0],"precision") == 0) {
286 if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
287 precision = utils::numeric(FLERR,arg[1],false,lmp);
288 if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) &&
289 (fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) &&
290 (fabs(precision-100000.0) > EPS) &&
291 (fabs(precision-1000000.0) > EPS))
292 error->all(FLERR,"Illegal dump_modify command");
293 return 2;
294 } else if (strcmp(arg[0],"sfactor") == 0) {
295 if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
296 sfactor = utils::numeric(FLERR,arg[1],false,lmp);
297 if (sfactor <= 0.0)
298 error->all(FLERR,"Illegal dump_modify sfactor value (must be > 0.0)");
299 return 2;
300 } else if (strcmp(arg[0],"tfactor") == 0) {
301 if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
302 tfactor = utils::numeric(FLERR,arg[1],false,lmp);
303 if (tfactor <= 0.0)
304 error->all(FLERR,"Illegal dump_modify tfactor value (must be > 0.0)");
305 return 2;
306 }
307 return 0;
308 }
309
310 /* ----------------------------------------------------------------------
311 return # of bytes of allocated memory in buf and global coords array
312 ------------------------------------------------------------------------- */
313
memory_usage()314 double DumpXTC::memory_usage()
315 {
316 double bytes = Dump::memory_usage();
317 bytes += memory->usage(coords,natoms*3);
318 return bytes;
319 }
320
321 /* ---------------------------------------------------------------------- */
322
write_frame()323 void DumpXTC::write_frame()
324 {
325 xdr3dfcoord(&xd,coords,&natoms,&precision);
326 }
327
328 // ----------------------------------------------------------------------
329 // C functions that create GROMOS-compatible XDR files
330 // open-source
331 // (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
332 // ----------------------------------------------------------------------
333
334 /*____________________________________________________________________________
335 |
336 | Below are the routines to be used by C programmers. Use the 'normal'
337 | xdr routines to write integers, floats, etc (see man xdr)
338 |
339 | int xdropen(XDR *xdrs, const char *filename, const char *type)
340 | This will open the file with the given filename and the
341 | given mode. You should pass it an allocated XDR struct
342 | in xdrs, to be used in all other calls to xdr routines.
343 | Mode is 'w' to create, or update an file, and for all
344 | other values of mode the file is opened for reading.
345 | You need to call xdrclose to flush the output and close
346 | the file.
347 |
348 | Note that you should not use xdrstdio_create, which
349 | comes with the standard xdr library.
350 |
351 | int xdrclose(XDR *xdrs)
352 | Flush the data to the file, and close the file;
353 | You should not use xdr_destroy (which comes standard
354 | with the xdr libraries).
355 |
356 | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
357 | This is \fInot\fR a standard xdr routine. I named it this
358 | way, because it invites people to use the other xdr
359 | routines.
360 |
361 | (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
362 */
363
364 #define MAXID 20
365 static FILE *xdrfiles[MAXID];
366 static XDR *xdridptr[MAXID];
367 static char xdrmodes[MAXID];
368 static int *ip = nullptr;
369 static int *buf = nullptr;
370
371 /*___________________________________________________________________________
372 |
373 | what follows are the C routines for opening, closing xdr streams
374 | and the routine to read/write compressed coordinates together
375 | with some routines to assist in this task (those are marked
376 | static and cannot be called from user programs)
377 */
378 #define MAXABS INT_MAX-2
379
380 #ifndef SQR
381 #define SQR(x) ((x)*(x))
382 #endif
383 static int magicints[] = {
384 0, 0, 0, 0, 0, 0, 0, 0, 0,
385 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
386 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
387 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
388 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
389 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127,
390 524287, 660561,
391 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
392 4194304, 5284491, 6658042,
393 8388607, 10568983, 13316085, 16777216 };
394
395 #define FIRSTIDX 9
396 /* note that magicints[FIRSTIDX-1] == 0 */
397 #define LASTIDX (sizeof(magicints) / sizeof(int) - 1 )
398
399 /*__________________________________________________________________________
400 |
401 | xdropen - open xdr file
402 |
403 | This versions differs from xdrstdio_create, because I need to know
404 | the state of the file (read or write) so I can use xdr3dfcoord
405 | in eigther read or write mode, and the file descriptor
406 | so I can close the file (something xdr_destroy doesn't do).
407 |
408 */
409
xdropen(XDR * xdrs,const char * filename,const char * type)410 int xdropen(XDR *xdrs, const char *filename, const char *type)
411 {
412 static int init_done = 0;
413 enum xdr_op lmode;
414 int xdrid;
415
416 if (init_done == 0) {
417 for (xdrid = 1; xdrid < MAXID; xdrid++) {
418 xdridptr[xdrid] = nullptr;
419 }
420 init_done = 1;
421 }
422 xdrid = 1;
423 while (xdrid < MAXID && xdridptr[xdrid] != nullptr) {
424 xdrid++;
425 }
426 if (xdrid == MAXID) {
427 return 0;
428 }
429 if (*type == 'w' || *type == 'W') {
430 type = (char *) "w+";
431 lmode = XDR_ENCODE;
432 } else {
433 type = (char *) "r";
434 lmode = XDR_DECODE;
435 }
436 xdrfiles[xdrid] = fopen(filename, type);
437 if (xdrfiles[xdrid] == nullptr) {
438 xdrs = nullptr;
439 return 0;
440 }
441 xdrmodes[xdrid] = *type;
442
443 /* next test isn't useful in the case of C language
444 * but is used for the Fortran interface
445 * (C users are expected to pass the address of an already allocated
446 * XDR staructure)
447 */
448 if (xdrs == nullptr) {
449 xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
450 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
451 } else {
452 xdridptr[xdrid] = xdrs;
453 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
454 }
455 return xdrid;
456 }
457
458 /*_________________________________________________________________________
459 |
460 | xdrclose - close a xdr file
461 |
462 | This will flush the xdr buffers, and destroy the xdr stream.
463 | It also closes the associated file descriptor (this is *not*
464 | done by xdr_destroy).
465 |
466 */
467
xdrclose(XDR * xdrs)468 int xdrclose(XDR *xdrs)
469 {
470 int xdrid;
471
472 if (xdrs == nullptr) {
473 fprintf(stderr, "xdrclose: passed a NULL pointer\n");
474 exit(1);
475 }
476 for (xdrid = 1; xdrid < MAXID; xdrid++) {
477 if (xdridptr[xdrid] == xdrs) {
478
479 xdr_destroy(xdrs);
480 fclose(xdrfiles[xdrid]);
481 xdridptr[xdrid] = nullptr;
482 return 1;
483 }
484 }
485 fprintf(stderr, "xdrclose: no such open xdr file\n");
486 exit(1);
487 return 1;
488 }
489
490 /*_________________________________________________________________________
491 |
492 | xdrfreebuf - free the buffers used by xdr3dfcoord
493 |
494 */
xdrfreebuf()495 void xdrfreebuf()
496 {
497 if (ip) free(ip);
498 if (buf) free(buf);
499 ip = nullptr;
500 buf = nullptr;
501 }
502
503
504 /*____________________________________________________________________________
505 |
506 | sendbits - encode num into buf using the specified number of bits
507 |
508 | This routines appends the value of num to the bits already present in
509 | the array buf. You need to give it the number of bits to use and you
510 | better make sure that this number of bits is enough to hold the value
511 | Also num must be positive.
512 |
513 */
514
sendbits(int buf[],int num_of_bits,int num)515 static void sendbits(int buf[], int num_of_bits, int num)
516 {
517 unsigned int cnt, lastbyte;
518 int lastbits;
519 unsigned char * cbuf;
520
521 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
522 cnt = (unsigned int) buf[0];
523 lastbits = buf[1];
524 lastbyte =(unsigned int) buf[2];
525 while (num_of_bits >= 8) {
526 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
527 cbuf[cnt++] = lastbyte >> lastbits;
528 num_of_bits -= 8;
529 }
530 if (num_of_bits > 0) {
531 lastbyte = (lastbyte << num_of_bits) | num;
532 lastbits += num_of_bits;
533 if (lastbits >= 8) {
534 lastbits -= 8;
535 cbuf[cnt++] = lastbyte >> lastbits;
536 }
537 }
538 buf[0] = cnt;
539 buf[1] = lastbits;
540 buf[2] = lastbyte;
541 if (lastbits>0) {
542 cbuf[cnt] = lastbyte << (8 - lastbits);
543 }
544 }
545
546 /*_________________________________________________________________________
547 |
548 | sizeofint - calculate bitsize of an integer
549 |
550 | return the number of bits needed to store an integer with given max size
551 |
552 */
553
sizeofint(const int size)554 static int sizeofint(const int size)
555 {
556 unsigned int num = 1;
557 int num_of_bits = 0;
558
559 while (size >= (int) num && num_of_bits < 32) {
560 num_of_bits++;
561 num <<= 1;
562 }
563 return num_of_bits;
564 }
565
566 /*___________________________________________________________________________
567 |
568 | sizeofints - calculate 'bitsize' of compressed ints
569 |
570 | given the number of small unsigned integers and the maximum value
571 | return the number of bits needed to read or write them with the
572 | routines receiveints and sendints. You need this parameter when
573 | calling these routines. Note that for many calls I can use
574 | the variable 'smallidx' which is exactly the number of bits, and
575 | So I don't need to call 'sizeofints for those calls.
576 */
577
sizeofints(const int num_of_ints,unsigned int sizes[])578 static int sizeofints( const int num_of_ints, unsigned int sizes[])
579 {
580 int i, num;
581 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
582 num_of_bytes = 1;
583 bytes[0] = 1;
584 num_of_bits = 0;
585 for (i=0; i < num_of_ints; i++) {
586 tmp = 0;
587 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
588 tmp = bytes[bytecnt] * sizes[i] + tmp;
589 bytes[bytecnt] = tmp & 0xff;
590 tmp >>= 8;
591 }
592 while (tmp != 0) {
593 bytes[bytecnt++] = tmp & 0xff;
594 tmp >>= 8;
595 }
596 num_of_bytes = bytecnt;
597 }
598 num = 1;
599 num_of_bytes--;
600 while ((int)bytes[num_of_bytes] >= num) {
601 num_of_bits++;
602 num *= 2;
603 }
604 return num_of_bits + num_of_bytes * 8;
605 }
606
607 /*____________________________________________________________________________
608 |
609 | sendints - send a small set of small integers in compressed
610 |
611 | this routine is used internally by xdr3dfcoord, to send a set of
612 | small integers to the buffer.
613 | Multiplication with fixed (specified maximum ) sizes is used to get
614 | to one big, multibyte integer. Although the routine could be
615 | modified to handle sizes bigger than 16777216, or more than just
616 | a few integers, this is not done, because the gain in compression
617 | isn't worth the effort. Note that overflowing the multiplication
618 | or the byte buffer (32 bytes) is unchecked and causes bad results.
619 |
620 */
621
sendints(int buf[],const int num_of_ints,const int num_of_bits,unsigned int sizes[],unsigned int nums[])622 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
623 unsigned int sizes[], unsigned int nums[])
624 {
625 int i;
626 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
627
628 tmp = nums[0];
629 num_of_bytes = 0;
630 do {
631 bytes[num_of_bytes++] = tmp & 0xff;
632 tmp >>= 8;
633 } while (tmp != 0);
634
635 for (i = 1; i < num_of_ints; i++) {
636 if (nums[i] >= sizes[i]) {
637 fprintf(stderr,"major breakdown in sendints num %d doesn't "
638 "match size %d\n", nums[i], sizes[i]);
639 exit(1);
640 }
641 /* use one step multiply */
642 tmp = nums[i];
643 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
644 tmp = bytes[bytecnt] * sizes[i] + tmp;
645 bytes[bytecnt] = tmp & 0xff;
646 tmp >>= 8;
647 }
648 while (tmp != 0) {
649 bytes[bytecnt++] = tmp & 0xff;
650 tmp >>= 8;
651 }
652 num_of_bytes = bytecnt;
653 }
654 if (num_of_bits >= (int)num_of_bytes * 8) {
655 for (i = 0; i < (int)num_of_bytes; i++) {
656 sendbits(buf, 8, bytes[i]);
657 }
658 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
659 } else {
660 for (i = 0; i < (int)num_of_bytes-1; i++) {
661 sendbits(buf, 8, bytes[i]);
662 }
663 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
664 }
665 }
666
667 /*___________________________________________________________________________
668 |
669 | receivebits - decode number from buf using specified number of bits
670 |
671 | extract the number of bits from the array buf and construct an integer
672 | from it. Return that value.
673 |
674 */
675
receivebits(int buf[],int num_of_bits)676 static int receivebits(int buf[], int num_of_bits)
677 {
678 int cnt, num;
679 unsigned int lastbits, lastbyte;
680 unsigned char * cbuf;
681 int mask = (1 << num_of_bits) -1;
682
683 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
684 cnt = buf[0];
685 lastbits = (unsigned int) buf[1];
686 lastbyte = (unsigned int) buf[2];
687
688 num = 0;
689 while (num_of_bits >= 8) {
690 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
691 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
692 num_of_bits -=8;
693 }
694 if (num_of_bits > 0) {
695 if ((int)lastbits < num_of_bits) {
696 lastbits += 8;
697 lastbyte = (lastbyte << 8) | cbuf[cnt++];
698 }
699 lastbits -= num_of_bits;
700 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
701 }
702 num &= mask;
703 buf[0] = cnt;
704 buf[1] = lastbits;
705 buf[2] = lastbyte;
706 return num;
707 }
708
709 /*____________________________________________________________________________
710 |
711 | receiveints - decode 'small' integers from the buf array
712 |
713 | this routine is the inverse from sendints() and decodes the small integers
714 | written to buf by calculating the remainder and doing divisions with
715 | the given sizes[]. You need to specify the total number of bits to be
716 | used from buf in num_of_bits.
717 |
718 */
719
receiveints(int buf[],const int num_of_ints,int num_of_bits,unsigned int sizes[],int nums[])720 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
721 unsigned int sizes[], int nums[])
722 {
723 int bytes[32];
724 int i, j, num_of_bytes, p, num;
725
726 bytes[1] = bytes[2] = bytes[3] = 0;
727 num_of_bytes = 0;
728 while (num_of_bits > 8) {
729 bytes[num_of_bytes++] = receivebits(buf, 8);
730 num_of_bits -= 8;
731 }
732 if (num_of_bits > 0) {
733 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
734 }
735 for (i = num_of_ints-1; i > 0; i--) {
736 num = 0;
737 for (j = num_of_bytes-1; j >=0; j--) {
738 num = (num << 8) | bytes[j];
739 p = num / sizes[i];
740 bytes[j] = p;
741 num = num - p * sizes[i];
742 }
743 nums[i] = num;
744 }
745 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
746 }
747
748 /*____________________________________________________________________________
749 |
750 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
751 |
752 | this routine reads or writes (depending on how you opened the file with
753 | xdropen() ) a large number of 3d coordinates (stored in *fp).
754 | The number of coordinates triplets to write is given by *size. On
755 | read this number may be zero, in which case it reads as many as were written
756 | or it may specify the number if triplets to read (which should match the
757 | number written).
758 | Compression is achieved by first converting all floating numbers to integer
759 | using multiplication by *precision and rounding to the nearest integer.
760 | Then the minimum and maximum value are calculated to determine the range.
761 | The limited range of integers so found, is used to compress the coordinates.
762 | In addition the differences between successive coordinates is calculated.
763 | If the difference happens to be 'small' then only the difference is saved,
764 | compressing the data even more. The notion of 'small' is changed dynamically
765 | and is enlarged or reduced whenever needed or possible.
766 | Extra compression is achieved in the case of GROMOS and coordinates of
767 | water molecules. GROMOS first writes out the Oxygen position, followed by
768 | the two hydrogens. In order to make the differences smaller (and thereby
769 | compression the data better) the order is changed into first one hydrogen
770 | then the oxygen, followed by the other hydrogen. This is rather special, but
771 | it shouldn't harm in the general case.
772 |
773 */
774
xdr3dfcoord(XDR * xdrs,float * fp,int * size,float * precision)775 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
776 {
777 static int oldsize;
778
779 int minint[3], maxint[3], mindiff, *lip, diff;
780 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
781 int minidx, maxidx;
782 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
783 int flag, k;
784 int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
785 float *lfp, lf;
786 int tmp, *thiscoord, prevcoord[3];
787 unsigned int tmpcoord[30];
788
789 int bufsize, xdrid, lsize;
790 unsigned int bitsize;
791 float inv_precision;
792 int errval = 1;
793
794 /* find out if xdrs is opened for reading or for writing */
795 xdrid = 0;
796 while (xdridptr[xdrid] != xdrs) {
797 xdrid++;
798 if (xdrid >= MAXID) {
799 fprintf(stderr, "xdr error. no open xdr stream\n");
800 exit (1);
801 }
802 }
803 if (xdrmodes[xdrid] == 'w') {
804
805 /* xdrs is open for writing */
806
807 if (xdr_int(xdrs, size) == 0)
808 return 0;
809 size3 = *size * 3;
810 /* when the number of coordinates is small, don't try to compress; just
811 * write them as floats using xdr_vector
812 */
813 if (*size <= 9) {
814 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
815 (xdrproc_t)xdr_float));
816 }
817
818 xdr_float(xdrs, precision);
819 if (ip == nullptr) {
820 ip = (int *) malloc(size3 * sizeof(*ip));
821 if (ip == nullptr) {
822 fprintf(stderr,"malloc failed\n");
823 exit(1);
824 }
825 bufsize = (int) (size3 * 1.2);
826 buf = (int *) malloc(bufsize * sizeof(*buf));
827 if (buf == nullptr) {
828 fprintf(stderr,"malloc failed\n");
829 exit(1);
830 }
831 oldsize = *size;
832 } else if (*size > oldsize) {
833 ip = (int *) realloc(ip, size3 * sizeof(*ip));
834 if (ip == nullptr) {
835 fprintf(stderr,"malloc failed\n");
836 exit(1);
837 }
838 bufsize = (int) (size3 * 1.2);
839 buf = (int *) realloc(buf, bufsize * sizeof(*buf));
840 if (buf == nullptr) {
841 fprintf(stderr,"malloc failed\n");
842 exit(1);
843 }
844 oldsize = *size;
845 }
846 /* buf[0-2] are special and do not contain actual data */
847 buf[0] = buf[1] = buf[2] = 0;
848 minint[0] = minint[1] = minint[2] = INT_MAX;
849 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
850 prevrun = -1;
851 lfp = fp;
852 lip = ip;
853 mindiff = INT_MAX;
854 oldlint1 = oldlint2 = oldlint3 = 0;
855 while (lfp < fp + size3) {
856 /* find nearest integer */
857 if (*lfp >= 0.0)
858 lf = *lfp * *precision + 0.5;
859 else
860 lf = *lfp * *precision - 0.5;
861 if (fabs(lf) > MAXABS) {
862 /* scaling would cause overflow */
863 errval = 0;
864 }
865 lint1 = (int) lf;
866 if (lint1 < minint[0]) minint[0] = lint1;
867 if (lint1 > maxint[0]) maxint[0] = lint1;
868 *lip++ = lint1;
869 lfp++;
870 if (*lfp >= 0.0)
871 lf = *lfp * *precision + 0.5;
872 else
873 lf = *lfp * *precision - 0.5;
874 if (fabs(lf) > MAXABS) {
875 /* scaling would cause overflow */
876 errval = 0;
877 }
878 lint2 = (int) lf;
879 if (lint2 < minint[1]) minint[1] = lint2;
880 if (lint2 > maxint[1]) maxint[1] = lint2;
881 *lip++ = lint2;
882 lfp++;
883 if (*lfp >= 0.0)
884 lf = *lfp * *precision + 0.5;
885 else
886 lf = *lfp * *precision - 0.5;
887 if (fabs(lf) > MAXABS) {
888 /* scaling would cause overflow */
889 errval = 0;
890 }
891 lint3 = (int) lf;
892 if (lint3 < minint[2]) minint[2] = lint3;
893 if (lint3 > maxint[2]) maxint[2] = lint3;
894 *lip++ = lint3;
895 lfp++;
896 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
897 if (diff < mindiff && lfp > fp + 3)
898 mindiff = diff;
899 oldlint1 = lint1;
900 oldlint2 = lint2;
901 oldlint3 = lint3;
902 }
903 xdr_int(xdrs, &(minint[0]));
904 xdr_int(xdrs, &(minint[1]));
905 xdr_int(xdrs, &(minint[2]));
906
907 xdr_int(xdrs, &(maxint[0]));
908 xdr_int(xdrs, &(maxint[1]));
909 xdr_int(xdrs, &(maxint[2]));
910
911 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
912 (float)maxint[1] - (float)minint[1] >= MAXABS ||
913 (float)maxint[2] - (float)minint[2] >= MAXABS) {
914 /* turning value in unsigned by subtracting minint
915 * would cause overflow
916 */
917 errval = 0;
918 }
919 sizeint[0] = maxint[0] - minint[0]+1;
920 sizeint[1] = maxint[1] - minint[1]+1;
921 sizeint[2] = maxint[2] - minint[2]+1;
922
923 /* check if one of the sizes is to big to be multiplied */
924 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
925 bitsizeint[0] = sizeofint(sizeint[0]);
926 bitsizeint[1] = sizeofint(sizeint[1]);
927 bitsizeint[2] = sizeofint(sizeint[2]);
928 bitsize = 0; /* flag the use of large sizes */
929 } else {
930 bitsize = sizeofints(3, sizeint);
931 }
932 lip = ip;
933 luip = (unsigned int *) ip;
934 smallidx = FIRSTIDX;
935 while (smallidx < (int)LASTIDX && magicints[smallidx] < mindiff) {
936 smallidx++;
937 }
938 xdr_int(xdrs, &smallidx);
939 maxidx = MYMIN((int)LASTIDX, smallidx + 8) ;
940 minidx = maxidx - 8; /* often this equal smallidx */
941 smaller = magicints[MYMAX(FIRSTIDX, smallidx-1)] / 2;
942 small = magicints[smallidx] / 2;
943 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
944 larger = magicints[maxidx] / 2;
945 i = 0;
946 while (i < *size) {
947 is_small = 0;
948 thiscoord = (int *)(luip) + i * 3;
949 if (smallidx < maxidx && i >= 1 &&
950 abs(thiscoord[0] - prevcoord[0]) < larger &&
951 abs(thiscoord[1] - prevcoord[1]) < larger &&
952 abs(thiscoord[2] - prevcoord[2]) < larger) {
953 is_smaller = 1;
954 } else if (smallidx > minidx) {
955 is_smaller = -1;
956 } else {
957 is_smaller = 0;
958 }
959 if (i + 1 < *size) {
960 if (abs(thiscoord[0] - thiscoord[3]) < small &&
961 abs(thiscoord[1] - thiscoord[4]) < small &&
962 abs(thiscoord[2] - thiscoord[5]) < small) {
963 /* interchange first with second atom for better
964 * compression of water molecules
965 */
966 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
967 thiscoord[3] = tmp;
968 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
969 thiscoord[4] = tmp;
970 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
971 thiscoord[5] = tmp;
972 is_small = 1;
973 }
974
975 }
976 tmpcoord[0] = thiscoord[0] - minint[0];
977 tmpcoord[1] = thiscoord[1] - minint[1];
978 tmpcoord[2] = thiscoord[2] - minint[2];
979 if (bitsize == 0) {
980 sendbits(buf, bitsizeint[0], tmpcoord[0]);
981 sendbits(buf, bitsizeint[1], tmpcoord[1]);
982 sendbits(buf, bitsizeint[2], tmpcoord[2]);
983 } else {
984 sendints(buf, 3, bitsize, sizeint, tmpcoord);
985 }
986 prevcoord[0] = thiscoord[0];
987 prevcoord[1] = thiscoord[1];
988 prevcoord[2] = thiscoord[2];
989 thiscoord = thiscoord + 3;
990 i++;
991
992 run = 0;
993 if (is_small == 0 && is_smaller == -1)
994 is_smaller = 0;
995 while (is_small && run < 8*3) {
996 if (is_smaller == -1 && (SQR(thiscoord[0] - prevcoord[0]) +
997 SQR(thiscoord[1] - prevcoord[1]) +
998 SQR(thiscoord[2] - prevcoord[2]) >=
999 smaller * smaller)) {
1000 is_smaller = 0;
1001 }
1002
1003 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1004 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1005 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1006
1007 prevcoord[0] = thiscoord[0];
1008 prevcoord[1] = thiscoord[1];
1009 prevcoord[2] = thiscoord[2];
1010
1011 i++;
1012 thiscoord = thiscoord + 3;
1013 is_small = 0;
1014 if (i < *size &&
1015 abs(thiscoord[0] - prevcoord[0]) < small &&
1016 abs(thiscoord[1] - prevcoord[1]) < small &&
1017 abs(thiscoord[2] - prevcoord[2]) < small) {
1018 is_small = 1;
1019 }
1020 }
1021 if (run != prevrun || is_smaller != 0) {
1022 prevrun = run;
1023 sendbits(buf, 1, 1); /* flag the change in run-length */
1024 sendbits(buf, 5, run+is_smaller+1);
1025 } else {
1026 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1027 }
1028 for (k=0; k < run; k+=3) {
1029 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1030 }
1031 if (is_smaller != 0) {
1032 smallidx += is_smaller;
1033 if (is_smaller < 0) {
1034 small = smaller;
1035 smaller = magicints[smallidx-1] / 2;
1036 } else {
1037 smaller = small;
1038 small = magicints[smallidx] / 2;
1039 }
1040 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1041 }
1042 }
1043 if (buf[1] != 0) buf[0]++;;
1044 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1045 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1046 } else {
1047
1048 /* xdrs is open for reading */
1049
1050 if (xdr_int(xdrs, &lsize) == 0)
1051 return 0;
1052 if (*size != 0 && lsize != *size) {
1053 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1054 "%d arg vs %d in file", *size, lsize);
1055 }
1056 *size = lsize;
1057 size3 = *size * 3;
1058 if (*size <= 9) {
1059 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1060 (xdrproc_t)xdr_float));
1061 }
1062 xdr_float(xdrs, precision);
1063 if (ip == nullptr) {
1064 ip = (int *) malloc(size3 * sizeof(*ip));
1065 if (ip == nullptr) {
1066 fprintf(stderr,"malloc failed\n");
1067 exit(1);
1068 }
1069 bufsize = (int) (size3 * 1.2);
1070 buf = (int *) malloc(bufsize * sizeof(*buf));
1071 if (buf == nullptr) {
1072 fprintf(stderr,"malloc failed\n");
1073 exit(1);
1074 }
1075 oldsize = *size;
1076 } else if (*size > oldsize) {
1077 ip = (int *)realloc(ip, size3 * sizeof(*ip));
1078 if (ip == nullptr) {
1079 fprintf(stderr,"malloc failed\n");
1080 exit(1);
1081 }
1082 bufsize = (int) (size3 * 1.2);
1083 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1084 if (buf == nullptr) {
1085 fprintf(stderr,"malloc failed\n");
1086 exit(1);
1087 }
1088 oldsize = *size;
1089 }
1090 buf[0] = buf[1] = buf[2] = 0;
1091
1092 xdr_int(xdrs, &(minint[0]));
1093 xdr_int(xdrs, &(minint[1]));
1094 xdr_int(xdrs, &(minint[2]));
1095
1096 xdr_int(xdrs, &(maxint[0]));
1097 xdr_int(xdrs, &(maxint[1]));
1098 xdr_int(xdrs, &(maxint[2]));
1099
1100 sizeint[0] = maxint[0] - minint[0]+1;
1101 sizeint[1] = maxint[1] - minint[1]+1;
1102 sizeint[2] = maxint[2] - minint[2]+1;
1103
1104 /* check if one of the sizes is to big to be multiplied */
1105 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1106 bitsizeint[0] = sizeofint(sizeint[0]);
1107 bitsizeint[1] = sizeofint(sizeint[1]);
1108 bitsizeint[2] = sizeofint(sizeint[2]);
1109 bitsize = 0; /* flag the use of large sizes */
1110 } else {
1111 bitsize = sizeofints(3, sizeint);
1112 }
1113
1114 xdr_int(xdrs, &smallidx);
1115 maxidx = MYMIN((int)LASTIDX, smallidx + 8) ;
1116 minidx = maxidx - 8; /* often this equal smallidx */
1117 smaller = magicints[MYMAX(FIRSTIDX, smallidx-1)] / 2;
1118 small = magicints[smallidx] / 2;
1119 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1120 larger = magicints[maxidx];
1121
1122 /* buf[0] holds the length in bytes */
1123
1124 if (xdr_int(xdrs, &(buf[0])) == 0)
1125 return 0;
1126 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1127 return 0;
1128 buf[0] = buf[1] = buf[2] = 0;
1129
1130 lfp = fp;
1131 inv_precision = 1.0 / * precision;
1132 run = 0;
1133 i = 0;
1134 lip = ip;
1135 while (i < lsize) {
1136 thiscoord = (int *)(lip) + i * 3;
1137
1138 if (bitsize == 0) {
1139 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1140 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1141 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1142 } else {
1143 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1144 }
1145
1146 i++;
1147 thiscoord[0] += minint[0];
1148 thiscoord[1] += minint[1];
1149 thiscoord[2] += minint[2];
1150
1151 prevcoord[0] = thiscoord[0];
1152 prevcoord[1] = thiscoord[1];
1153 prevcoord[2] = thiscoord[2];
1154
1155
1156 flag = receivebits(buf, 1);
1157 is_smaller = 0;
1158 if (flag == 1) {
1159 run = receivebits(buf, 5);
1160 is_smaller = run % 3;
1161 run -= is_smaller;
1162 is_smaller--;
1163 }
1164 if (run > 0) {
1165 thiscoord += 3;
1166 for (k = 0; k < run; k+=3) {
1167 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1168 i++;
1169 thiscoord[0] += prevcoord[0] - small;
1170 thiscoord[1] += prevcoord[1] - small;
1171 thiscoord[2] += prevcoord[2] - small;
1172 if (k == 0) {
1173 /* interchange first with second atom for better
1174 * compression of water molecules
1175 */
1176 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1177 prevcoord[0] = tmp;
1178 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1179 prevcoord[1] = tmp;
1180 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1181 prevcoord[2] = tmp;
1182 *lfp++ = prevcoord[0] * inv_precision;
1183 *lfp++ = prevcoord[1] * inv_precision;
1184 *lfp++ = prevcoord[2] * inv_precision;
1185 } else {
1186 prevcoord[0] = thiscoord[0];
1187 prevcoord[1] = thiscoord[1];
1188 prevcoord[2] = thiscoord[2];
1189 }
1190 *lfp++ = thiscoord[0] * inv_precision;
1191 *lfp++ = thiscoord[1] * inv_precision;
1192 *lfp++ = thiscoord[2] * inv_precision;
1193 }
1194 } else {
1195 *lfp++ = thiscoord[0] * inv_precision;
1196 *lfp++ = thiscoord[1] * inv_precision;
1197 *lfp++ = thiscoord[2] * inv_precision;
1198 }
1199 smallidx += is_smaller;
1200 if (is_smaller < 0) {
1201 small = smaller;
1202 if (smallidx > FIRSTIDX) {
1203 smaller = magicints[smallidx - 1] /2;
1204 } else {
1205 smaller = 0;
1206 }
1207 } else if (is_smaller > 0) {
1208 smaller = small;
1209 small = magicints[smallidx] / 2;
1210 }
1211 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1212 }
1213 }
1214 return 1;
1215 }
1216