From fe66246c6bb6957c7ecb058e0ac23806e3d43059 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 17 Jan 2008 20:22:20 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1380 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/XTC/Install.csh | 6 + src/XTC/dump_xtc.cpp | 81 +++-- src/XTC/dump_xtc.h | 9 +- src/XTC/xdr_compat.cpp | 717 +++++++++++++++++++++++++++++++++++++++ src/XTC/xdr_compat.h | 222 ++++++++++++ src/dump_dcd.cpp | 91 +++-- src/dump_dcd.h | 3 +- src/style_user_ackland.h | 20 ++ src/style_user_ewaldn.h | 30 ++ 9 files changed, 1126 insertions(+), 53 deletions(-) create mode 100644 src/XTC/xdr_compat.cpp create mode 100644 src/XTC/xdr_compat.h diff --git a/src/XTC/Install.csh b/src/XTC/Install.csh index 319ec056ed..bb1573f031 100644 --- a/src/XTC/Install.csh +++ b/src/XTC/Install.csh @@ -8,6 +8,9 @@ if ($1 == 1) then cp dump_xtc.h .. + cp xdr_compat.cpp .. + cp xdr_compat.h .. + else if ($1 == 0) then rm ../style_xtc.h @@ -17,4 +20,7 @@ else if ($1 == 0) then rm ../dump_xtc.h + rm ../xdr_compat.cpp + rm ../xdr_compat.h + endif diff --git a/src/XTC/dump_xtc.cpp b/src/XTC/dump_xtc.cpp index f81f7382f8..9a66cb5901 100644 --- a/src/XTC/dump_xtc.cpp +++ b/src/XTC/dump_xtc.cpp @@ -23,8 +23,6 @@ #include "stdlib.h" #include "string.h" #include "limits.h" -#include "rpc/rpc.h" -#include "rpc/xdr.h" #include "dump_xtc.h" #include "domain.h" #include "atom.h" @@ -42,15 +40,21 @@ using namespace LAMMPS_NS; #define MYMAX(a,b) ((a) > (b) ? (a) : (b)) int xdropen(XDR *, const char *, const char *); -int xdrclose(XDR *) ; +int xdrclose(XDR *); void xdrfreebuf(); int xdr3dfcoord(XDR *, float *, int *, float *); +// include XDR compatibility code + +#ifdef LAMMPS_USE_XDR_COMPAT +#include "xdr_compat.c" +#endif + /* ---------------------------------------------------------------------- */ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { - if (narg != 5 && narg != 6) error->all("Illegal dump xtc command"); + if (narg != 5) error->all("Illegal dump xtc command"); if (igroup != group->find("all")) error->all("Dump xtc must use group all"); if (binary || compressed || multifile || multiproc) error->all("Invalid dump xtc filename"); @@ -58,17 +62,8 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) size_one = 4; format_default = NULL; flush_flag = 0; - - // parse extra argument - + unwrap_flag = 0; precision = 1000.0; - if (narg == 6) { - precision = atof(arg[5]); - 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("Illegal dump xtc command"); - } // allocate global array for atom coords @@ -79,6 +74,7 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); // sfactor = conversion of coords to XTC units + // GROMACS standard is nanometers, not Angstroms sfactor = 0.1; if (strcmp(update->unit_style,"lj") == 0) sfactor = 1.0; @@ -92,6 +88,7 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) DumpXTC::~DumpXTC() { memory->sfree(coords); + if (me == 0) { xdrclose(&xd); xdrfreebuf(); @@ -107,6 +104,29 @@ void DumpXTC::init() if (flush_flag) error->all("Cannot set dump_modify flush for dump xtc"); } +/* ---------------------------------------------------------------------- */ + +int DumpXTC::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all("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("Illegal dump_modify command"); + return 2; + } else if (strcmp(arg[0],"precision") == 0) { + if (narg < 2) error->all("Illegal dump_modify command"); + precision = atof(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("Illegal dump_modify command"); + return 2; + } + return 0; +} + /* ---------------------------------------------------------------------- return # of bytes of allocated memory in buf and global coords array ------------------------------------------------------------------------- */ @@ -154,9 +174,9 @@ void DumpXTC::write_header(int n) // cell basis vectors float zero = 0.0; - float xdim = domain->boxhi[0] - domain->boxlo[0]; - float ydim = domain->boxhi[1] - domain->boxlo[1]; - float zdim = domain->boxhi[2] - domain->boxlo[2]; + 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); @@ -176,16 +196,31 @@ int DumpXTC::pack() { int *tag = atom->tag; double **x = atom->x; + int *image = atom->image; int nlocal = atom->nlocal; // assume group all, so no need to perform mask check int m = 0; - for (int i = 0; i < nlocal; i++) { - buf[m++] = tag[i]; - buf[m++] = sfactor*x[i][0]; - buf[m++] = sfactor*x[i][1]; - buf[m++] = sfactor*x[i][2]; + if (unwrap_flag == 1) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + for (int i = 0; i < nlocal; i++) { + buf[m++] = tag[i]; + buf[m++] = sfactor*(x[i][0] + ((image[i] & 1023) - 512) * xprd); + buf[m++] = sfactor*(x[i][1] + ((image[i] >> 10 & 1023) - 512) * yprd); + buf[m++] = sfactor*(x[i][2] + ((image[i] >> 20) - 512) * zprd); + } + + } else { + for (int i = 0; i < nlocal; i++) { + buf[m++] = tag[i]; + buf[m++] = sfactor*x[i][0]; + buf[m++] = sfactor*x[i][1]; + buf[m++] = sfactor*x[i][2]; + } } return m; @@ -264,7 +299,6 @@ void DumpXTC::write_frame() static FILE *xdrfiles[MAXID]; static XDR *xdridptr[MAXID]; static char xdrmodes[MAXID]; -static unsigned int cnt; static int *ip = NULL; static int *buf = NULL; @@ -384,6 +418,7 @@ int xdrclose(XDR *xdrs) } fprintf(stderr, "xdrclose: no such open xdr file\n"); exit(1); + return 1; } /*_________________________________________________________________________ diff --git a/src/XTC/dump_xtc.h b/src/XTC/dump_xtc.h index 7e6e0571f6..eabbb2b3ff 100644 --- a/src/XTC/dump_xtc.h +++ b/src/XTC/dump_xtc.h @@ -15,8 +15,13 @@ #define 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 { @@ -29,11 +34,13 @@ class DumpXTC : public Dump { private: int natoms,ntotal; - float precision; + int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no + float precision; // user-adjustable precision setting float *coords; double sfactor; XDR xd; + int modify_param(int, char **); void openfile(); void write_header(int); int count(); diff --git a/src/XTC/xdr_compat.cpp b/src/XTC/xdr_compat.cpp new file mode 100644 index 0000000000..866df9a9ea --- /dev/null +++ b/src/XTC/xdr_compat.cpp @@ -0,0 +1,717 @@ +#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 + +#endif diff --git a/src/XTC/xdr_compat.h b/src/XTC/xdr_compat.h new file mode 100644 index 0000000000..78863c2130 --- /dev/null +++ b/src/XTC/xdr_compat.h @@ -0,0 +1,222 @@ +#ifndef XDR_COMPAT_H +#define 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 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 + */ + +/* + * 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; + +/* + * 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 */ + diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index ac4faaa044..46570b415a 100644 --- a/src/dump_dcd.cpp +++ b/src/dump_dcd.cpp @@ -44,8 +44,14 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : - Dump(lmp, narg, arg) +static inline void fwrite_int32(FILE* fd, uint32_t i) +{ + fwrite(&i,sizeof(uint32_t),1,fd); +} + +/* ---------------------------------------------------------------------- */ + +DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { if (narg != 5) error->all("Illegal dump dcd command"); if (igroup != group->find("all")) error->all("Dump dcd must use group all"); @@ -53,8 +59,9 @@ DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : error->all("Invalid dump dcd filename"); size_one = 4; + unwrap_flag = 0; format_default = NULL; - + // allocate global array for atom coords natoms = static_cast (atom->natoms); @@ -83,6 +90,9 @@ DumpDCD::~DumpDCD() void DumpDCD::init() { + if (unwrap_flag == 1 && domain->triclinic) + error->all("Dump dcd cannot dump unwrapped coords with triclinic box"); + // check that dump frequency has not changed if (nevery_save == 0) { @@ -99,6 +109,20 @@ void DumpDCD::init() } } +/* ---------------------------------------------------------------------- */ + +int DumpDCD::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all("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("Illegal dump_modify command"); + return 2; + } + return 0; +} + /* ---------------------------------------------------------------------- return # of bytes of allocated memory in buf and global coords array ------------------------------------------------------------------------- */ @@ -156,9 +180,10 @@ void DumpDCD::write_header(int n) } uint32_t out_integer = 48; - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); fwrite(dim,out_integer,1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); + if (flush_flag) fflush(fp); } /* ---------------------------------------------------------------------- */ @@ -174,16 +199,31 @@ int DumpDCD::pack() { int *tag = atom->tag; double **x = atom->x; + int *image = atom->image; int nlocal = atom->nlocal; // assume group all, so no need to perform mask check int m = 0; - for (int i = 0; i < nlocal; i++) { - buf[m++] = tag[i]; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; + if (unwrap_flag) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + for (int i = 0; i < nlocal; i++) { + buf[m++] = tag[i]; + buf[m++] = x[i][0] + ((image[i] & 1023) - 512) * xprd; + buf[m++] = x[i][1] + ((image[i] >> 10 & 1023) - 512) * yprd; + buf[m++] = x[i][2] + ((image[i] >> 20) - 512) * zprd; + } + + } else { + for (int i = 0; i < nlocal; i++) { + buf[m++] = tag[i]; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + } } return m; @@ -221,25 +261,25 @@ void DumpDCD::write_frame() // write coords uint32_t out_integer = natoms*sizeof(float); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); fwrite(xf,out_integer,1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); + fwrite_int32(fp,out_integer); fwrite(yf,out_integer,1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); + fwrite_int32(fp,out_integer); fwrite(zf,out_integer,1,fp); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); // update NFILE and NSTEP fields in DCD header nframes++; out_integer = nframes; fseek(fp,NFILE_POS,SEEK_SET); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); out_integer = update->ntimestep; fseek(fp,NSTEP_POS,SEEK_SET); - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); fseek(fp,0,SEEK_END); } @@ -252,10 +292,9 @@ void DumpDCD::write_dcd_header(const char *remarks) char title_string[200]; time_t cur_time; struct tm *tmbuf; - char time_str[81]; out_integer = 84; - fwrite(&out_integer,sizeof(uint32_t),1,fp); + fwrite_int32(fp,out_integer); strcpy(title_string,"CORD"); fwrite(title_string,4,1,fp); fwrite_int32(fp,0); // NFILE = # of snapshots in file @@ -287,16 +326,12 @@ void DumpDCD::write_dcd_header(const char *remarks) fwrite(title_string,80,1,fp); cur_time=time(NULL); tmbuf=localtime(&cur_time); - strftime(time_str,80,"REMARKS Created %d %B,%Y at %R",tmbuf); - fwrite(time_str,80,1,fp); + memset(title_string,' ',81); + strftime(title_string,80,"REMARKS Created %d %B,%Y at %R",tmbuf); + fwrite(title_string,80,1,fp); fwrite_int32(fp,164); fwrite_int32(fp,4); fwrite_int32(fp,natoms); // number of atoms in each snapshot fwrite_int32(fp,4); -} - -/* ---------------------------------------------------------------------- */ - -void DumpDCD::fwrite_int32(FILE* fp, uint32_t i) { - fwrite(&i,4,1,fp); + if (flush_flag) fflush(fp); } diff --git a/src/dump_dcd.h b/src/dump_dcd.h index db0eb59be9..1e1411467f 100644 --- a/src/dump_dcd.h +++ b/src/dump_dcd.h @@ -30,6 +30,7 @@ class DumpDCD : public Dump { private: int natoms,ntotal,headerflag,nevery_save,nframes; float *coords,*xf,*yf,*zf; + int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no void openfile(); void write_header(int); @@ -37,9 +38,9 @@ class DumpDCD : public Dump { int pack(); void write_data(int, double *); + int modify_param(int, char **); void write_frame(); void write_dcd_header(const char *); - void fwrite_int32(FILE *, uint32_t); }; } diff --git a/src/style_user_ackland.h b/src/style_user_ackland.h index e69de29bb2..6e7483a9f7 100644 --- a/src/style_user_ackland.h +++ b/src/style_user_ackland.h @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + 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 ComputeInclude +#include "compute_ackland_atom.h" +#endif + +#ifdef ComputeClass +ComputeStyle(ackland/atom,ComputeAcklandAtom) +#endif diff --git a/src/style_user_ewaldn.h b/src/style_user_ewaldn.h index e69de29bb2..3eafa50744 100644 --- a/src/style_user_ewaldn.h +++ b/src/style_user_ewaldn.h @@ -0,0 +1,30 @@ +/* ---------------------------------------------------------------------- + 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 KSpaceInclude +#include "ewald_n.h" +#endif + +#ifdef KSpaceClass +KSpaceStyle(ewald/n,EwaldN) +#endif + +#ifdef PairInclude +#include "pair_buck_coul.h" +#include "pair_lj_coul.h" +#endif + +#ifdef PairClass +PairStyle(buck/coul,PairBuckCoul) +PairStyle(lj/coul,PairLJCoul) +#endif