Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78835216
dump_xtc.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 23, 06:07
Size
36 KB
Mime Type
text/x-c
Expires
Sun, Aug 25, 06:07 (2 d)
Engine
blob
Format
Raw Data
Handle
19617207
Attached To
rLAMMPS lammps
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;
}
Event Timeline
Log In to Comment