From 19122f203ea31a0776f195adb61ef30eb18b2cfb Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 28 Apr 2016 14:49:49 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14902 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MISC/dump_xtc.cpp | 1214 +++++++++++++++++++++++++++++++++++++++ src/MISC/dump_xtc.h | 118 ++++ src/MISC/xdr_compat.cpp | 722 +++++++++++++++++++++++ src/MISC/xdr_compat.h | 228 ++++++++ 4 files changed, 2282 insertions(+) create mode 100644 src/MISC/dump_xtc.cpp create mode 100644 src/MISC/dump_xtc.h create mode 100644 src/MISC/xdr_compat.cpp create mode 100644 src/MISC/xdr_compat.h diff --git a/src/MISC/dump_xtc.cpp b/src/MISC/dump_xtc.cpp new file mode 100644 index 0000000000..cd68e969af --- /dev/null +++ b/src/MISC/dump_xtc.cpp @@ -0,0 +1,1214 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U) + open-source XDR routines from + Frans van Hoesel (http://md.chem.rug.nl/hoesel) + are included in this file + Axel Kohlmeyer (Temple U) + port to platforms without XDR support + added support for unwrapped trajectories + support for groups +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include +#include "dump_xtc.h" +#include "domain.h" +#include "atom.h" +#include "update.h" +#include "group.h" +#include "output.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define EPS 1e-5 +#define XTC_MAGIC 1995 + +#define MYMIN(a,b) ((a) < (b) ? (a) : (b)) +#define MYMAX(a,b) ((a) > (b) ? (a) : (b)) + +int xdropen(XDR *, const char *, const char *); +int xdrclose(XDR *); +void xdrfreebuf(); +int xdr3dfcoord(XDR *, float *, int *, float *); + +/* ---------------------------------------------------------------------- */ + +DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) +{ + if (narg != 5) error->all(FLERR,"Illegal dump xtc command"); + if (binary || compressed || multifile || multiproc) + error->all(FLERR,"Invalid dump xtc filename"); + + size_one = 3; + sort_flag = 1; + sortcol = 0; + format_default = NULL; + flush_flag = 0; + unwrap_flag = 0; + precision = 1000.0; + + // allocate global array for atom coords + + bigint n = group->count(igroup); + if (n > static_cast(MAXSMALLINT/3/sizeof(float))) + error->all(FLERR,"Too many atoms for dump xtc"); + natoms = static_cast (n); + + memory->create(coords,3*natoms,"dump:coords"); + + // sfactor = conversion of coords to XTC units + // tfactor = conversion of simulation time to XTC units + // GROMACS standard is nanometers and picoseconds + + sfactor = 0.1 / force->angstrom; + tfactor = 0.001 / force->femtosecond; + + // in reduced units we do not scale anything + if (strcmp(update->unit_style,"lj") == 0) { + sfactor = tfactor = 1.0; + if (comm->me == 0) + error->warning(FLERR,"No automatic unit conversion to XTC file " + "format conventions possible for units lj"); + } + + openfile(); + nevery_save = 0; + ntotal = 0; +} + +/* ---------------------------------------------------------------------- */ + +DumpXTC::~DumpXTC() +{ + memory->destroy(coords); + + if (me == 0) { + xdrclose(&xd); + xdrfreebuf(); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::init_style() +{ + if (sort_flag == 0 || sortcol != 0) + error->all(FLERR,"Dump xtc requires sorting by atom ID"); + + // check that flush_flag is not set since dump::write() will use it + + if (flush_flag) error->all(FLERR,"Cannot set dump_modify flush for dump xtc"); + + // check that dump frequency has not changed and is not a variable + + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + if (output->every_dump[idump] == 0) + error->all(FLERR,"Cannot use variable every setting for dump xtc"); + + if (nevery_save == 0) nevery_save = output->every_dump[idump]; + else if (nevery_save != output->every_dump[idump]) + error->all(FLERR,"Cannot change dump_modify every for dump xtc"); +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::openfile() +{ + // XTC maintains it's own XDR file ptr + // set fp to NULL so parent dump class will not use it + + fp = NULL; + if (me == 0) + if (xdropen(&xd,filename,"w") == 0) error->one(FLERR,"Cannot open dump file"); +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::write_header(bigint nbig) +{ + if (nbig > MAXSMALLINT) error->all(FLERR,"Too many atoms for dump xtc"); + int n = nbig; + if (update->ntimestep > MAXSMALLINT) + error->one(FLERR,"Too big a timestep for dump xtc"); + int ntimestep = update->ntimestep; + + // all procs realloc coords if total count grew + + if (n != natoms) { + natoms = n; + memory->destroy(coords); + memory->create(coords,3*natoms,"dump:coords"); + } + + // only proc 0 writes header + + if (me != 0) return; + + int tmp = XTC_MAGIC; + xdr_int(&xd,&tmp); + xdr_int(&xd,&n); + xdr_int(&xd,&ntimestep); + float time_value = ntimestep * tfactor * update->dt; + xdr_float(&xd,&time_value); + + // cell basis vectors + if (domain->triclinic) { + float zero = 0.0; + float xdim = sfactor * (domain->boxhi[0] - domain->boxlo[0]); + float ydim = sfactor * (domain->boxhi[1] - domain->boxlo[1]); + float zdim = sfactor * (domain->boxhi[2] - domain->boxlo[2]); + float xy = sfactor * domain->xy; + float xz = sfactor * domain->xz; + float yz = sfactor * domain->yz; + + xdr_float(&xd,&xdim); xdr_float(&xd,&zero); xdr_float(&xd,&zero); + xdr_float(&xd,&xy ); xdr_float(&xd,&ydim); xdr_float(&xd,&zero); + xdr_float(&xd,&xz ); xdr_float(&xd,&yz ); xdr_float(&xd,&zdim); + } else { + float zero = 0.0; + float xdim = sfactor * (domain->boxhi[0] - domain->boxlo[0]); + float ydim = sfactor * (domain->boxhi[1] - domain->boxlo[1]); + float zdim = sfactor * (domain->boxhi[2] - domain->boxlo[2]); + + xdr_float(&xd,&xdim); xdr_float(&xd,&zero); xdr_float(&xd,&zero); + xdr_float(&xd,&zero); xdr_float(&xd,&ydim); xdr_float(&xd,&zero); + xdr_float(&xd,&zero); xdr_float(&xd,&zero); xdr_float(&xd,&zdim); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::pack(tagint *ids) +{ + int m,n; + + tagint *tag = atom->tag; + double **x = atom->x; + imageint *image = atom->image; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + m = n = 0; + if (unwrap_flag == 1) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xy = domain->xy; + double xz = domain->xz; + double yz = domain->yz; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + int ix = (image[i] & IMGMASK) - IMGMAX; + int iy = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + int iz = (image[i] >> IMG2BITS) - IMGMAX; + + if (domain->triclinic) { + buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz); + buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz); + buf[m++] = sfactor * (x[i][2] + iz * zprd); + } else { + buf[m++] = sfactor * (x[i][0] + ix * xprd); + buf[m++] = sfactor * (x[i][1] + iy * yprd); + buf[m++] = sfactor * (x[i][2] + iz * zprd); + } + ids[n++] = tag[i]; + } + + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + buf[m++] = sfactor*x[i][0]; + buf[m++] = sfactor*x[i][1]; + buf[m++] = sfactor*x[i][2]; + ids[n++] = tag[i]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::write_data(int n, double *mybuf) +{ + // copy buf atom coords into global array + + int m = 0; + int k = 3*ntotal; + for (int i = 0; i < n; i++) { + coords[k++] = mybuf[m++]; + coords[k++] = mybuf[m++]; + coords[k++] = mybuf[m++]; + ntotal++; + } + + // if last chunk of atoms in this snapshot, write global arrays to file + + if (ntotal == natoms) { + write_frame(); + ntotal = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpXTC::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; + else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; + else error->all(FLERR,"Illegal dump_modify command"); + return 2; + } else if (strcmp(arg[0],"precision") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + precision = force->numeric(FLERR,arg[1]); + if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) && + (fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) && + (fabs(precision-100000.0) > EPS) && + (fabs(precision-1000000.0) > EPS)) + error->all(FLERR,"Illegal dump_modify command"); + return 2; + } else if (strcmp(arg[0],"sfactor") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + sfactor = force->numeric(FLERR,arg[1]); + if (sfactor <= 0.0) + error->all(FLERR,"Illegal dump_modify sfactor value (must be > 0.0)"); + return 2; + } else if (strcmp(arg[0],"tfactor") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + tfactor = force->numeric(FLERR,arg[1]); + if (tfactor <= 0.0) + error->all(FLERR,"Illegal dump_modify tfactor value (must be > 0.0)"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory in buf and global coords array +------------------------------------------------------------------------- */ + +bigint DumpXTC::memory_usage() +{ + bigint bytes = Dump::memory_usage(); + bytes += memory->usage(coords,natoms*3); + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +void DumpXTC::write_frame() +{ + xdr3dfcoord(&xd,coords,&natoms,&precision); +} + +// ---------------------------------------------------------------------- +// C functions that create GROMOS-compatible XDR files +// open-source +// (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl +// ---------------------------------------------------------------------- + +/*____________________________________________________________________________ + | + | Below are the routines to be used by C programmers. Use the 'normal' + | xdr routines to write integers, floats, etc (see man xdr) + | + | int xdropen(XDR *xdrs, const char *filename, const char *type) + | This will open the file with the given filename and the + | given mode. You should pass it an allocated XDR struct + | in xdrs, to be used in all other calls to xdr routines. + | Mode is 'w' to create, or update an file, and for all + | other values of mode the file is opened for reading. + | You need to call xdrclose to flush the output and close + | the file. + | + | Note that you should not use xdrstdio_create, which + | comes with the standard xdr library. + | + | int xdrclose(XDR *xdrs) + | Flush the data to the file, and close the file; + | You should not use xdr_destroy (which comes standard + | with the xdr libraries). + | + | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) + | This is \fInot\fR a standard xdr routine. I named it this + | way, because it invites people to use the other xdr + | routines. + | + | (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl +*/ + +#define MAXID 20 +static FILE *xdrfiles[MAXID]; +static XDR *xdridptr[MAXID]; +static char xdrmodes[MAXID]; +static int *ip = NULL; +static int *buf = NULL; + +/*___________________________________________________________________________ + | + | what follows are the C routines for opening, closing xdr streams + | and the routine to read/write compressed coordinates together + | with some routines to assist in this task (those are marked + | static and cannot be called from user programs) +*/ +#define MAXABS INT_MAX-2 + +#ifndef SQR +#define SQR(x) ((x)*(x)) +#endif +static int magicints[] = { + 0, 0, 0, 0, 0, 0, 0, 0, 0, + 8, 10, 12, 16, 20, 25, 32, 40, 50, 64, + 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, + 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501, + 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536, + 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, + 524287, 660561, + 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, + 4194304, 5284491, 6658042, + 8388607, 10568983, 13316085, 16777216 }; + +#define FIRSTIDX 9 +/* note that magicints[FIRSTIDX-1] == 0 */ +#define LASTIDX (sizeof(magicints) / sizeof(int) - 1 ) + +/*__________________________________________________________________________ + | + | xdropen - open xdr file + | + | This versions differs from xdrstdio_create, because I need to know + | the state of the file (read or write) so I can use xdr3dfcoord + | in eigther read or write mode, and the file descriptor + | so I can close the file (something xdr_destroy doesn't do). + | +*/ + +int xdropen(XDR *xdrs, const char *filename, const char *type) +{ + static int init_done = 0; + enum xdr_op lmode; + int xdrid; + + if (init_done == 0) { + for (xdrid = 1; xdrid < MAXID; xdrid++) { + xdridptr[xdrid] = NULL; + } + init_done = 1; + } + xdrid = 1; + while (xdrid < MAXID && xdridptr[xdrid] != NULL) { + xdrid++; + } + if (xdrid == MAXID) { + return 0; + } + if (*type == 'w' || *type == 'W') { + type = (char *) "w+"; + lmode = XDR_ENCODE; + } else { + type = (char *) "r"; + lmode = XDR_DECODE; + } + xdrfiles[xdrid] = fopen(filename, type); + if (xdrfiles[xdrid] == NULL) { + xdrs = NULL; + return 0; + } + xdrmodes[xdrid] = *type; + + /* next test isn't usefull in the case of C language + * but is used for the Fortran interface + * (C users are expected to pass the address of an already allocated + * XDR staructure) + */ + if (xdrs == NULL) { + xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR)); + xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode); + } else { + xdridptr[xdrid] = xdrs; + xdrstdio_create(xdrs, xdrfiles[xdrid], lmode); + } + return xdrid; +} + +/*_________________________________________________________________________ + | + | xdrclose - close a xdr file + | + | This will flush the xdr buffers, and destroy the xdr stream. + | It also closes the associated file descriptor (this is *not* + | done by xdr_destroy). + | +*/ + +int xdrclose(XDR *xdrs) +{ + int xdrid; + + if (xdrs == NULL) { + fprintf(stderr, "xdrclose: passed a NULL pointer\n"); + exit(1); + } + for (xdrid = 1; xdrid < MAXID; xdrid++) { + if (xdridptr[xdrid] == xdrs) { + + xdr_destroy(xdrs); + fclose(xdrfiles[xdrid]); + xdridptr[xdrid] = NULL; + return 1; + } + } + fprintf(stderr, "xdrclose: no such open xdr file\n"); + exit(1); + return 1; +} + +/*_________________________________________________________________________ + | + | xdrfreebuf - free the buffers used by xdr3dfcoord + | +*/ +void xdrfreebuf() +{ + if (ip) free(ip); + if (buf) free(buf); + ip = NULL; + buf = NULL; +} + + +/*____________________________________________________________________________ + | + | sendbits - encode num into buf using the specified number of bits + | + | This routines appends the value of num to the bits already present in + | the array buf. You need to give it the number of bits to use and you + | better make sure that this number of bits is enough to hold the value + | Also num must be positive. + | +*/ + +static void sendbits(int buf[], int num_of_bits, int num) +{ + unsigned int cnt, lastbyte; + int lastbits; + unsigned char * cbuf; + + cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); + cnt = (unsigned int) buf[0]; + lastbits = buf[1]; + lastbyte =(unsigned int) buf[2]; + while (num_of_bits >= 8) { + lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/); + cbuf[cnt++] = lastbyte >> lastbits; + num_of_bits -= 8; + } + if (num_of_bits > 0) { + lastbyte = (lastbyte << num_of_bits) | num; + lastbits += num_of_bits; + if (lastbits >= 8) { + lastbits -= 8; + cbuf[cnt++] = lastbyte >> lastbits; + } + } + buf[0] = cnt; + buf[1] = lastbits; + buf[2] = lastbyte; + if (lastbits>0) { + cbuf[cnt] = lastbyte << (8 - lastbits); + } +} + +/*_________________________________________________________________________ + | + | sizeofint - calculate bitsize of an integer + | + | return the number of bits needed to store an integer with given max size + | +*/ + +static int sizeofint(const int size) +{ + unsigned int num = 1; + int num_of_bits = 0; + + while (size >= num && num_of_bits < 32) { + num_of_bits++; + num <<= 1; + } + return num_of_bits; +} + +/*___________________________________________________________________________ + | + | sizeofints - calculate 'bitsize' of compressed ints + | + | given the number of small unsigned integers and the maximum value + | return the number of bits needed to read or write them with the + | routines receiveints and sendints. You need this parameter when + | calling these routines. Note that for many calls I can use + | the variable 'smallidx' which is exactly the number of bits, and + | So I don't need to call 'sizeofints for those calls. +*/ + +static int sizeofints( const int num_of_ints, unsigned int sizes[]) +{ + int i, num; + unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp; + num_of_bytes = 1; + bytes[0] = 1; + num_of_bits = 0; + for (i=0; i < num_of_ints; i++) { + tmp = 0; + for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { + tmp = bytes[bytecnt] * sizes[i] + tmp; + bytes[bytecnt] = tmp & 0xff; + tmp >>= 8; + } + while (tmp != 0) { + bytes[bytecnt++] = tmp & 0xff; + tmp >>= 8; + } + num_of_bytes = bytecnt; + } + num = 1; + num_of_bytes--; + while (bytes[num_of_bytes] >= num) { + num_of_bits++; + num *= 2; + } + return num_of_bits + num_of_bytes * 8; +} + +/*____________________________________________________________________________ + | + | sendints - send a small set of small integers in compressed + | + | this routine is used internally by xdr3dfcoord, to send a set of + | small integers to the buffer. + | Multiplication with fixed (specified maximum ) sizes is used to get + | to one big, multibyte integer. Allthough the routine could be + | modified to handle sizes bigger than 16777216, or more than just + | a few integers, this is not done, because the gain in compression + | isn't worth the effort. Note that overflowing the multiplication + | or the byte buffer (32 bytes) is unchecked and causes bad results. + | + */ + +static void sendints(int buf[], const int num_of_ints, const int num_of_bits, + unsigned int sizes[], unsigned int nums[]) +{ + int i; + unsigned int bytes[32], num_of_bytes, bytecnt, tmp; + + tmp = nums[0]; + num_of_bytes = 0; + do { + bytes[num_of_bytes++] = tmp & 0xff; + tmp >>= 8; + } while (tmp != 0); + + for (i = 1; i < num_of_ints; i++) { + if (nums[i] >= sizes[i]) { + fprintf(stderr,"major breakdown in sendints num %d doesn't " + "match size %d\n", nums[i], sizes[i]); + exit(1); + } + /* use one step multiply */ + tmp = nums[i]; + for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { + tmp = bytes[bytecnt] * sizes[i] + tmp; + bytes[bytecnt] = tmp & 0xff; + tmp >>= 8; + } + while (tmp != 0) { + bytes[bytecnt++] = tmp & 0xff; + tmp >>= 8; + } + num_of_bytes = bytecnt; + } + if (num_of_bits >= num_of_bytes * 8) { + for (i = 0; i < num_of_bytes; i++) { + sendbits(buf, 8, bytes[i]); + } + sendbits(buf, num_of_bits - num_of_bytes * 8, 0); + } else { + for (i = 0; i < num_of_bytes-1; i++) { + sendbits(buf, 8, bytes[i]); + } + sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]); + } +} + +/*___________________________________________________________________________ + | + | receivebits - decode number from buf using specified number of bits + | + | extract the number of bits from the array buf and construct an integer + | from it. Return that value. + | +*/ + +static int receivebits(int buf[], int num_of_bits) +{ + int cnt, num; + unsigned int lastbits, lastbyte; + unsigned char * cbuf; + int mask = (1 << num_of_bits) -1; + + cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); + cnt = buf[0]; + lastbits = (unsigned int) buf[1]; + lastbyte = (unsigned int) buf[2]; + + num = 0; + while (num_of_bits >= 8) { + lastbyte = ( lastbyte << 8 ) | cbuf[cnt++]; + num |= (lastbyte >> lastbits) << (num_of_bits - 8); + num_of_bits -=8; + } + if (num_of_bits > 0) { + if (lastbits < num_of_bits) { + lastbits += 8; + lastbyte = (lastbyte << 8) | cbuf[cnt++]; + } + lastbits -= num_of_bits; + num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1); + } + num &= mask; + buf[0] = cnt; + buf[1] = lastbits; + buf[2] = lastbyte; + return num; +} + +/*____________________________________________________________________________ + | + | receiveints - decode 'small' integers from the buf array + | + | this routine is the inverse from sendints() and decodes the small integers + | written to buf by calculating the remainder and doing divisions with + | the given sizes[]. You need to specify the total number of bits to be + | used from buf in num_of_bits. + | +*/ + +static void receiveints(int buf[], const int num_of_ints, int num_of_bits, + unsigned int sizes[], int nums[]) +{ + int bytes[32]; + int i, j, num_of_bytes, p, num; + + bytes[1] = bytes[2] = bytes[3] = 0; + num_of_bytes = 0; + while (num_of_bits > 8) { + bytes[num_of_bytes++] = receivebits(buf, 8); + num_of_bits -= 8; + } + if (num_of_bits > 0) { + bytes[num_of_bytes++] = receivebits(buf, num_of_bits); + } + for (i = num_of_ints-1; i > 0; i--) { + num = 0; + for (j = num_of_bytes-1; j >=0; j--) { + num = (num << 8) | bytes[j]; + p = num / sizes[i]; + bytes[j] = p; + num = num - p * sizes[i]; + } + nums[i] = num; + } + nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24); +} + +/*____________________________________________________________________________ + | + | xdr3dfcoord - read or write compressed 3d coordinates to xdr file. + | + | this routine reads or writes (depending on how you opened the file with + | xdropen() ) a large number of 3d coordinates (stored in *fp). + | The number of coordinates triplets to write is given by *size. On + | read this number may be zero, in which case it reads as many as were written + | or it may specify the number if triplets to read (which should match the + | number written). + | Compression is achieved by first converting all floating numbers to integer + | using multiplication by *precision and rounding to the nearest integer. + | Then the minimum and maximum value are calculated to determine the range. + | The limited range of integers so found, is used to compress the coordinates. + | In addition the differences between succesive coordinates is calculated. + | If the difference happens to be 'small' then only the difference is saved, + | compressing the data even more. The notion of 'small' is changed dynamically + | and is enlarged or reduced whenever needed or possible. + | Extra compression is achieved in the case of GROMOS and coordinates of + | water molecules. GROMOS first writes out the Oxygen position, followed by + | the two hydrogens. In order to make the differences smaller (and thereby + | compression the data better) the order is changed into first one hydrogen + | then the oxygen, followed by the other hydrogen. This is rather special, but + | it shouldn't harm in the general case. + | + */ + +int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) +{ + static int oldsize; + + int minint[3], maxint[3], mindiff, *lip, diff; + int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx; + int minidx, maxidx; + unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip; + int flag, k; + int small, smaller, larger, i, is_small, is_smaller, run, prevrun; + float *lfp, lf; + int tmp, *thiscoord, prevcoord[3]; + unsigned int tmpcoord[30]; + + int bufsize, xdrid, lsize; + unsigned int bitsize; + float inv_precision; + int errval = 1; + + /* find out if xdrs is opened for reading or for writing */ + xdrid = 0; + while (xdridptr[xdrid] != xdrs) { + xdrid++; + if (xdrid >= MAXID) { + fprintf(stderr, "xdr error. no open xdr stream\n"); + exit (1); + } + } + if (xdrmodes[xdrid] == 'w') { + + /* xdrs is open for writing */ + + if (xdr_int(xdrs, size) == 0) + return 0; + size3 = *size * 3; + /* when the number of coordinates is small, don't try to compress; just + * write them as floats using xdr_vector + */ + if (*size <= 9 ) { + return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), + (xdrproc_t)xdr_float)); + } + + xdr_float(xdrs, precision); + if (ip == NULL) { + ip = (int *) malloc(size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = (int) (size3 * 1.2); + buf = (int *) malloc(bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } else if (*size > oldsize) { + ip = (int *) realloc(ip, size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = (int) (size3 * 1.2); + buf = (int *) realloc(buf, bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } + /* buf[0-2] are special and do not contain actual data */ + buf[0] = buf[1] = buf[2] = 0; + minint[0] = minint[1] = minint[2] = INT_MAX; + maxint[0] = maxint[1] = maxint[2] = INT_MIN; + prevrun = -1; + lfp = fp; + lip = ip; + mindiff = INT_MAX; + oldlint1 = oldlint2 = oldlint3 = 0; + while(lfp < fp + size3 ) { + /* find nearest integer */ + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint1 = (int) lf; + if (lint1 < minint[0]) minint[0] = lint1; + if (lint1 > maxint[0]) maxint[0] = lint1; + *lip++ = lint1; + lfp++; + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint2 = (int) lf; + if (lint2 < minint[1]) minint[1] = lint2; + if (lint2 > maxint[1]) maxint[1] = lint2; + *lip++ = lint2; + lfp++; + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint3 = (int) lf; + if (lint3 < minint[2]) minint[2] = lint3; + if (lint3 > maxint[2]) maxint[2] = lint3; + *lip++ = lint3; + lfp++; + diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3); + if (diff < mindiff && lfp > fp + 3) + mindiff = diff; + oldlint1 = lint1; + oldlint2 = lint2; + oldlint3 = lint3; + } + xdr_int(xdrs, &(minint[0])); + xdr_int(xdrs, &(minint[1])); + xdr_int(xdrs, &(minint[2])); + + xdr_int(xdrs, &(maxint[0])); + xdr_int(xdrs, &(maxint[1])); + xdr_int(xdrs, &(maxint[2])); + + if ((float)maxint[0] - (float)minint[0] >= MAXABS || + (float)maxint[1] - (float)minint[1] >= MAXABS || + (float)maxint[2] - (float)minint[2] >= MAXABS) { + /* turning value in unsigned by subtracting minint + * would cause overflow + */ + errval = 0; + } + sizeint[0] = maxint[0] - minint[0]+1; + sizeint[1] = maxint[1] - minint[1]+1; + sizeint[2] = maxint[2] - minint[2]+1; + + /* check if one of the sizes is to big to be multiplied */ + if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { + bitsizeint[0] = sizeofint(sizeint[0]); + bitsizeint[1] = sizeofint(sizeint[1]); + bitsizeint[2] = sizeofint(sizeint[2]); + bitsize = 0; /* flag the use of large sizes */ + } else { + bitsize = sizeofints(3, sizeint); + } + lip = ip; + luip = (unsigned int *) ip; + smallidx = FIRSTIDX; + while (smallidx < LASTIDX && magicints[smallidx] < mindiff) { + smallidx++; + } + xdr_int(xdrs, &smallidx); + maxidx = MYMIN(LASTIDX, smallidx + 8) ; + minidx = maxidx - 8; /* often this equal smallidx */ + smaller = magicints[MYMAX(FIRSTIDX, smallidx-1)] / 2; + small = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + larger = magicints[maxidx] / 2; + i = 0; + while (i < *size) { + is_small = 0; + thiscoord = (int *)(luip) + i * 3; + if (smallidx < maxidx && i >= 1 && + abs(thiscoord[0] - prevcoord[0]) < larger && + abs(thiscoord[1] - prevcoord[1]) < larger && + abs(thiscoord[2] - prevcoord[2]) < larger) { + is_smaller = 1; + } else if (smallidx > minidx) { + is_smaller = -1; + } else { + is_smaller = 0; + } + if (i + 1 < *size) { + if (abs(thiscoord[0] - thiscoord[3]) < small && + abs(thiscoord[1] - thiscoord[4]) < small && + abs(thiscoord[2] - thiscoord[5]) < small) { + /* interchange first with second atom for better + * compression of water molecules + */ + tmp = thiscoord[0]; thiscoord[0] = thiscoord[3]; + thiscoord[3] = tmp; + tmp = thiscoord[1]; thiscoord[1] = thiscoord[4]; + thiscoord[4] = tmp; + tmp = thiscoord[2]; thiscoord[2] = thiscoord[5]; + thiscoord[5] = tmp; + is_small = 1; + } + + } + tmpcoord[0] = thiscoord[0] - minint[0]; + tmpcoord[1] = thiscoord[1] - minint[1]; + tmpcoord[2] = thiscoord[2] - minint[2]; + if (bitsize == 0) { + sendbits(buf, bitsizeint[0], tmpcoord[0]); + sendbits(buf, bitsizeint[1], tmpcoord[1]); + sendbits(buf, bitsizeint[2], tmpcoord[2]); + } else { + sendints(buf, 3, bitsize, sizeint, tmpcoord); + } + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + thiscoord = thiscoord + 3; + i++; + + run = 0; + if (is_small == 0 && is_smaller == -1) + is_smaller = 0; + while (is_small && run < 8*3) { + if (is_smaller == -1 && (SQR(thiscoord[0] - prevcoord[0]) + + SQR(thiscoord[1] - prevcoord[1]) + + SQR(thiscoord[2] - prevcoord[2]) >= + smaller * smaller)) { + is_smaller = 0; + } + + tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small; + tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small; + tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small; + + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + + i++; + thiscoord = thiscoord + 3; + is_small = 0; + if (i < *size && + abs(thiscoord[0] - prevcoord[0]) < small && + abs(thiscoord[1] - prevcoord[1]) < small && + abs(thiscoord[2] - prevcoord[2]) < small) { + is_small = 1; + } + } + if (run != prevrun || is_smaller != 0) { + prevrun = run; + sendbits(buf, 1, 1); /* flag the change in run-length */ + sendbits(buf, 5, run+is_smaller+1); + } else { + sendbits(buf, 1, 0); /* flag the fact that runlength did not change */ + } + for (k=0; k < run; k+=3) { + sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]); + } + if (is_smaller != 0) { + smallidx += is_smaller; + if (is_smaller < 0) { + small = smaller; + smaller = magicints[smallidx-1] / 2; + } else { + smaller = small; + small = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + } + } + if (buf[1] != 0) buf[0]++;; + xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */ + return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0])); + } else { + + /* xdrs is open for reading */ + + if (xdr_int(xdrs, &lsize) == 0) + return 0; + if (*size != 0 && lsize != *size) { + fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; " + "%d arg vs %d in file", *size, lsize); + } + *size = lsize; + size3 = *size * 3; + if (*size <= 9) { + return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), + (xdrproc_t)xdr_float)); + } + xdr_float(xdrs, precision); + if (ip == NULL) { + ip = (int *) malloc(size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = (int) (size3 * 1.2); + buf = (int *) malloc(bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } else if (*size > oldsize) { + ip = (int *)realloc(ip, size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = (int) (size3 * 1.2); + buf = (int *)realloc(buf, bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } + buf[0] = buf[1] = buf[2] = 0; + + xdr_int(xdrs, &(minint[0])); + xdr_int(xdrs, &(minint[1])); + xdr_int(xdrs, &(minint[2])); + + xdr_int(xdrs, &(maxint[0])); + xdr_int(xdrs, &(maxint[1])); + xdr_int(xdrs, &(maxint[2])); + + sizeint[0] = maxint[0] - minint[0]+1; + sizeint[1] = maxint[1] - minint[1]+1; + sizeint[2] = maxint[2] - minint[2]+1; + + /* check if one of the sizes is to big to be multiplied */ + if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { + bitsizeint[0] = sizeofint(sizeint[0]); + bitsizeint[1] = sizeofint(sizeint[1]); + bitsizeint[2] = sizeofint(sizeint[2]); + bitsize = 0; /* flag the use of large sizes */ + } else { + bitsize = sizeofints(3, sizeint); + } + + xdr_int(xdrs, &smallidx); + maxidx = MYMIN(LASTIDX, smallidx + 8) ; + minidx = maxidx - 8; /* often this equal smallidx */ + smaller = magicints[MYMAX(FIRSTIDX, smallidx-1)] / 2; + small = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + larger = magicints[maxidx]; + + /* buf[0] holds the length in bytes */ + + if (xdr_int(xdrs, &(buf[0])) == 0) + return 0; + if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0) + return 0; + buf[0] = buf[1] = buf[2] = 0; + + lfp = fp; + inv_precision = 1.0 / * precision; + run = 0; + i = 0; + lip = ip; + while ( i < lsize ) { + thiscoord = (int *)(lip) + i * 3; + + if (bitsize == 0) { + thiscoord[0] = receivebits(buf, bitsizeint[0]); + thiscoord[1] = receivebits(buf, bitsizeint[1]); + thiscoord[2] = receivebits(buf, bitsizeint[2]); + } else { + receiveints(buf, 3, bitsize, sizeint, thiscoord); + } + + i++; + thiscoord[0] += minint[0]; + thiscoord[1] += minint[1]; + thiscoord[2] += minint[2]; + + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + + + flag = receivebits(buf, 1); + is_smaller = 0; + if (flag == 1) { + run = receivebits(buf, 5); + is_smaller = run % 3; + run -= is_smaller; + is_smaller--; + } + if (run > 0) { + thiscoord += 3; + for (k = 0; k < run; k+=3) { + receiveints(buf, 3, smallidx, sizesmall, thiscoord); + i++; + thiscoord[0] += prevcoord[0] - small; + thiscoord[1] += prevcoord[1] - small; + thiscoord[2] += prevcoord[2] - small; + if (k == 0) { + /* interchange first with second atom for better + * compression of water molecules + */ + tmp = thiscoord[0]; thiscoord[0] = prevcoord[0]; + prevcoord[0] = tmp; + tmp = thiscoord[1]; thiscoord[1] = prevcoord[1]; + prevcoord[1] = tmp; + tmp = thiscoord[2]; thiscoord[2] = prevcoord[2]; + prevcoord[2] = tmp; + *lfp++ = prevcoord[0] * inv_precision; + *lfp++ = prevcoord[1] * inv_precision; + *lfp++ = prevcoord[2] * inv_precision; + } else { + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + } + *lfp++ = thiscoord[0] * inv_precision; + *lfp++ = thiscoord[1] * inv_precision; + *lfp++ = thiscoord[2] * inv_precision; + } + } else { + *lfp++ = thiscoord[0] * inv_precision; + *lfp++ = thiscoord[1] * inv_precision; + *lfp++ = thiscoord[2] * inv_precision; + } + smallidx += is_smaller; + if (is_smaller < 0) { + small = smaller; + if (smallidx > FIRSTIDX) { + smaller = magicints[smallidx - 1] /2; + } else { + smaller = 0; + } + } else if (is_smaller > 0) { + smaller = small; + small = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + } + } + return 1; +} diff --git a/src/MISC/dump_xtc.h b/src/MISC/dump_xtc.h new file mode 100644 index 0000000000..9be66ca21d --- /dev/null +++ b/src/MISC/dump_xtc.h @@ -0,0 +1,118 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DUMP_CLASS + +DumpStyle(xtc,DumpXTC) + +#else + +#ifndef LMP_DUMP_XTC_H +#define LMP_DUMP_XTC_H + +#include "dump.h" + +#ifdef LAMMPS_XDR +#include "xdr_compat.h" +#else +#include "rpc/rpc.h" +#include "rpc/xdr.h" +#endif + +namespace LAMMPS_NS { + +class DumpXTC : public Dump { + public: + DumpXTC(class LAMMPS *, int, char**); + virtual ~DumpXTC(); + + private: + int natoms,ntotal; + int nevery_save; + int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no + float precision; // user-adjustable precision setting + float *coords; + double sfactor,tfactor; // scaling factors for positions and time unit + XDR xd; + + void init_style(); + int modify_param(int, char **); + void openfile(); + void write_header(bigint); + void pack(tagint *); + void write_data(int, double *); + bigint memory_usage(); + + void write_frame(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Invalid dump xtc filename + +Filenames used with the dump xtc style cannot be binary or compressed +or cause multiple files to be written. + +E: Too many atoms for dump xtc + +The system size must fit in a 32-bit integer to use this dump +style. + +W: No automatic unit conversion to XTC file format conventions possible for units lj + +This means no scaling will be performed. + +E: Dump xtc requires sorting by atom ID + +Use the dump_modify sort command to enable this. + +E: Cannot set dump_modify flush for dump xtc + +Self-explanatory. + +E: Cannot use variable every setting for dump xtc + +The format of this file requires snapshots at regular intervals. + +E: Cannot change dump_modify every for dump xtc + +The frequency of writing dump xtc snapshots cannot be changed. + +E: Cannot open dump file + +Self-explanatory. + +E: Too big a timestep for dump xtc + +The timestep must fit in a 32-bit integer to use this dump style. + +E: Illegal dump_modify sfactor value (must be > 0.0) + +Self-explanatory. + +E: Illegal dump_modify tfactor value (must be > 0.0) + +Self-explanatory. + +*/ diff --git a/src/MISC/xdr_compat.cpp b/src/MISC/xdr_compat.cpp new file mode 100644 index 0000000000..12bc777b90 --- /dev/null +++ b/src/MISC/xdr_compat.cpp @@ -0,0 +1,722 @@ +#ifdef LAMMPS_XDR + +#include +#include +#include +#include "xdr_compat.h" + +/* This file is needed for systems, that do not provide XDR support + * in their system libraries. It was written for windows, but will + * most probably work on other platforms too. better make sure you + * test that the xtc files produced are ok before using it. + * + * It is also needed on BG/L and Cray XT3/XT4 as we don't have + * XDR support in the lightweight kernel runtimes either. + * + * This file contains the definitions for Sun External Data + * Representation (XDR) headers and routines. + * + * Although the rest of LAMPPS is GPL, you can copy and use the XDR + * routines in any way you want as long as you obey Sun's license: + * + * Sun RPC is a product of Sun Microsystems, Inc. and is provided for + * unrestricted use provided that this legend is included on all tape + * media and as a part of the software program in whole or part. Users + * may copy or modify Sun RPC without charge, but are not authorized + * to license or distribute it to anyone else except as part of a product or + * program developed by the user. + * + * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE + * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR + * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE. + * + * Sun RPC is provided with no support and without any obligation on the + * part of Sun Microsystems, Inc. to assist in its use, correction, + * modification or enhancement. + * + * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE + * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC + * OR ANY PART THEREOF. + * + * In no event will Sun Microsystems, Inc. be liable for any lost revenue + * or profits or other special, indirect and consequential damages, even if + * Sun has been advised of the possibility of such damages. + * + * Sun Microsystems, Inc. + * 2550 Garcia Avenue + * Mountain View, California 94043 + */ + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * for unit alignment + */ +static char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0}; + +static xdr_uint32_t xdr_swapbytes(xdr_uint32_t x) +{ + xdr_uint32_t y; + int i; + char *px=(char *)&x; + char *py=(char *)&y; + + for(i=0;i<4;i++) + py[i]=px[3-i]; + + return y; +} + +static xdr_uint32_t xdr_htonl(xdr_uint32_t x) +{ + short s=0x0F00; + if( *((char *)&s)==(char)0x0F) { + /* bigendian, do nothing */ + return x; + } else { + /* smallendian,swap bytes */ + return xdr_swapbytes(x); + } +} + +static xdr_uint32_t xdr_ntohl(xdr_uint32_t x) +{ + short s=0x0F00; + if( *((char *)&s)==(char)0x0F) { + /* bigendian, do nothing */ + return x; + } else { + /* smallendian, swap bytes */ + return xdr_swapbytes(x); + } +} + + +/* + * Free a data structure using XDR + * Not a filter, but a convenient utility nonetheless + */ +void +xdr_free (xdrproc_t proc, char *objp) +{ + XDR x; + + x.x_op = XDR_FREE; + (*proc) (&x, objp); +} + +/* + * XDR nothing + */ +bool_t +xdr_void (void) +{ + return TRUE; +} + +/* + * XDR integers + */ +bool_t +xdr_int (XDR *xdrs, int *ip) +{ + xdr_int32_t l; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + l = (xdr_int32_t) (*ip); + return xdr_putint32 (xdrs, &l); + + case XDR_DECODE: + if (!xdr_getint32 (xdrs, &l)) + { + return FALSE; + } + *ip = (int) l; + + case XDR_FREE: + return TRUE; + } + return FALSE; +} + + +/* + * XDR unsigned integers + */ +bool_t +xdr_u_int (XDR *xdrs, unsigned int *up) +{ + xdr_uint32_t l; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + l = (xdr_uint32_t) (*up); + return xdr_putuint32 (xdrs, &l); + + case XDR_DECODE: + if (!xdr_getuint32 (xdrs, &l)) + { + return FALSE; + } + *up = (unsigned int) l; + + case XDR_FREE: + return TRUE; + } + return FALSE; +} + + + + +/* + * XDR short integers + */ +bool_t +xdr_short (XDR *xdrs, short *sp) +{ + xdr_int32_t l; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + l = (xdr_int32_t) *sp; + return xdr_putint32 (xdrs, &l); + + case XDR_DECODE: + if (!xdr_getint32 (xdrs, &l)) + { + return FALSE; + } + *sp = (short) l; + return TRUE; + + case XDR_FREE: + return TRUE; + } + return FALSE; +} + + +/* + * XDR unsigned short integers + */ +bool_t +xdr_u_short (XDR *xdrs, unsigned short *usp) +{ + xdr_uint32_t l; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + l = (xdr_uint32_t) *usp; + return xdr_putuint32 (xdrs, &l); + + case XDR_DECODE: + if (!xdr_getuint32 (xdrs, &l)) + { + return FALSE; + } + *usp = (unsigned short) l; + return TRUE; + + case XDR_FREE: + return TRUE; + } + return FALSE; +} + + +/* + * XDR a char + */ +bool_t +xdr_char (XDR *xdrs, char *cp) +{ + int i; + + i = (*cp); + if (!xdr_int (xdrs, &i)) + { + return FALSE; + } + *cp = i; + return TRUE; +} + +/* + * XDR an unsigned char + */ +bool_t +xdr_u_char (XDR *xdrs, unsigned char *cp) +{ + unsigned int u; + + u = (*cp); + if (!xdr_u_int (xdrs, &u)) + { + return FALSE; + } + *cp = u; + return TRUE; +} + +/* + * XDR booleans + */ +bool_t +xdr_bool (XDR *xdrs, int *bp) +{ +#define XDR_FALSE ((xdr_int32_t) 0) +#define XDR_TRUE ((xdr_int32_t) 1) + + xdr_int32_t lb; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + lb = *bp ? XDR_TRUE : XDR_FALSE; + return xdr_putint32 (xdrs, &lb); + + case XDR_DECODE: + if (!xdr_getint32 (xdrs, &lb)) + { + return FALSE; + } + *bp = (lb == XDR_FALSE) ? FALSE : TRUE; + return TRUE; + + case XDR_FREE: + return TRUE; + } + return FALSE; +#undef XDR_FALSE +#undef XDR_TRUE +} + + + +/* + * XDR opaque data + * Allows the specification of a fixed size sequence of opaque bytes. + * cp points to the opaque object and cnt gives the byte length. + */ +bool_t +xdr_opaque (XDR *xdrs, char *cp, unsigned int cnt) +{ + unsigned int rndup; + static char crud[BYTES_PER_XDR_UNIT]; + + /* + * if no data we are done + */ + if (cnt == 0) + return TRUE; + + /* + * round byte count to full xdr units + */ + rndup = cnt % BYTES_PER_XDR_UNIT; + if (rndup > 0) + rndup = BYTES_PER_XDR_UNIT - rndup; + + switch (xdrs->x_op) + { + case XDR_DECODE: + if (!xdr_getbytes (xdrs, cp, cnt)) + { + return FALSE; + } + if (rndup == 0) + return TRUE; + return xdr_getbytes (xdrs, (char *)crud, rndup); + + case XDR_ENCODE: + if (!xdr_putbytes (xdrs, cp, cnt)) + { + return FALSE; + } + if (rndup == 0) + return TRUE; + return xdr_putbytes (xdrs, xdr_zero, rndup); + + case XDR_FREE: + return TRUE; + } + return FALSE; +} + + +/* + * XDR null terminated ASCII strings + * xdr_string deals with "C strings" - arrays of bytes that are + * terminated by a NULL character. The parameter cpp references a + * pointer to storage; If the pointer is null, then the necessary + * storage is allocated. The last parameter is the max allowed length + * of the string as specified by a protocol. + */ +bool_t +xdr_string (XDR *xdrs, char **cpp, unsigned int maxsize) +{ + char *sp = *cpp; /* sp is the actual string pointer */ + unsigned int size = 0; + unsigned int nodesize = 0; + + /* + * first deal with the length since xdr strings are counted-strings + */ + switch (xdrs->x_op) + { + case XDR_FREE: + if (sp == NULL) + { + return TRUE; /* already free */ + } + /* fall through... */ + case XDR_ENCODE: + if (sp == NULL) + return FALSE; + size = strlen (sp); + break; + case XDR_DECODE: + break; + } + + if (!xdr_u_int (xdrs, &size)) + { + return FALSE; + } + if (size > maxsize) + { + return FALSE; + } + nodesize = size + 1; + + /* + * now deal with the actual bytes + */ + switch (xdrs->x_op) + { + case XDR_DECODE: + if (nodesize == 0) + { + return TRUE; + } + if (sp == NULL) + *cpp = sp = (char *) malloc (nodesize); + if (sp == NULL) + { + (void) fputs ("xdr_string: out of memory\n", stderr); + return FALSE; + } + sp[size] = 0; + /* fall into ... */ + + case XDR_ENCODE: + return xdr_opaque (xdrs, sp, size); + + case XDR_FREE: + free (sp); + *cpp = NULL; + return TRUE; + } + return FALSE; +} + + + +/* Floating-point stuff */ + +bool_t +xdr_float(XDR *xdrs, float *fp) +{ + xdr_int32_t tmp; + + switch (xdrs->x_op) { + + case XDR_ENCODE: + tmp = *(xdr_int32_t *)fp; + return (xdr_putint32(xdrs, &tmp)); + + break; + + case XDR_DECODE: + if (xdr_getint32(xdrs, &tmp)) { + *(xdr_int32_t *)fp = tmp; + return (TRUE); + } + + break; + + case XDR_FREE: + return (TRUE); + } + return (FALSE); +} + + +bool_t +xdr_double(XDR *xdrs, double *dp) +{ + + /* Windows and some other systems dont define double-precision + * word order in the header files, so unfortunately we have + * to calculate it! + */ + static int LSW=-1; /* Least significant fp word */ + int *ip; + xdr_int32_t tmp[2]; + + if(LSW<0) { + double x=0.987654321; /* Just a number */ + + /* Possible representations in IEEE double precision: + * (S=small endian, B=big endian) + * + * Byte order, Word order, Hex + * S S b8 56 0e 3c dd 9a ef 3f + * B S 3c 0e 56 b8 3f ef 9a dd + * S B dd 9a ef 3f b8 56 0e 3c + * B B 3f ef 9a dd 3c 0e 56 b8 + */ + + unsigned char ix = *((char *)&x); + + if(ix==0xdd || ix==0x3f) + LSW=1; /* Big endian word order */ + else if(ix==0xb8 || ix==0x3c) + LSW=0; /* Small endian word order */ + else { /* Catch strange errors */ + printf("Error when detecting floating-point word order.\n" + "Do you have a non-IEEE system?\n" + "If possible, use the XDR libraries provided with your system,\n" + "instead of the Gromacs fallback XDR source.\n"); + exit(0); + } + } + + switch (xdrs->x_op) { + + case XDR_ENCODE: + ip = (int *)dp; + tmp[0] = ip[!LSW]; + tmp[1] = ip[LSW]; + return (xdr_putint32(xdrs, tmp) && + xdr_putint32(xdrs, tmp+1)); + + break; + + case XDR_DECODE: + ip = (int *)dp; + if (xdr_getint32(xdrs, tmp+!LSW) && + xdr_getint32(xdrs, tmp+LSW)) { + ip[0] = tmp[0]; + ip[1] = tmp[1]; + return (TRUE); + } + + break; + + case XDR_FREE: + return (TRUE); + } + return (FALSE); +} + + +/* Array routines */ + +/* + * xdr_vector(): + * + * XDR a fixed length array. Unlike variable-length arrays, + * the storage of fixed length arrays is static and unfreeable. + * > basep: base of the array + * > size: size of the array + * > elemsize: size of each element + * > xdr_elem: routine to XDR each element + */ +bool_t +xdr_vector (XDR *xdrs, char *basep, unsigned int nelem, + unsigned int elemsize, xdrproc_t xdr_elem) +{ +#define LASTUNSIGNED ((unsigned int)0-1) + unsigned int i; + char *elptr; + + elptr = basep; + for (i = 0; i < nelem; i++) + { + if (!(*xdr_elem) (xdrs, elptr, LASTUNSIGNED)) + { + return FALSE; + } + elptr += elemsize; + } + return TRUE; +#undef LASTUNSIGNED +} + + + +static bool_t xdrstdio_getbytes (XDR *, char *, unsigned int); +static bool_t xdrstdio_putbytes (XDR *, char *, unsigned int); +static unsigned int xdrstdio_getpos (XDR *); +static bool_t xdrstdio_setpos (XDR *, unsigned int); +static xdr_int32_t *xdrstdio_inline (XDR *, int); +static void xdrstdio_destroy (XDR *); +static bool_t xdrstdio_getint32 (XDR *, xdr_int32_t *); +static bool_t xdrstdio_putint32 (XDR *, xdr_int32_t *); +static bool_t xdrstdio_getuint32 (XDR *, xdr_uint32_t *); +static bool_t xdrstdio_putuint32 (XDR *, xdr_uint32_t *); + +/* + * Ops vector for stdio type XDR + */ +static const struct xdr_ops xdrstdio_ops = +{ + xdrstdio_getbytes, /* deserialize counted bytes */ + xdrstdio_putbytes, /* serialize counted bytes */ + xdrstdio_getpos, /* get offset in the stream */ + xdrstdio_setpos, /* set offset in the stream */ + xdrstdio_inline, /* prime stream for inline macros */ + xdrstdio_destroy, /* destroy stream */ + xdrstdio_getint32, /* deserialize a int */ + xdrstdio_putint32, /* serialize a int */ + xdrstdio_getuint32, /* deserialize a int */ + xdrstdio_putuint32 /* serialize a int */ +}; + +/* + * Initialize a stdio xdr stream. + * Sets the xdr stream handle xdrs for use on the stream file. + * Operation flag is set to op. + */ +void +xdrstdio_create (XDR *xdrs, FILE *file, enum xdr_op op) +{ + xdrs->x_op = op; + /* We have to add the const since the `struct xdr_ops' in `struct XDR' + is not `const'. */ + xdrs->x_ops = (struct xdr_ops *) &xdrstdio_ops; + xdrs->x_private = (char *) file; + xdrs->x_handy = 0; + xdrs->x_base = 0; +} + +/* + * Destroy a stdio xdr stream. + * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create. + */ +static void +xdrstdio_destroy (XDR *xdrs) +{ + (void) fflush ((FILE *) xdrs->x_private); + /* xx should we close the file ?? */ +} + + +static bool_t +xdrstdio_getbytes (XDR *xdrs, char *addr, unsigned int len) +{ + if ((len != 0) && (fread (addr, (int) len, 1, + (FILE *) xdrs->x_private) != 1)) + return FALSE; + return TRUE; +} + +static bool_t +xdrstdio_putbytes (XDR *xdrs, char *addr, unsigned int len) +{ + if ((len != 0) && (fwrite (addr, (int) len, 1, + (FILE *) xdrs->x_private) != 1)) + return FALSE; + return TRUE; +} + +static unsigned int +xdrstdio_getpos (XDR *xdrs) +{ + return (unsigned int) ftell ((FILE *) xdrs->x_private); +} + +static bool_t +xdrstdio_setpos (XDR *xdrs, unsigned int pos) +{ + return fseek ((FILE *) xdrs->x_private, (xdr_int32_t) pos, 0) < 0 ? FALSE : TRUE; +} + +static xdr_int32_t * +xdrstdio_inline (XDR *xdrs, int len) +{ + /* + * Must do some work to implement this: must insure + * enough data in the underlying stdio buffer, + * that the buffer is aligned so that we can indirect through a + * long *, and stuff this pointer in xdrs->x_buf. Doing + * a fread or fwrite to a scratch buffer would defeat + * most of the gains to be had here and require storage + * management on this buffer, so we don't do this. + */ + return NULL; +} + +static bool_t +xdrstdio_getint32 (XDR *xdrs, xdr_int32_t *ip) +{ + xdr_int32_t mycopy; + + if (fread ((char *) &mycopy, 4, 1, (FILE *) xdrs->x_private) != 1) + return FALSE; + *ip = xdr_ntohl (mycopy); + return TRUE; +} + +static bool_t +xdrstdio_putint32 (XDR *xdrs, xdr_int32_t *ip) +{ + xdr_int32_t mycopy = xdr_htonl (*ip); + + ip = &mycopy; + if (fwrite ((char *) ip, 4, 1, (FILE *) xdrs->x_private) != 1) + return FALSE; + return TRUE; +} + +static bool_t +xdrstdio_getuint32 (XDR *xdrs, xdr_uint32_t *ip) +{ + xdr_uint32_t mycopy; + + if (fread ((char *) &mycopy, 4, 1, (FILE *) xdrs->x_private) != 1) + return FALSE; + *ip = xdr_ntohl (mycopy); + return TRUE; +} + +static bool_t +xdrstdio_putuint32 (XDR *xdrs, xdr_uint32_t *ip) +{ + xdr_uint32_t mycopy = xdr_htonl (*ip); + + ip = &mycopy; + if (fwrite ((char *) ip, 4, 1, (FILE *) xdrs->x_private) != 1) + return FALSE; + return TRUE; +} + +#ifdef __cplusplus +} +#endif + +#else +/* satisfy compilers that do not like to compile empty files. */ +static void i_am_a_dummy_subroutine(void) { + return; +} +#endif diff --git a/src/MISC/xdr_compat.h b/src/MISC/xdr_compat.h new file mode 100644 index 0000000000..ac318aa58c --- /dev/null +++ b/src/MISC/xdr_compat.h @@ -0,0 +1,228 @@ +#ifndef LMP_XDR_COMPAT_H +#define LMP_XDR_COMPAT_H + +#include +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * This file is needed for systems, that do not provide XDR support + * in their system libraries. It was written for windows, but will + * most probably work on other platforms too. better make sure you + * test that the xtc files produced are ok before using it. + * + * It is also needed on BG/L, BG/P and Cray XT3/XT4/XT5 as we don't + * have XDR support in the lightweight kernel runtimes either. + * + * This file contains the definitions for Sun External Data + * Representation (XDR) headers and routines. + * + * Although the rest of LAMPPS is GPL, you can copy and use the XDR + * routines in any way you want as long as you obey Sun's license: + * Sun RPC is a product of Sun Microsystems, Inc. and is provided for + * unrestricted use provided that this legend is included on all tape + * media and as a part of the software program in whole or part. Users + * may copy or modify Sun RPC without charge, but are not authorized + * to license or distribute it to anyone else except as part of a product or + * program developed by the user. + * + * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE + * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR + * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE. + * + * Sun RPC is provided with no support and without any obligation on the + * part of Sun Microsystems, Inc. to assist in its use, correction, + * modification or enhancement. + * + * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE + * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC + * OR ANY PART THEREOF. + * + * In no event will Sun Microsystems, Inc. be liable for any lost revenue + * or profits or other special, indirect and consequential damages, even if + * Sun has been advised of the possibility of such damages. + * + * Sun Microsystems, Inc. + * 2550 Garcia Avenue + * Mountain View, California 94043 + */ + +/* + * Xdr operations. XDR_ENCODE causes the type to be encoded into the + * stream. XDR_DECODE causes the type to be extracted from the stream. + * XDR_FREE can be used to release the space allocated by an + * XDR_DECODE request. + */ + +typedef int bool_t; + +#if defined(__MINGW32__) +typedef char * caddr_t; +typedef unsigned int u_int; +#endif + +/* + * Aninteger type that is 32 bits wide. Check if int, + * long or short is 32 bits and die if none of them is :-) + */ +#if (INT_MAX == 2147483647) + typedef int xdr_int32_t; + typedef unsigned int xdr_uint32_t; +#elif (LONG_MAX == 2147483647L) + typedef long xdr_int32_t; + typedef unsigned long xdr_uint32_t; +#elif (SHRT_MAX == 2147483647) + typedef short xdr_int32_t; + typedef unsigned short xdr_uint32_t; +#else +# error ERROR: No 32 bit wide integer type found! +#endif + +enum xdr_op { + XDR_ENCODE = 0, + XDR_DECODE = 1, + XDR_FREE = 2 +}; + +#ifndef FALSE +# define FALSE (0) +#endif +#ifndef TRUE +# define TRUE (1) +#endif + +#define BYTES_PER_XDR_UNIT (4) +/* Macro to round up to units of 4. */ +#define XDR_RNDUP(x) (((x) + BYTES_PER_XDR_UNIT - 1) & ~(BYTES_PER_XDR_UNIT - 1)) + + +/* + * The XDR handle. + * Contains operation which is being applied to the stream, + * an operations vector for the particular implementation (e.g. see xdr_mem.c), + * and two private fields for the use of the particular implementation. + */ + +typedef struct XDR XDR; +struct XDR + { + enum xdr_op x_op; /* operation; fast additional param */ + struct xdr_ops *x_ops; + char *x_public; /* users' data */ + char *x_private; /* pointer to private data */ + char *x_base; /* private used for position info */ + int x_handy; /* extra private word */ + }; + +struct xdr_ops + { + bool_t (*x_getbytes) (XDR *__xdrs, char *__addr, unsigned int __len); + /* get some bytes from " */ + bool_t (*x_putbytes) (XDR *__xdrs, char *__addr, unsigned int __len); + /* put some bytes to " */ + unsigned int (*x_getpostn) (XDR *__xdrs); + /* returns bytes off from beginning */ + bool_t (*x_setpostn) (XDR *__xdrs, unsigned int __pos); + /* lets you reposition the stream */ + xdr_int32_t *(*x_inline) (XDR *__xdrs, int __len); + /* buf quick ptr to buffered data */ + void (*x_destroy) (XDR *__xdrs); + /* free privates of this xdr_stream */ + bool_t (*x_getint32) (XDR *__xdrs, xdr_int32_t *__ip); + /* get a int from underlying stream */ + bool_t (*x_putint32) (XDR *__xdrs, xdr_int32_t *__ip); + /* put a int to " */ + bool_t (*x_getuint32) (XDR *__xdrs, xdr_uint32_t *__ip); + /* get a unsigned int from underlying stream */ + bool_t (*x_putuint32) (XDR *__xdrs, xdr_uint32_t *__ip); + /* put a int to " */ +}; + +/* + * A xdrproc_t exists for each data type which is to be encoded or decoded. + * + * The second argument to the xdrproc_t is a pointer to an opaque pointer. + * The opaque pointer generally points to a structure of the data type + * to be decoded. If this pointer is 0, then the type routines should + * allocate dynamic storage of the appropriate size and return it. + */ + +typedef bool_t (*xdrproc_t) (XDR *, void *,...); + +/* + * Operations defined on a XDR handle + * + * XDR *xdrs; + * xdr_int32_t *int32p; + * long *longp; + * char *addr; + * unsigned int len; + * unsigned int pos; + */ + + +#define xdr_getint32(xdrs, int32p) \ + (*(xdrs)->x_ops->x_getint32)(xdrs, int32p) + +#define xdr_putint32(xdrs, int32p) \ + (*(xdrs)->x_ops->x_putint32)(xdrs, int32p) + +#define xdr_getuint32(xdrs, uint32p) \ + (*(xdrs)->x_ops->x_getuint32)(xdrs, uint32p) + +#define xdr_putuint32(xdrs, uint32p) \ + (*(xdrs)->x_ops->x_putuint32)(xdrs, uint32p) + +#define xdr_getbytes(xdrs, addr, len) \ + (*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len) + +#define xdr_putbytes(xdrs, addr, len) \ + (*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len) + +#define xdr_getpos(xdrs) \ + (*(xdrs)->x_ops->x_getpostn)(xdrs) + +#define xdr_setpos(xdrs, pos) \ + (*(xdrs)->x_ops->x_setpostn)(xdrs, pos) + +#define xdr_inline(xdrs, len) \ + (*(xdrs)->x_ops->x_inline)(xdrs, len) + +#define xdr_destroy(xdrs) \ + do { \ + if ((xdrs)->x_ops->x_destroy) \ + (*(xdrs)->x_ops->x_destroy)(xdrs); \ + } while (0) + + +extern bool_t xdr_int (XDR *__xdrs, int *__ip); +extern bool_t xdr_u_int (XDR *__xdrs, unsigned int *__ip); +extern bool_t xdr_short (XDR *__xdrs, short *__ip); +extern bool_t xdr_u_short (XDR *__xdrs, unsigned short *__ip); +extern bool_t xdr_bool (XDR *__xdrs, int *__bp); +extern bool_t xdr_opaque (XDR *__xdrs, char *__cp, unsigned int __cnt); +extern bool_t xdr_string (XDR *__xdrs, char **__cpp, unsigned int __maxsize); +extern bool_t xdr_char (XDR *__xdrs, char *__cp); +extern bool_t xdr_u_char (XDR *__xdrs, unsigned char *__cp); +extern bool_t xdr_vector (XDR *__xdrs, char *__basep, unsigned int __nelem, + unsigned int __elemsize, xdrproc_t __xdr_elem); +extern bool_t xdr_float (XDR *__xdrs, float *__fp); +extern bool_t xdr_double (XDR *__xdrs, double *__dp); +extern void xdrstdio_create (XDR *__xdrs, FILE *__file, enum xdr_op __xop); + +/* free memory buffers for xdr */ +extern void xdr_free (xdrproc_t __proc, char *__objp); + +#ifdef __cplusplus +} +#endif + +#endif /* XDR_COMPAT_H */ + +/* ERROR/WARNING messages: + +*/