Homec4science
Diffusion lammps (master) Browse

lammps/src/USER-MISC/dump_xtc.cpp
0d4d950ce0c2master

Last Change
Raw File

dump_xtc.cpp
View Options

/* ----------------------------------------------------------------------
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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#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<bigint>(MAXSMALLINT/3/sizeof(float)))
error->all(FLERR,"Too many atoms for dump xtc");
natoms = static_cast<int> (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;
}
c4science · Help