commit 43d2d74468b8ee17668a2cf1500761ade16e2905 Author: daniel Date: Thu Apr 16 13:03:33 2020 +0200 c_abi diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..96ef6c0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target +Cargo.lock diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..f1c63cd --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "xdrfile" +version = "0.1.0" +authors = ["Daniel Bauer "] +license = "GNU" +edition = "2018" +description = "Wrapper around the gromacs libxdrfile library. Can be used to read and write gromacs trajectories in xtc and trr format." +build = "build.rs" + +[dependencies] +failure = "0.1" +lazy-init = "0.3" + +[dev-dependencies] +tempfile = "3.1.0" + +[build-dependencies] +cc = { version = "1.0", features = ["parallel" ]} diff --git a/build.rs b/build.rs new file mode 100644 index 0000000..59de77d --- /dev/null +++ b/build.rs @@ -0,0 +1,17 @@ +extern crate cc; + +use std::fs; + +fn main() { + // This builds gromacs' xdrfile library + let source_files: Vec<_> = fs::read_dir("external/xdrfile/src") + .unwrap() + .map(|f| f.unwrap()) + .map(|f| f.path()) + .collect(); + cc::Build::new() + .files(source_files) + .include("external/xdrfile/include") + .warnings(false) + .compile("libxdrfile.a") +} diff --git a/external/xdrfile/README b/external/xdrfile/README new file mode 100644 index 0000000..d6d3f42 --- /dev/null +++ b/external/xdrfile/README @@ -0,0 +1,21 @@ +Low-level C libraries for manipulating GROMACS XTC and TRR files. +This code is is taken from mdtraj library which is available at + +https://github.com/mdtraj/mdtraj + +and originally derived from xdrfile-1.1.4, available at + +http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library +ftp://ftp.gromacs.org/pub/contrib/xdrfile-1.1.4.tar.gz + +These files are licensed under a BSD 2-clause license, and copyright +Erik Lindahl, David van der Spoel, Robert T. McGibbon. + +Changes from upsteam's xdrfile-1.1.4 include: + - More descriptive error strings printed inside xdrfile.c + - Fix for a segfault on malformed xtc files inside xdrfile_decompress_coord_float (https://github.com/mdtraj/mdtraj/pull/607, https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2014-September/007942.html) + - Addition of read_xtc_nframes function in xdrfile_xtc.c + - Addition of read_trr_nframes function in xdrfile_trr.c + - Bugfix in do_trnheader to return the appropriate error code when reading magic, and properly check the value of the magic in xdrfile_trr.c + - Bugfix of float exception (divide by zero) in xdrfile.c, see https://github.com/SimTk/mdtraj/issues/616 + - Implemented efficient seeking pattern inspired by xdrlib2 (part of MDAnalysis) \ No newline at end of file diff --git a/external/xdrfile/include/ms_stdint.h b/external/xdrfile/include/ms_stdint.h new file mode 100644 index 0000000..cb2acd9 --- /dev/null +++ b/external/xdrfile/include/ms_stdint.h @@ -0,0 +1,259 @@ +// ISO C9x compliant stdint.h for Microsoft Visual Studio +// Based on ISO/IEC 9899:TC2 Committee draft (May 6, 2005) WG14/N1124 +// +// Copyright (c) 2006-2013 Alexander Chemeris +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the product nor the names of its contributors may +// be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED +// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO +// EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef _MSC_VER // [ +#error "Use this header only with Microsoft Visual C++ compilers!" +#endif // _MSC_VER ] + +#ifndef _MSC_STDINT_H_ // [ +#define _MSC_STDINT_H_ + +#if _MSC_VER > 1000 +#pragma once +#endif + +#if _MSC_VER >= 1600 // [ +#include +#else // ] _MSC_VER >= 1600 [ + +#include + +// For Visual Studio 6 in C++ mode and for many Visual Studio versions when +// compiling for ARM we should wrap include with 'extern "C++" {}' +// or compiler give many errors like this: +// error C2733: second C linkage of overloaded function 'wmemchr' not allowed +#ifdef __cplusplus +extern "C" { +#endif +# include +#ifdef __cplusplus +} +#endif + +// Define _W64 macros to mark types changing their size, like intptr_t. +#ifndef _W64 +# if !defined(__midl) && (defined(_X86_) || defined(_M_IX86)) && _MSC_VER >= 1300 +# define _W64 __w64 +# else +# define _W64 +# endif +#endif + + +// 7.18.1 Integer types + +// 7.18.1.1 Exact-width integer types + +// Visual Studio 6 and Embedded Visual C++ 4 doesn't +// realize that, e.g. char has the same size as __int8 +// so we give up on __intX for them. +#if (_MSC_VER < 1300) + typedef signed char int8_t; + typedef signed short int16_t; + typedef signed int int32_t; + typedef unsigned char uint8_t; + typedef unsigned short uint16_t; + typedef unsigned int uint32_t; +#else + typedef signed __int8 int8_t; + typedef signed __int16 int16_t; + typedef signed __int32 int32_t; + typedef unsigned __int8 uint8_t; + typedef unsigned __int16 uint16_t; + typedef unsigned __int32 uint32_t; +#endif +typedef signed __int64 int64_t; +typedef unsigned __int64 uint64_t; + + +// 7.18.1.2 Minimum-width integer types +typedef int8_t int_least8_t; +typedef int16_t int_least16_t; +typedef int32_t int_least32_t; +typedef int64_t int_least64_t; +typedef uint8_t uint_least8_t; +typedef uint16_t uint_least16_t; +typedef uint32_t uint_least32_t; +typedef uint64_t uint_least64_t; + +// 7.18.1.3 Fastest minimum-width integer types +typedef int8_t int_fast8_t; +typedef int16_t int_fast16_t; +typedef int32_t int_fast32_t; +typedef int64_t int_fast64_t; +typedef uint8_t uint_fast8_t; +typedef uint16_t uint_fast16_t; +typedef uint32_t uint_fast32_t; +typedef uint64_t uint_fast64_t; + +// 7.18.1.4 Integer types capable of holding object pointers +#ifdef _WIN64 // [ + typedef signed __int64 intptr_t; + typedef unsigned __int64 uintptr_t; +#else // _WIN64 ][ + typedef _W64 signed int intptr_t; + typedef _W64 unsigned int uintptr_t; +#endif // _WIN64 ] + +// 7.18.1.5 Greatest-width integer types +typedef int64_t intmax_t; +typedef uint64_t uintmax_t; + + +// 7.18.2 Limits of specified-width integer types + +#if !defined(__cplusplus) || defined(__STDC_LIMIT_MACROS) // [ See footnote 220 at page 257 and footnote 221 at page 259 + +// 7.18.2.1 Limits of exact-width integer types +#define INT8_MIN ((int8_t)_I8_MIN) +#define INT8_MAX _I8_MAX +#define INT16_MIN ((int16_t)_I16_MIN) +#define INT16_MAX _I16_MAX +#define INT32_MIN ((int32_t)_I32_MIN) +#define INT32_MAX _I32_MAX +#define INT64_MIN ((int64_t)_I64_MIN) +#define INT64_MAX _I64_MAX +#define UINT8_MAX _UI8_MAX +#define UINT16_MAX _UI16_MAX +#define UINT32_MAX _UI32_MAX +#define UINT64_MAX _UI64_MAX + +// 7.18.2.2 Limits of minimum-width integer types +#define INT_LEAST8_MIN INT8_MIN +#define INT_LEAST8_MAX INT8_MAX +#define INT_LEAST16_MIN INT16_MIN +#define INT_LEAST16_MAX INT16_MAX +#define INT_LEAST32_MIN INT32_MIN +#define INT_LEAST32_MAX INT32_MAX +#define INT_LEAST64_MIN INT64_MIN +#define INT_LEAST64_MAX INT64_MAX +#define UINT_LEAST8_MAX UINT8_MAX +#define UINT_LEAST16_MAX UINT16_MAX +#define UINT_LEAST32_MAX UINT32_MAX +#define UINT_LEAST64_MAX UINT64_MAX + +// 7.18.2.3 Limits of fastest minimum-width integer types +#define INT_FAST8_MIN INT8_MIN +#define INT_FAST8_MAX INT8_MAX +#define INT_FAST16_MIN INT16_MIN +#define INT_FAST16_MAX INT16_MAX +#define INT_FAST32_MIN INT32_MIN +#define INT_FAST32_MAX INT32_MAX +#define INT_FAST64_MIN INT64_MIN +#define INT_FAST64_MAX INT64_MAX +#define UINT_FAST8_MAX UINT8_MAX +#define UINT_FAST16_MAX UINT16_MAX +#define UINT_FAST32_MAX UINT32_MAX +#define UINT_FAST64_MAX UINT64_MAX + +// 7.18.2.4 Limits of integer types capable of holding object pointers +#ifdef _WIN64 // [ +# define INTPTR_MIN INT64_MIN +# define INTPTR_MAX INT64_MAX +# define UINTPTR_MAX UINT64_MAX +#else // _WIN64 ][ +# define INTPTR_MIN INT32_MIN +# define INTPTR_MAX INT32_MAX +# define UINTPTR_MAX UINT32_MAX +#endif // _WIN64 ] + +// 7.18.2.5 Limits of greatest-width integer types +#define INTMAX_MIN INT64_MIN +#define INTMAX_MAX INT64_MAX +#define UINTMAX_MAX UINT64_MAX + +// 7.18.3 Limits of other integer types + +#ifdef _WIN64 // [ +# define PTRDIFF_MIN _I64_MIN +# define PTRDIFF_MAX _I64_MAX +#else // _WIN64 ][ +# define PTRDIFF_MIN _I32_MIN +# define PTRDIFF_MAX _I32_MAX +#endif // _WIN64 ] + +#define SIG_ATOMIC_MIN INT_MIN +#define SIG_ATOMIC_MAX INT_MAX + +#ifndef SIZE_MAX // [ +# ifdef _WIN64 // [ +# define SIZE_MAX _UI64_MAX +# else // _WIN64 ][ +# define SIZE_MAX _UI32_MAX +# endif // _WIN64 ] +#endif // SIZE_MAX ] + +// WCHAR_MIN and WCHAR_MAX are also defined in +#ifndef WCHAR_MIN // [ +# define WCHAR_MIN 0 +#endif // WCHAR_MIN ] +#ifndef WCHAR_MAX // [ +# define WCHAR_MAX _UI16_MAX +#endif // WCHAR_MAX ] + +#define WINT_MIN 0 +#define WINT_MAX _UI16_MAX + +#endif // __STDC_LIMIT_MACROS ] + + +// 7.18.4 Limits of other integer types + +#if !defined(__cplusplus) || defined(__STDC_CONSTANT_MACROS) // [ See footnote 224 at page 260 + +// 7.18.4.1 Macros for minimum-width integer constants + +#define INT8_C(val) val##i8 +#define INT16_C(val) val##i16 +#define INT32_C(val) val##i32 +#define INT64_C(val) val##i64 + +#define UINT8_C(val) val##ui8 +#define UINT16_C(val) val##ui16 +#define UINT32_C(val) val##ui32 +#define UINT64_C(val) val##ui64 + +// 7.18.4.2 Macros for greatest-width integer constants +// These #ifndef's are needed to prevent collisions with . +// Check out Issue 9 for the details. +#ifndef INTMAX_C // [ +# define INTMAX_C INT64_C +#endif // INTMAX_C ] +#ifndef UINTMAX_C // [ +# define UINTMAX_C UINT64_C +#endif // UINTMAX_C ] + +#endif // __STDC_CONSTANT_MACROS ] + +#endif // _MSC_VER >= 1600 ] + +#endif // _MSC_STDINT_H_ ] diff --git a/external/xdrfile/include/trr_header.h b/external/xdrfile/include/trr_header.h new file mode 100644 index 0000000..dedfc77 --- /dev/null +++ b/external/xdrfile/include/trr_header.h @@ -0,0 +1,30 @@ +#ifndef _trr_header_h_ +#define _trr_header_h_ + +typedef struct /* This struct describes the order and the */ +/* sizes of the structs in a trjfile, sizes are given in bytes. */ +{ + int bDouble; /* Double precision? */ + int ir_size; /* Backward compatibility */ + int e_size; /* Backward compatibility */ + int box_size; /* Non zero if a box is present */ + int vir_size; /* Backward compatibility */ + int pres_size; /* Backward compatibility */ + int top_size; /* Backward compatibility */ + int sym_size; /* Backward compatibility */ + int x_size; /* Non zero if coordinates are present */ + int v_size; /* Non zero if velocities are present */ + int f_size; /* Non zero if forces are present */ + + int natoms; /* The total number of atoms */ + int step; /* Current step number */ + int nre; /* Backward compatibility */ + float tf; /* Current time */ + float lambdaf; /* Current value of lambda */ + double td; /* Current time */ + double lambdad; /* Current value of lambda */ +} t_trnheader; + +int do_trnheader(XDRFILE *xd, char bRead, t_trnheader *sh); + +#endif diff --git a/external/xdrfile/include/xdr_seek.h b/external/xdrfile/include/xdr_seek.h new file mode 100644 index 0000000..9cadbbb --- /dev/null +++ b/external/xdrfile/include/xdr_seek.h @@ -0,0 +1,17 @@ +#ifndef _xdr_seek_h +#define _xdr_seek_h + +// for int64_t on older M$ Visual Studio +#if _MSC_VER && _MSVC_VER < 1600 && !__INTEL_COMPILER + #include "ms_stdint.h" +#else + #include +#endif + +#include "xdrfile.h" + +int64_t xdr_tell(XDRFILE *xd); +int xdr_seek(XDRFILE *xd, int64_t pos, int whence); +int xdr_flush(XDRFILE* xd); + +#endif diff --git a/external/xdrfile/include/xdrfile.h b/external/xdrfile/include/xdrfile.h new file mode 100644 index 0000000..5272f7d --- /dev/null +++ b/external/xdrfile/include/xdrfile.h @@ -0,0 +1,632 @@ + /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/*! \file xdrfile.h + * \brief Interface to read/write portabile binary files using XDR. + * + * This file provides an interface to read & write portably binary files, + * using XDR - the external data representation standard defined in RFC 1014. + * + * There are several advantages to the XDR approach: + * + * -# It is portable. And not just portable between big/small integer endian, + * but truly portable if you have system XDR routines. For example: + * - It doesn't matter if the character representation is ASCII or EBCDIC. + * - Some systems are small endian but use big endian order of the two + * dword in a double precision floating-point variable. The system XDR + * libraries will read/write this correctly. + * - Some systems (VAX...) don't use IEEE floating point. Their system + * XDR libraries will convert to/from this automatically. + * -# XDR libraries are required for NFS and lots of other network functions. + * This means there isn't a single Unix-like system that doesn't have them. + * -# There is NO extra metadata whatsoever, and we write plain XDR files. + * If you write a float, it will take exactly 4 bytes in the file. + * (All basic datatypes are 4 bytes, double fp 8 bytes). + * -# You can read/write the files by calling the system XDR routines directly + * too - you don't have to use the routines defined in this file. + * -# It is no problem if your system doesn't have XDR libraries (MS Windows). + * We have written our own versions of the necessary routines that work if + * your system uses ASCII for strings and IEEE floating-point. All types + * of byte and dword endian for integer and floating-point are supported. + * -# You can use these routines for any type of data, but since we designed + * them for Gromacs we also provide a special routine to write coordinates + * with (adjustable) lossy compression. The default precision will give you + * three decimals guaranteed accuracy, and reduces the filesize to 1/10th + * of normal binary data. + * + * We do not support getting or setting positions in XDR files, since it can + * break in horrible ways for large (64-bit) files, resulting in silent data + * corruption. Note that it works great to open/read/write 64-bit files if + * your system supports it; it is just the random access we cannot trust! + * + * We also provide wrapper routines so this module can be used from FORTRAN - + * see the file xdrfile_fortran.txt in the Gromacs distribution for + * documentation on the FORTRAN interface! + */ + +#ifndef _XDRFILE_H_ +#define _XDRFILE_H_ + + +#ifdef __cplusplus +extern "C" +{ +#endif + + /*! \brief Abstract datatype for an portable binary file handle + * + * This datatype essentially works just like the standard FILE type in C. + * The actual contents is hidden in the implementation, so you can only + * define pointers to it, for use with the xdrfile routines. + * + * If you \a really need to see the definition it is in xdrfile.c, but you + * cannot access elements of the structure outside that file. + * + * \warning The implementation is completely different from the C standard + * library FILE, so don't even think about using an XDRFILE pointer as an + * argument to a routine that needs a standard FILE pointer. + */ + typedef struct XDRFILE XDRFILE; + + enum { exdrOK, exdrHEADER, exdrSTRING, exdrDOUBLE, + exdrINT, exdrFLOAT, exdrUINT, exdr3DX, exdrCLOSE, exdrMAGIC, + exdrNOMEM, exdrENDOFFILE, exdrFILENOTFOUND, exdrNR }; + + extern char *exdr_message[exdrNR]; + +#define DIM 3 + typedef float matrix[DIM][DIM]; + typedef float rvec[DIM]; + typedef int mybool; + + + /*! \brief Open a portable binary file, just like fopen() + * + * Use this routine much like calls to the standard library function + * fopen(). The only difference is that the returned pointer should only + * be used with routines defined in this header. + * + * \param path Full or relative path (including name) of the file + * \param mode "r" for reading, "w" for writing, "a" for append. + * + * \return Pointer to abstract xdr file datatype, or NULL if an error occurs. + * + */ + XDRFILE * + xdrfile_open (const char * path, + const char * mode); + + + /*! \brief Close a previously opened portable binary file, just like fclose() + * + * Use this routine much like calls to the standard library function + * fopen(). The only difference is that it is used for an XDRFILE handle + * instead of a FILE handle. + * + * \param xfp Pointer to an abstract XDRFILE datatype + * + * \return 0 on success, non-zero on error. + */ + int + xdrfile_close (XDRFILE * xfp); + + + + + /*! \brief Read one or more \a char type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of characters to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of characters read + */ + int + xdrfile_read_char(char * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a characters type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of characters to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of characters written + */ + int + xdrfile_write_char(char * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read one or more \a unsigned \a char type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of unsigned characters to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned characters read + */ + int + xdrfile_read_uchar(unsigned char * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a unsigned \a characters type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of unsigned characters to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned characters written + */ + int + xdrfile_write_uchar(unsigned char * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read one or more \a short type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of shorts to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of shorts read + */ + int + xdrfile_read_short(short * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a short type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of shorts to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of shorts written + */ + int + xdrfile_write_short(short * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read one or more \a unsigned \a short type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of unsigned shorts to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned shorts read + */ + int + xdrfile_read_ushort(unsigned short * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a unsigned \a short type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of unsigned shorts to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned shorts written + */ + int + xdrfile_write_ushort(unsigned short * ptr, + int ndata, + XDRFILE * xfp); + + + /*! \brief Read one or more \a integer type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of integers to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of integers read + * + * The integer data type is assumed to be less than or equal to 32 bits. + * + * We do not provide any routines for reading/writing 64-bit integers, since + * - Not all XDR implementations support it + * - Not all machines have 64-bit integers + * + * Split your 64-bit data into two 32-bit integers for portability! + */ + int + xdrfile_read_int(int * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a integer type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of integers to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of integers written + * + * The integer data type is assumed to be less than or equal to 32 bits. + * + * We do not provide any routines for reading/writing 64-bit integers, since + * - Not all XDR implementations support it + * - Not all machines have 64-bit integers + * + * Split your 64-bit data into two 32-bit integers for portability! + */ + int + xdrfile_write_int(int * ptr, + int ndata, + XDRFILE * xfp); + + /*! \brief Read one or more \a unsigned \a integers type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of unsigned integers to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned integers read + * + * The integer data type is assumed to be less than or equal to 32 bits. + * + * We do not provide any routines for reading/writing 64-bit integers, since + * - Not all XDR implementations support it + * - Not all machines have 64-bit integers + * + * Split your 64-bit data into two 32-bit integers for portability! + */ + int + xdrfile_read_uint(unsigned int * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a unsigned \a integer type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of unsigned integers to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of unsigned integers written + * + * The integer data type is assumed to be less than or equal to 32 bits. + * + * We do not provide any routines for reading/writing 64-bit integers, since + * - Not all XDR implementations support it + * - Not all machines have 64-bit integers + * + * Split your 64-bit data into two 32-bit integers for portability! + */ + int + xdrfile_write_uint(unsigned int * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read one or more \a float type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of floats to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of floats read + */ + int + xdrfile_read_float(float * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a float type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of floats to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of floats written + */ + int + xdrfile_write_float(float * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read one or more \a double type variable(s) + * + * \param ptr Pointer to memory where data should be written + * \param ndata Number of doubles to read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of doubles read + */ + int + xdrfile_read_double(double * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Write one or more \a double type variable(s) + * + * \param ptr Pointer to memory where data should be read + * \param ndata Number of double to write. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of doubles written + */ + int + xdrfile_write_double(double * ptr, + int ndata, + XDRFILE * xfp); + + + + /*! \brief Read a string (array of characters) + * + * \param ptr Pointer to memory where data should be written + * \param maxlen Maximum length of string. If no end-of-string is encountered, + * one byte less than this is read and end-of-string appended. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of characters read, including end-of-string + */ + int + xdrfile_read_string(char * ptr, + int maxlen, + XDRFILE * xfp); + + + + /*! \brief Write a string (array of characters) + * + * \param ptr Pointer to memory where data should be read + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of characters written, including end-of-string + */ + int + xdrfile_write_string(char * ptr, + XDRFILE * xfp); + + + + /*! \brief Read raw bytes from file (unknown datatype) + * + * \param ptr Pointer to memory where data should be written + * \param nbytes Number of bytes to read. No conversion whatsoever is done. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of bytes read from file + */ + int + xdrfile_read_opaque(char * ptr, + int nbytes, + XDRFILE * xfp); + + + + + /*! \brief Write raw bytes to file (unknown datatype) + * + * \param ptr Pointer to memory where data should be read + * \param nbytes Number of bytes to write. No conversion whatsoever is done. + * \param xfp Handle to portable binary file, created with xdrfile_open() + * + * \return Number of bytes written to file + */ + int + xdrfile_write_opaque(char * ptr, + int nbytes, + XDRFILE * xfp); + + + + + + + /*! \brief Compress coordiates in a float array to XDR file + * + * This routine will perform \a lossy compression on the three-dimensional + * coordinate data data specified and store it in the XDR file. + * + * The lossy part of the compression consists of multiplying each + * coordinate with the precision argument and then rounding to integers. + * We suggest a default value of 1000.0, which means you are guaranteed + * three decimals of accuracy. The only limitation is that scaled coordinates + * must still fit in an integer variable, so if the precision is 1000.0 the + * coordinate magnitudes must be less than +-2e6. + * + * \param ptr Pointer to coordinates to compress (length 3*ncoord) + * \param ncoord Number of coordinate triplets in data + * \param precision Scaling factor for lossy compression. If it is <=0, + * the default value of 1000.0 is used. + * \param xfp Handle to portably binary file + * + * \return Number of coordinate triplets written. + * IMPORTANT: Check that this is equal to ncoord - if it is + * negative, an error occured. This should not happen with + * normal data, but if your coordinates are NaN or very + * large (>1e6) it is not possible to use the compression. + * + * \warning The compression algorithm is not part of the XDR standard, + * and very complicated, so you will need this xdrfile module + * to read it later. + */ + int + xdrfile_compress_coord_float(float * ptr, + int ncoord, + float precision, + XDRFILE * xfp); + + + + + /*! \brief Decompress coordiates from XDR file to array of floats + * + * This routine will decompress three-dimensional coordinate data previously + * stored in an XDR file and store it in the specified array of floats. + * + * The precision used during the earlier compression is read from the file + * and returned - you cannot adjust the accuracy at this stage. + * + * \param ptr Pointer to coordinates to compress (length>= 3*ncoord) + * \param ncoord Max number of coordinate triplets to read on input, actual + * number of coordinate triplets read on return. If this + * is smaller than the number of coordinates in the frame an + * error will occur. + * \param precision The precision used in the previous compression will be + * written to this variable on return. + * \param xfp Handle to portably binary file + * + * \return Number of coordinate triplets read. If this is negative, + * an error occured. + * + * \warning Since we cannot count on being able to set/get the + * position of large files (>2Gb), it is not possible to + * recover from errors by re-reading the frame if the + * storage area you provided was too small. To avoid this + * from happening, we recommend that you store the number of + * coordinates triplet as an integer either in a header or + * just before the compressed coordinate data, so you can + * read it first and allocated enough memory. + */ + int + xdrfile_decompress_coord_float(float * ptr, + int * ncoord, + float * precision, + XDRFILE * xfp); + + + + + /*! \brief Compress coordiates in a double array to XDR file + * + * This routine will perform \a lossy compression on the three-dimensional + * coordinate data data specified and store it in the XDR file. Double will + * NOT give you any extra precision since the coordinates are compressed. This + * routine just avoids allocating a temporary array of floats. + * + * The lossy part of the compression consists of multiplying each + * coordinate with the precision argument and then rounding to integers. + * We suggest a default value of 1000.0, which means you are guaranteed + * three decimals of accuracy. The only limitation is that scaled coordinates + * must still fit in an integer variable, so if the precision is 1000.0 the + * coordinate magnitudes must be less than +-2e6. + * + * \param ptr Pointer to coordinates to compress (length 3*ncoord) + * \param ncoord Number of coordinate triplets in data + * \param precision Scaling factor for lossy compression. If it is <=0, the + * default value of 1000.0 is used. + * \param xfp Handle to portably binary file + * + * \return Number of coordinate triplets written. + * IMPORTANT: Check that this is equal to ncoord - if it is + * negative, an error occured. This should not happen with + * normal data, but if your coordinates are NaN or very + * large (>1e6) it is not possible to use the compression. + * + * \warning The compression algorithm is not part of the XDR standard, + * and very complicated, so you will need this xdrfile module + * to read it later. + */ + int + xdrfile_compress_coord_double(double * ptr, + int ncoord, + double precision, + XDRFILE * xfp); + + + + + /*! \brief Decompress coordiates from XDR file to array of doubles + * + * This routine will decompress three-dimensional coordinate data previously + * stored in an XDR file and store it in the specified array of doubles. + * Double will NOT give you any extra precision since the coordinates are + * compressed. This routine just avoids allocating a temporary array of floats. + * + * The precision used during the earlier compression is read from the file + * and returned - you cannot adjust the accuracy at this stage. + * + * \param ptr Pointer to coordinates to compress (length>= 3*ncoord) + * \param ncoord Max number of coordinate triplets to read on input, actual + * number of coordinate triplets read on return. If this + * is smaller than the number of coordinates in the frame an + * error will occur. + * \param precision The precision used in the previous compression will be + * written to this variable on return. + * \param xfp Handle to portably binary file + * + * \return Number of coordinate triplets read. If this is negative, + * an error occured. + * + * \warning Since we cannot count on being able to set/get the + * position of large files (>2Gb), it is not possible to + * recover from errors by re-reading the frame if the + * storage area you provided was too small. To avoid this + * from happening, we recommend that you store the number of + * coordinates triplet as an integer either in a header or + * just before the compressed coordinate data, so you can + * read it first and allocated enough memory. + */ + int + xdrfile_decompress_coord_double(double * ptr, + int * ncoord, + double * precision, + XDRFILE * xfp); + + + +#ifdef __cplusplus +} +#endif + +#endif /* _XDRFILE_H_ */ + diff --git a/external/xdrfile/include/xdrfile_trr.h b/external/xdrfile/include/xdrfile_trr.h new file mode 100644 index 0000000..14fe105 --- /dev/null +++ b/external/xdrfile/include/xdrfile_trr.h @@ -0,0 +1,61 @@ + /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef _xdrfile_trr_h +#define _xdrfile_trr_h + +#ifdef __cplusplus +extern "C" { +#endif + +#include "xdrfile.h" + + /* All functions return exdrOK if successful. + * (error codes defined in xdrfile.h). + */ + + /* This function returns the number of atoms in the xtc file in *natoms */ + extern int read_trr_natoms(char *fn,int *natoms); + extern int read_trr_nframes(char* fn, unsigned long *nframes); + + /* Read one frame of an open xtc file. If either of x,v,f,box are + NULL the arrays will be read from the file but not used. */ + extern int read_trr(XDRFILE *xd,int natoms,int *step,float *t,float *lambda, + matrix box,rvec *x,rvec *v,rvec *f); + + /* Write a frame to xtc file */ + extern int write_trr(XDRFILE *xd,int natoms,int step,float t,float lambda, + matrix box,rvec *x,rvec *v,rvec *f); + + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/external/xdrfile/include/xdrfile_xtc.h b/external/xdrfile/include/xdrfile_xtc.h new file mode 100644 index 0000000..ff6f4bd --- /dev/null +++ b/external/xdrfile/include/xdrfile_xtc.h @@ -0,0 +1,61 @@ + /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef _xdrfile_xtc_h +#define _xdrfile_xtc_h + +#ifdef __cplusplus +extern "C" { +#endif + +#include "xdrfile.h" + + /* All functions return exdrOK if succesfull. + * (error codes defined in xdrfile.h). + */ + + /* This function returns the number of atoms in the xtc file in *natoms */ + extern int read_xtc_natoms(char *fn,int *natoms); + + int read_xtc_nframes(char* fn, unsigned long *nframes); + + /* Read one frame of an open xtc file */ + extern int read_xtc(XDRFILE *xd,int natoms,int *step,float *time, + matrix box,rvec *x,float *prec); + + /* Write a frame to xtc file */ + extern int write_xtc(XDRFILE *xd, + int natoms,int step,float time, + matrix box,rvec *x,float prec); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/external/xdrfile/src/xdr_seek.c b/external/xdrfile/src/xdr_seek.c new file mode 100644 index 0000000..17602ae --- /dev/null +++ b/external/xdrfile/src/xdr_seek.c @@ -0,0 +1,55 @@ +/* 64 bit fileseek operations */ +#define _FILE_OFFSETS_BITS 64 +#include "xdr_seek.h" +#include + +/// copied from xtcfile.c (version 1.1.4) +struct XDRFILE +{ + FILE * fp; /**< pointer to standard C library file handle */ + void * /*this used to be (XDR*) */ xdr; /**< pointer to corresponding XDR handle */ + char mode; /**< r=read, w=write, a=append */ + int * buf1; /**< Buffer for internal use */ + int buf1size; /**< Current allocated length of buf1 */ + int * buf2; /**< Buffer for internal use */ + int buf2size; /**< Current allocated length of buf2 */ +}; +//// end of copied + +int64_t xdr_tell(XDRFILE *xd) +{ + FILE* fptr = xd->fp; + +#ifndef _WIN32 + // use posix 64 bit ftell version + return ftello(fptr); +#elif defined(_MSVC_VER) && !__INTEL_COMPILER + return _ftelli64(fptr); +#else + return ftell(fptr); +#endif +} + +int xdr_seek(XDRFILE *xd, int64_t pos, int whence) +{ + int result = 1; + FILE* fptr = xd->fp; + +#ifndef _WIN32 + // use posix 64 bit ftell version + result = fseeko(fptr, pos, whence) < 0 ? exdrNR : exdrOK; +#elif _MSVC_VER && !__INTEL_COMPILER + result = _fseeki64(fptr, pos, whence) < 0 ? exdrNR : exdrOK; +#else + result = fseek(fptr, pos, whence) < 0 ? exdrNR : exdrOK; +#endif + if (result != exdrOK) + return result; + + return exdrOK; +} + +int xdr_flush(XDRFILE* xdr) +{ + return fflush(xdr->fp); +} diff --git a/external/xdrfile/src/xdrfile.c b/external/xdrfile/src/xdrfile.c new file mode 100644 index 0000000..24285f3 --- /dev/null +++ b/external/xdrfile/src/xdrfile.c @@ -0,0 +1,2634 @@ + /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* Get HAVE_RPC_XDR_H, F77_FUNC from config.h if available */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#define _FILE_OFFSET_BITS 64 + +/* get fixed-width types if we are using ANSI C99 */ +#ifdef HAVE_STDINT_H +# include +#elif (defined HAVE_INTTYPES_H) +# include +#endif + + +#ifdef HAVE_RPC_XDR_H +# include +# include +#endif + +#include "xdrfile.h" + +/* Default FORTRAN name mangling is: lower case name, append underscore */ +#ifndef F77_FUNC +#define F77_FUNC(name,NAME) name ## _ +#endif + +char *exdr_message[exdrNR] = { + "OK", + "Header", + "String", + "Double", + "Integer", + "Float", + "Unsigned integer", + "Compressed 3D coordinate", + "Closing file", + "Magic number", + "Not enough memory", + "End of file", + "File not found" +}; + +/* + * Declare our own XDR routines statically if no libraries are present. + * Actual implementation is at the end of this file. + * + * We don't want the low-level XDR implementation as part of the Gromacs + * documentation, so skip it for doxygen too... + */ +#if (!defined HAVE_RPC_XDR_H && !defined DOXYGEN) + +enum xdr_op +{ + XDR_ENCODE = 0, + XDR_DECODE = 1, + XDR_FREE = 2 +}; + + +/* We need integer types that are guaranteed to be 4 bytes wide. + * If ANSI C99 headers were included they are already defined + * as int32_t and uint32_t. Check, and if not define them ourselves. + * Since it is just our workaround for missing ANSI C99 types, avoid adding + * it to the doxygen documentation. + */ +#if !(defined INT32_MAX || defined DOXYGEN) +# if (INT_MAX == 2147483647) +# define int32_t int +# define uint32_t unsigned int +# define INT32_MAX 2147483647 +# elif (LONG_MAX == 2147483647) +# define int32_t long +# define uint32_t unsigned long +# define INT32_MAX 2147483647L +# else +# error ERROR: No 32 bit wide integer type found! +# error Use system XDR libraries instead, or update xdrfile.c +# endif +#endif + +typedef struct XDR XDR; + +struct XDR +{ + enum xdr_op x_op; + struct xdr_ops + { + int (*x_getlong) (XDR *__xdrs, int32_t *__lp); + int (*x_putlong) (XDR *__xdrs, int32_t *__lp); + int (*x_getbytes) (XDR *__xdrs, char *__addr, unsigned int __len); + int (*x_putbytes) (XDR *__xdrs, char *__addr, unsigned int __len); + /* two next routines are not 64-bit IO safe - don't use! */ + unsigned int (*x_getpostn) (XDR *__xdrs); + int (*x_setpostn) (XDR *__xdrs, unsigned int __pos); + void (*x_destroy) (XDR *__xdrs); + } + *x_ops; + char *x_private; +}; + +static int xdr_char (XDR *xdrs, char *ip); +static int xdr_u_char (XDR *xdrs, unsigned char *ip); +static int xdr_short (XDR *xdrs, short *ip); +static int xdr_u_short (XDR *xdrs, unsigned short *ip); +static int xdr_int (XDR *xdrs, int *ip); +static int xdr_u_int (XDR *xdrs, unsigned int *ip); +static int xdr_float (XDR *xdrs, float *ip); +static int xdr_double (XDR *xdrs, double *ip); +static int xdr_string (XDR *xdrs, char **ip, unsigned int maxsize); +static int xdr_opaque (XDR *xdrs, char *cp, unsigned int cnt); +static void xdrstdio_create (XDR *xdrs, FILE *fp, enum xdr_op xop); + +#define xdr_getpos(xdrs) \ + (*(xdrs)->x_ops->x_getpostn)(xdrs) +#define xdr_setpos(xdrs, pos) \ + (*(xdrs)->x_ops->x_setpostn)(xdrs, pos) +#define xdr_destroy(xdrs) \ + do { \ + if ((xdrs)->x_ops->x_destroy) \ + (*(xdrs)->x_ops->x_destroy)(xdrs); \ + } while (0) +#endif /* end of our own XDR declarations */ + + + + + +/** Contents of the abstract XDRFILE data structure. + * + * @internal + * + * This structure is used to provide an XDR file interface that is + * virtual identical to the standard UNIX fopen/fread/fwrite/fclose. + */ +struct XDRFILE +{ + FILE * fp; /**< pointer to standard C library file handle */ + XDR * xdr; /**< pointer to corresponding XDR handle */ + char mode; /**< r=read, w=write, a=append */ + int * buf1; /**< Buffer for internal use */ + int buf1size; /**< Current allocated length of buf1 */ + int * buf2; /**< Buffer for internal use */ + int buf2size; /**< Current allocated length of buf2 */ +}; + + + + +/************************************************************* + * Implementation of higher-level routines to read/write * + * portable data based on the XDR standard. These should be * + * called from C - see further down for Fortran77 wrappers. * + *************************************************************/ + +XDRFILE * +xdrfile_open(const char *path, const char *mode) +{ + char newmode[5]; + enum xdr_op xdrmode; + XDRFILE *xfp; + + /* make sure XDR files are opened in binary mode... */ + if(*mode=='w' || *mode=='W') + { + sprintf(newmode,"wb+"); + xdrmode=XDR_ENCODE; + } else if(*mode == 'a' || *mode == 'A') + { + sprintf(newmode,"ab+"); + xdrmode = XDR_ENCODE; + } else if(*mode == 'r' || *mode == 'R') + { + sprintf(newmode,"rb"); + xdrmode = XDR_DECODE; + } else /* cannot determine mode */ + return NULL; + + if((xfp=(XDRFILE *)malloc(sizeof(XDRFILE)))==NULL) + return NULL; + if((xfp->fp=fopen(path,newmode))==NULL) + { + free(xfp); + return NULL; + } + if((xfp->xdr=(XDR *)malloc(sizeof(XDR)))==NULL) + { + fclose(xfp->fp); + free(xfp); + return NULL; + } + xfp->mode=*mode; + xdrstdio_create((XDR *)(xfp->xdr),xfp->fp,xdrmode); + xfp->buf1 = xfp->buf2 = NULL; + xfp->buf1size = xfp->buf2size = 0; + return xfp; +} + +int +xdrfile_close(XDRFILE *xfp) +{ + int ret=exdrCLOSE; + if(xfp) + { + /* flush and destroy XDR stream */ + if(xfp->xdr) + xdr_destroy((XDR *)(xfp->xdr)); + free(xfp->xdr); + /* close the file */ + ret=fclose(xfp->fp); + if(xfp->buf1size) + free(xfp->buf1); + if(xfp->buf2size) + free(xfp->buf2); + free(xfp); + } + return ret; /* return 0 if ok */ +} + + + +int +xdrfile_read_int(int *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_int(int *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + + +int +xdrfile_read_uint(unsigned int *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_uint(unsigned int *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_read_char(char *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_char(char *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + + +int +xdrfile_read_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_read_short(short *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_short(short *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + + +int +xdrfile_read_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + + return i; +} + +int +xdrfile_write_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_read_float(float *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_write_float(float *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_read_double(double *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_write_double(double *ptr, int ndata, XDRFILE* xfp) +{ + int i=0; + /* read write is encoded in the XDR struct */ + while(ixdr),ptr+i)) + i++; + return i; +} + +int +xdrfile_read_string(char *ptr, int maxlen, XDRFILE* xfp) +{ + int i; + if(xdr_string((XDR *)(xfp->xdr),&ptr,maxlen)) { + i=0; + while(ixdr),&ptr,len)) + return len; + else + return 0; +} + + +int +xdrfile_read_opaque(char *ptr, int cnt, XDRFILE* xfp) +{ + if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt)) + return cnt; + else + return 0; +} + + +int +xdrfile_write_opaque(char *ptr, int cnt, XDRFILE* xfp) +{ + if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt)) + return cnt; + else + return 0; +} + + +/* Internal support routines for reading/writing compressed coordinates + * sizeofint - calculate smallest number of bits necessary + * to represent a certain integer. + */ +static int +sizeofint(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 a number of small unsigned integers and the maximum value + * return the number of bits needed to read or write them with the + * routines encodeints/decodeints. You need this parameter when + * calling those routines. + * (However, in some cases we can just use the variable 'smallidx' + * which is the exact number of bits, and them we dont need to call + * this routine). + */ +static int +sizeofints(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; + +} + + +/* + * encodebits - 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 had + * better make sure that this number of bits is enough to hold the value. + * Num must also be positive. + */ +static void +encodebits(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); + } +} + +/* + * encodeints - encode a small set of small integers in compressed format + * + * this routine is used internally by xdr3dfcoord, to encode a set of + * small integers to the buffer for writing to a file. + * 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 whould cause bad results. + * THese things are checked in the calling routines, so make sure not + * to remove those checks... + */ + +static void +encodeints(int buf[], int num_of_ints, 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,"(xdrfile error) major breakdown in encodeints - num %u doesn't " + "match size %u\n", nums[i], sizes[i]); + abort(); + } + /* 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++) + { + encodebits(buf, 8, bytes[i]); + } + encodebits(buf, num_of_bits - num_of_bytes * 8, 0); + } + else + { + for (i = 0; i < num_of_bytes-1; i++) + { + encodebits(buf, 8, bytes[i]); + } + encodebits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]); + } +} + + +/* + * decodebits - 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 +decodebits(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; +} + +/* + * decodeints - decode 'small' integers from the buf array + * + * this routine is the inverse from encodeints() 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 +decodeints(int buf[], 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++] = decodebits(buf, 8); + num_of_bits -= 8; + } + if (num_of_bits > 0) + { + bytes[num_of_bytes++] = decodebits(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); +} + + +static const 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(*magicints)) + +/* Compressed coordinate routines - modified from the original + * implementation by Frans v. Hoesel to make them threadsafe. + */ +int +xdrfile_decompress_coord_float(float *ptr, + int *size, + float *precision, + XDRFILE* xfp) +{ + int minint[3], maxint[3], *lip; + int smallidx, minidx, maxidx; + unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3; + int k, *buf1, *buf2, lsize, flag; + int smallnum, smaller, larger, i, is_smaller, run; + float *lfp, inv_precision; + int tmp, *thiscoord, prevcoord[3]; + unsigned int bitsize; + const float* ptrstart = ptr; + + bitsizeint[0] = 0; + bitsizeint[1] = 0; + bitsizeint[2] = 0; + + if(xfp==NULL || ptr==NULL) { + fprintf(stderr, "(xdrfile error) Null pointer issue\n"); + return -1; + } + tmp=xdrfile_read_int(&lsize,1,xfp); + if(tmp==0) { + fprintf(stderr, "(xdrfile error) Size could not be read\n"); + return -1; /* return if we could not read size */ + } + if (*size < lsize) + { + fprintf(stderr, "(xdrfile error) Requested to decompress %d coords, file contains %d\n", + *size, lsize); + return -1; + } + *size = lsize; + size3 = *size * 3; + if(size3>xfp->buf1size) + { + if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for decompressing coordinates.\n"); + return -1; + } + xfp->buf1size=size3; + xfp->buf2size=size3*1.2; + if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for decompressing coordinates.\n"); + return -1; + } + } + /* Dont bother with compression for three atoms or less */ + if(*size<=9) + { + return xdrfile_read_float(ptr,size3,xfp)/3; + /* return number of coords, not floats */ + } + /* Compression-time if we got here. Read precision first */ + xdrfile_read_float(precision,1,xfp); + + /* avoid repeated pointer dereferencing. */ + buf1=xfp->buf1; + buf2=xfp->buf2; + /* buf2[0-2] are special and do not contain actual data */ + buf2[0] = buf2[1] = buf2[2] = 0; + xdrfile_read_int(minint,3,xfp); + xdrfile_read_int(maxint,3,xfp); + + 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); + } + + if (xdrfile_read_int(&smallidx,1,xfp) == 0) { + fprintf(stderr,"(xdrfile error) Undocumented error 1"); + return 0; /* not sure what has happened here or why we return... */ + } + tmp=smallidx+8; + maxidx = (LASTIDXtmp) ? FIRSTIDX : tmp; + smaller = magicints[tmp] / 2; + smallnum = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + larger = magicints[maxidx]; + + /* buf2[0] holds the length in bytes */ + + if (xdrfile_read_int(buf2,1,xfp) == 0) { + fprintf(stderr, "(xdrfile error) Undocumented error 2"); + return 0; + } + if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0) { + fprintf(stderr, "(xdrfile error) Undocumented error 3"); + return 0; + } + buf2[0] = buf2[1] = buf2[2] = 0; + + lfp = ptr; + inv_precision = 1.0 / * precision; + run = 0; + i = 0; + lip = buf1; + while ( i < lsize ) + { + thiscoord = (int *)(lip) + i * 3; + + if (bitsize == 0) + { + thiscoord[0] = decodebits(buf2, bitsizeint[0]); + thiscoord[1] = decodebits(buf2, bitsizeint[1]); + thiscoord[2] = decodebits(buf2, bitsizeint[2]); + } + else + { + decodeints(buf2, 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 = decodebits(buf2, 1); + is_smaller = 0; + if (flag == 1) + { + run = decodebits(buf2, 5); + is_smaller = run % 3; + run -= is_smaller; + is_smaller--; + } + if ((lfp-ptrstart)+run > size3) + { + fprintf(stderr, "(xdrfile error) Buffer overrun during decompression.\n"); + return 0; + } + if (run > 0) + { + thiscoord += 3; + for (k = 0; k < run; k+=3) + { + decodeints(buf2, 3, smallidx, sizesmall, thiscoord); + i++; + thiscoord[0] += prevcoord[0] - smallnum; + thiscoord[1] += prevcoord[1] - smallnum; + thiscoord[2] += prevcoord[2] - smallnum; + 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) + { + smallnum = smaller; + + if (smallidx > FIRSTIDX) + { + smaller = magicints[smallidx - 1] /2; + } + else + { + smaller = 0; + } + } + else if (is_smaller > 0) + { + smaller = smallnum; + smallnum = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + if (sizesmall[0]==0 || sizesmall[1]==0 || sizesmall[2]==0) + { + fprintf(stderr, "(xdrfile error) Undefined error.\n"); + return 0; + } + } + return *size; +} + +int +xdrfile_compress_coord_float(float *ptr, + int size, + float precision, + XDRFILE* xfp) +{ + 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 k, *buf1, *buf2; + int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun; + float *lfp, lf; + int tmp, tmpsum, *thiscoord, prevcoord[3]; + unsigned int tmpcoord[30]; + int errval=1; + unsigned int bitsize; + + if(xfp==NULL) + return -1; + size3=3*size; + + bitsizeint[0] = 0; + bitsizeint[1] = 0; + bitsizeint[2] = 0; + + if(size3>xfp->buf1size) + { + if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for compressing coordinates.\n"); + return -1; + } + xfp->buf1size=size3; + xfp->buf2size=size3*1.2; + if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for compressing coordinates.\n"); + return -1; + } + } + if(xdrfile_write_int(&size,1,xfp)==0) + return -1; /* return if we could not write size */ + /* Dont bother with compression for three atoms or less */ + if(size<=9) + { + return xdrfile_write_float(ptr,size3,xfp)/3; + /* return number of coords, not floats */ + } + /* Compression-time if we got here. Write precision first */ + if (precision <= 0) + precision = 1000; + xdrfile_write_float(&precision,1,xfp); + /* avoid repeated pointer dereferencing. */ + buf1=xfp->buf1; + buf2=xfp->buf2; + /* buf2[0-2] are special and do not contain actual data */ + buf2[0] = buf2[1] = buf2[2] = 0; + minint[0] = minint[1] = minint[2] = INT_MAX; + maxint[0] = maxint[1] = maxint[2] = INT_MIN; + prevrun = -1; + lfp = ptr; + lip = buf1; + mindiff = INT_MAX; + oldlint1 = oldlint2 = oldlint3 = 0; + while(lfp < ptr + size3 ) + { + /* find nearest integer */ + if (*lfp >= 0.0) + lf = *lfp * precision + 0.5; + else + lf = *lfp * precision - 0.5; + if (fabs(lf) > INT_MAX-2) + { + /* scaling would cause overflow */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + errval=0; + } + lint1 = 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) > INT_MAX-2) + { + /* scaling would cause overflow */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + errval=0; + } + lint2 = 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) > INT_MAX-2) + { + errval=0; + } + lint3 = 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 > ptr + 3) + mindiff = diff; + oldlint1 = lint1; + oldlint2 = lint2; + oldlint3 = lint3; + } + xdrfile_write_int(minint,3,xfp); + xdrfile_write_int(maxint,3,xfp); + + if ((float)maxint[0] - (float)minint[0] >= INT_MAX-2 || + (float)maxint[1] - (float)minint[1] >= INT_MAX-2 || + (float)maxint[2] - (float)minint[2] >= INT_MAX-2) { + /* turning value in unsigned by subtracting minint + * would cause overflow + */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + 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 = buf1; + luip = (unsigned int *) buf1; + smallidx = FIRSTIDX; + while (smallidx < LASTIDX && magicints[smallidx] < mindiff) + { + smallidx++; + } + xdrfile_write_int(&smallidx,1,xfp); + tmp=smallidx+8; + maxidx = (LASTIDXtmp) ? FIRSTIDX : tmp; + smaller = magicints[tmp] / 2; + smallnum = 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]) < smallnum && + abs(thiscoord[1] - thiscoord[4]) < smallnum && + abs(thiscoord[2] - thiscoord[5]) < smallnum) + { + /* 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) + { + encodebits(buf2, bitsizeint[0], tmpcoord[0]); + encodebits(buf2, bitsizeint[1], tmpcoord[1]); + encodebits(buf2, bitsizeint[2], tmpcoord[2]); + } + else + { + encodeints(buf2, 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) + { + tmpsum=0; + for(j=0;j<3;j++) + { + tmp=thiscoord[j] - prevcoord[j]; + tmpsum+=tmp*tmp; + } + if (is_smaller == -1 && tmpsum >= smaller * smaller) + { + is_smaller = 0; + } + + tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum; + tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum; + tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum; + + 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]) < smallnum && + abs(thiscoord[1] - prevcoord[1]) < smallnum && + abs(thiscoord[2] - prevcoord[2]) < smallnum) + { + is_small = 1; + } + } + if (run != prevrun || is_smaller != 0) + { + prevrun = run; + encodebits(buf2, 1, 1); /* flag the change in run-length */ + encodebits(buf2, 5, run+is_smaller+1); + } + else + { + encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */ + } + for (k=0; k < run; k+=3) + { + encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]); + } + if (is_smaller != 0) + { + smallidx += is_smaller; + if (is_smaller < 0) + { + smallnum = smaller; + smaller = magicints[smallidx-1] / 2; + } + else + { + smaller = smallnum; + smallnum = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + } + } + if (buf2[1] != 0) buf2[0]++; + xdrfile_write_int(buf2,1,xfp); /* buf2[0] holds the length in bytes */ + tmp=xdrfile_write_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp); + if(tmp==(unsigned int)buf2[0]) + return size; + else + return -1; +} + + +int +xdrfile_decompress_coord_double(double *ptr, + int *size, + double *precision, + XDRFILE* xfp) +{ + int minint[3], maxint[3], *lip; + int smallidx, minidx, maxidx; + unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3; + int k, *buf1, *buf2, lsize, flag; + int smallnum, smaller, larger, i, is_smaller, run; + double *lfp, inv_precision; + float float_prec, tmpdata[30]; + int tmp, *thiscoord, prevcoord[3]; + unsigned int bitsize; + + bitsizeint[0] = 0; + bitsizeint[1] = 0; + bitsizeint[2] = 0; + + if(xfp==NULL || ptr==NULL) + return -1; + tmp=xdrfile_read_int(&lsize,1,xfp); + if(tmp==0) + return -1; /* return if we could not read size */ + if (*size < lsize) + { + fprintf(stderr, "(xdrfile error) Requested to decompress %d coords, file contains %d\n", + *size, lsize); + return -1; + } + *size = lsize; + size3 = *size * 3; + if(size3>xfp->buf1size) + { + if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for decompression coordinates.\n"); + return -1; + } + xfp->buf1size=size3; + xfp->buf2size=size3*1.2; + if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL) + { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for decompressing coordinates.\n"); + return -1; + } + } + /* Dont bother with compression for three atoms or less */ + if(*size<=9) + { + tmp=xdrfile_read_float(tmpdata,size3,xfp); + for(i=0;i<9*3;i++) + ptr[i]=tmpdata[i]; + return tmp/3; + /* return number of coords, not floats */ + } + /* Compression-time if we got here. Read precision first */ + xdrfile_read_float(&float_prec,1,xfp); + *precision=float_prec; + /* avoid repeated pointer dereferencing. */ + buf1=xfp->buf1; + buf2=xfp->buf2; + /* buf2[0-2] are special and do not contain actual data */ + buf2[0] = buf2[1] = buf2[2] = 0; + xdrfile_read_int(minint,3,xfp); + xdrfile_read_int(maxint,3,xfp); + + 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); + } + + if (xdrfile_read_int(&smallidx,1,xfp) == 0) + return 0; + tmp=smallidx+8; + maxidx = (LASTIDXtmp) ? FIRSTIDX : tmp; + smaller = magicints[tmp] / 2; + smallnum = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + larger = magicints[maxidx]; + + /* buf2[0] holds the length in bytes */ + + if (xdrfile_read_int(buf2,1,xfp) == 0) + return 0; + if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0) + return 0; + buf2[0] = buf2[1] = buf2[2] = 0; + + lfp = ptr; + inv_precision = 1.0 / * precision; + run = 0; + i = 0; + lip = buf1; + while ( i < lsize ) + { + thiscoord = (int *)(lip) + i * 3; + + if (bitsize == 0) + { + thiscoord[0] = decodebits(buf2, bitsizeint[0]); + thiscoord[1] = decodebits(buf2, bitsizeint[1]); + thiscoord[2] = decodebits(buf2, bitsizeint[2]); + } else { + decodeints(buf2, 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 = decodebits(buf2, 1); + is_smaller = 0; + if (flag == 1) + { + run = decodebits(buf2, 5); + is_smaller = run % 3; + run -= is_smaller; + is_smaller--; + } + if (run > 0) + { + thiscoord += 3; + for (k = 0; k < run; k+=3) + { + decodeints(buf2, 3, smallidx, sizesmall, thiscoord); + i++; + thiscoord[0] += prevcoord[0] - smallnum; + thiscoord[1] += prevcoord[1] - smallnum; + thiscoord[2] += prevcoord[2] - smallnum; + 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) { + smallnum = smaller; + if (smallidx > FIRSTIDX) { + smaller = magicints[smallidx - 1] /2; + } else { + smaller = 0; + } + } else if (is_smaller > 0) { + smaller = smallnum; + smallnum = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + } + return *size; +} + +int +xdrfile_compress_coord_double(double *ptr, + int size, + double precision, + XDRFILE* xfp) +{ + 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 k, *buf1, *buf2; + int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun; + double *lfp; + float float_prec, lf,tmpdata[30]; + int tmp, tmpsum, *thiscoord, prevcoord[3]; + unsigned int tmpcoord[30]; + int errval=1; + unsigned int bitsize; + + bitsizeint[0] = 0; + bitsizeint[1] = 0; + bitsizeint[2] = 0; + + if(xfp==NULL) + return -1; + size3=3*size; + if(size3>xfp->buf1size) { + if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL) { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for compressing coordinates.\n"); + return -1; + } + xfp->buf1size=size3; + xfp->buf2size=size3*1.2; + if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL) { + fprintf(stderr, "(xdrfile error) Cannot allocate memory for compressing coordinates.\n"); + return -1; + } + } + if(xdrfile_write_int(&size,1,xfp)==0) + return -1; /* return if we could not write size */ + /* Dont bother with compression for three atoms or less */ + if(size<=9) { + for(i=0;i<9*3;i++) + tmpdata[i]=ptr[i]; + return xdrfile_write_float(tmpdata,size3,xfp)/3; + /* return number of coords, not floats */ + } + /* Compression-time if we got here. Write precision first */ + if (precision <= 0) + precision = 1000; + float_prec=precision; + xdrfile_write_float(&float_prec,1,xfp); + /* avoid repeated pointer dereferencing. */ + buf1=xfp->buf1; + buf2=xfp->buf2; + /* buf2[0-2] are special and do not contain actual data */ + buf2[0] = buf2[1] = buf2[2] = 0; + minint[0] = minint[1] = minint[2] = INT_MAX; + maxint[0] = maxint[1] = maxint[2] = INT_MIN; + prevrun = -1; + lfp = ptr; + lip = buf1; + mindiff = INT_MAX; + oldlint1 = oldlint2 = oldlint3 = 0; + while(lfp < ptr + size3 ) { + /* find nearest integer */ + if (*lfp >= 0.0) + lf = (float)*lfp * float_prec + 0.5; + else + lf = (float)*lfp * float_prec - 0.5; + if (fabs(lf) > INT_MAX-2) { + /* scaling would cause overflow */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + errval=0; + } + lint1 = lf; + if (lint1 < minint[0]) minint[0] = lint1; + if (lint1 > maxint[0]) maxint[0] = lint1; + *lip++ = lint1; + lfp++; + if (*lfp >= 0.0) + lf = (float)*lfp * float_prec + 0.5; + else + lf = (float)*lfp * float_prec - 0.5; + if (fabs(lf) > INT_MAX-2) { + /* scaling would cause overflow */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + errval=0; + } + lint2 = lf; + if (lint2 < minint[1]) minint[1] = lint2; + if (lint2 > maxint[1]) maxint[1] = lint2; + *lip++ = lint2; + lfp++; + if (*lfp >= 0.0) + lf = (float)*lfp * float_prec + 0.5; + else + lf = (float)*lfp * float_prec - 0.5; + if (fabs(lf) > INT_MAX-2) { + errval=0; + } + lint3 = 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 > ptr + 3) + mindiff = diff; + oldlint1 = lint1; + oldlint2 = lint2; + oldlint3 = lint3; + } + xdrfile_write_int(minint,3,xfp); + xdrfile_write_int(maxint,3,xfp); + + if ((float)maxint[0] - (float)minint[0] >= INT_MAX-2 || + (float)maxint[1] - (float)minint[1] >= INT_MAX-2 || + (float)maxint[2] - (float)minint[2] >= INT_MAX-2) { + /* turning value in unsigned by subtracting minint + * would cause overflow + */ + fprintf(stderr, "(xdrfile error) Internal overflow compressing coordinates.\n"); + 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 = buf1; + luip = (unsigned int *) buf1; + smallidx = FIRSTIDX; + while (smallidx < LASTIDX && magicints[smallidx] < mindiff) { + smallidx++; + } + xdrfile_write_int(&smallidx,1,xfp); + tmp=smallidx+8; + maxidx = (LASTIDXtmp) ? FIRSTIDX : tmp; + smaller = magicints[tmp] / 2; + smallnum = 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]) < smallnum && + abs(thiscoord[1] - thiscoord[4]) < smallnum && + abs(thiscoord[2] - thiscoord[5]) < smallnum) { + /* 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) { + encodebits(buf2, bitsizeint[0], tmpcoord[0]); + encodebits(buf2, bitsizeint[1], tmpcoord[1]); + encodebits(buf2, bitsizeint[2], tmpcoord[2]); + } else { + encodeints(buf2, 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) { + tmpsum=0; + for(j=0;j<3;j++) { + tmp=thiscoord[j] - prevcoord[j]; + tmpsum+=tmp*tmp; + } + if (is_smaller == -1 && tmpsum >= smaller * smaller) { + is_smaller = 0; + } + + tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum; + tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum; + tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum; + + 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]) < smallnum && + abs(thiscoord[1] - prevcoord[1]) < smallnum && + abs(thiscoord[2] - prevcoord[2]) < smallnum) { + is_small = 1; + } + } + if (run != prevrun || is_smaller != 0) { + prevrun = run; + encodebits(buf2, 1, 1); /* flag the change in run-length */ + encodebits(buf2, 5, run+is_smaller+1); + } else { + encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */ + } + for (k=0; k < run; k+=3) { + encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]); + } + if (is_smaller != 0) { + smallidx += is_smaller; + if (is_smaller < 0) { + smallnum = smaller; + smaller = magicints[smallidx-1] / 2; + } else { + smaller = smallnum; + smallnum = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + } + } + if (buf2[1] != 0) buf2[0]++; + xdrfile_write_int(buf2,1,xfp); /* buf2[0] holds the length in bytes */ + tmp=xdrfile_write_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp); + if(tmp==(unsigned int)buf2[0]) + return size; + else + return -1; +} + + +/* Dont try do document Fortran interface, since + * Doxygen barfs at the F77_FUNC macro + */ +#ifndef DOXYGEN + +/************************************************************* + * Fortran77 interface for reading/writing portable data * + * The routine are not threadsafe when called from Fortran * + * (as they are when called from C) unless you compile with * + * this file with posix thread support. * + * Note that these are not multithread-safe. * + *************************************************************/ +#define MAX_FORTRAN_XDR 1024 +static XDRFILE *f77xdr[MAX_FORTRAN_XDR]; /* array of file handles */ +static int f77init = 1; /* zero array first time */ + +/* internal to this file: C<-->Fortran string conversion */ +static int ftocstr(char *dest, int dest_len, char *src, int src_len); +static int ctofstr(char *dest, int dest_len, char *src); + + +void +F77_FUNC(xdropen,XDROPEN)(int *fid, char *filename, char *mode, + int fn_len, int mode_len) +{ + char cfilename[512]; + char cmode[5]; + int i; + + /* zero array at first invocation */ + if(f77init) { + for(i=0;i= src && *p == ' ' ); + srclen = p - src + 1; + destlen--; + dest[0] = 0; + if (srclen > destlen) + return 1; + while (srclen--) + (*dest++ = *src++); + *dest = '\0'; + return 0; +} + + +static int ctofstr(char *dest, int destlen, char *src) +{ + while (destlen && *src) { + *dest++ = *src++; + destlen--; + } + while (destlen--) + *dest++ = ' '; + return 0; +} + + +void +F77_FUNC(xdrrstring,XDRRSTRING)(int *fid, char *str, int *ret, int len) +{ + char *cstr; + + if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) { + *ret = 0; + return; + } + if (ftocstr(cstr, len+1, str, len)) { + *ret = 0; + free(cstr); + return; + } + + *ret = xdrfile_read_string(cstr, len+1,f77xdr[*fid]); + ctofstr( str, len , cstr); + free(cstr); +} + +void +F77_FUNC(xdrwstring,XDRWSTRING)(int *fid, char *str, int *ret, int len) +{ + char *cstr; + + if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) { + *ret = 0; + return; + } + if (ftocstr(cstr, len+1, str, len)) { + *ret = 0; + free(cstr); + return; + } + + *ret = xdrfile_write_string(cstr, f77xdr[*fid]); + ctofstr( str, len , cstr); + free(cstr); +} + +void +F77_FUNC(xdrropaque,XDRROPAQUE)(int *fid, char *data, int *ndata, int *ret) +{ + *ret = xdrfile_read_opaque(data,*ndata,f77xdr[*fid]); +} + +void +F77_FUNC(xdrwopaque,XDRWOPAQUE)(int *fid, char *data, int *ndata, int *ret) +{ + *ret = xdrfile_write_opaque(data,*ndata,f77xdr[*fid]); +} + + +/* Write single-precision compressed 3d coordinates */ +void +F77_FUNC(xdrccs,XDRCCS)(int *fid, float *data, int *ncoord, + float *precision, int *ret) +{ + *ret = xdrfile_compress_coord_float(data,*ncoord,*precision,f77xdr[*fid]); +} + + +/* Read single-precision compressed 3d coordinates */ +void +F77_FUNC(xdrdcs,XDRDCS)(int *fid, float *data, int *ncoord, + float *precision, int *ret) +{ + *ret = xdrfile_decompress_coord_float(data,ncoord,precision,f77xdr[*fid]); +} + + +/* Write compressed 3d coordinates from double precision data */ +void +F77_FUNC(xdrccd,XDRCCD)(int *fid, double *data, int *ncoord, + double *precision, int *ret) +{ + *ret = xdrfile_compress_coord_double(data,*ncoord,*precision,f77xdr[*fid]); +} + +/* Read compressed 3d coordinates into double precision data */ +void +F77_FUNC(xddcd,XDRDCD)(int *fid, double *data, int *ncoord, + double *precision, int *ret) +{ + *ret = xdrfile_decompress_coord_double(data,ncoord,precision,f77xdr[*fid]); +} + + + + + + + +#endif /* DOXYGEN */ + +/************************************************************* + * End of higher-level routines - dont change things below! * + *************************************************************/ + + + + + + + + + + + + + + + + + + + + +/************************************************************* + * The rest of this file contains our own implementation * + * of the XDR calls in case you are compiling without them. * + * You do NOT want to change things here since it would make * + * things incompatible with the standard RPC/XDR routines. * + *************************************************************/ +#ifndef HAVE_RPC_XDR_H + +/* + * What follows is a modified version of the Sun XDR code. For reference + * we include their copyright and 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 + */ + +/* INT_MAX is defined in limits.h according to ANSI C */ +#if (INT_MAX > 2147483647) +# error Error: Cannot use builtin XDR support when size of int +# error is larger than 4 bytes. Use your system XDR libraries +# error instead, or modify the source code in xdrfile.c +#endif /* Check for 4 byte int type */ + + + + + +typedef int (*xdrproc_t) (XDR *, void *,...); + +#define xdr_getlong(xdrs, longp) \ + (*(xdrs)->x_ops->x_getlong)(xdrs, longp) +#define xdr_putlong(xdrs, longp) \ + (*(xdrs)->x_ops->x_putlong)(xdrs, longp) +#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 BYTES_PER_XDR_UNIT 4 +static char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0}; + +static int32_t +xdr_swapbytes(int32_t x) +{ + int32_t y,i; + char *px=(char *)&x; + char *py=(char *)&y; + + for(i=0;i<4;i++) + py[i]=px[3-i]; + + return y; +} + +static int32_t +xdr_htonl(int32_t x) +{ + int s=0x1234; + if( *((char *)&s)==(char)0x34) { + /* smallendian,swap bytes */ + return xdr_swapbytes(x); + } else { + /* bigendian, do nothing */ + return x; + } +} + +static int32_t +xdr_ntohl(int x) +{ + int s=0x1234; + if( *((char *)&s)==(char)0x34) { + /* smallendian, swap bytes */ + return xdr_swapbytes(x); + } else { + /* bigendian, do nothing */ + return x; + } +} + +static int +xdr_int (XDR *xdrs, int *ip) +{ + int32_t i32; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + i32 = (int32_t) *ip; + return xdr_putlong (xdrs, &i32); + + case XDR_DECODE: + if (!xdr_getlong (xdrs, &i32)) + { + return 0; + } + *ip = (int) i32; + case XDR_FREE: + return 1; + } + return 0; +} + +static int +xdr_u_int (XDR *xdrs, unsigned int *up) +{ + uint32_t ui32; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + ui32 = (uint32_t) * up; + return xdr_putlong (xdrs, (int32_t *)&ui32); + + case XDR_DECODE: + if (!xdr_getlong (xdrs, (int32_t *)&ui32)) + { + return 0; + } + *up = (uint32_t) ui32; + case XDR_FREE: + return 1; + } + return 0; +} + +static int +xdr_short (XDR *xdrs, short *sp) +{ + int32_t i32; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + i32 = (int32_t) *sp; + return xdr_putlong (xdrs, &i32); + + case XDR_DECODE: + if (!xdr_getlong (xdrs, &i32)) + { + return 0; + } + *sp = (short) i32; + return 1; + + case XDR_FREE: + return 1; + } + return 0; +} + +static int +xdr_u_short (XDR *xdrs, unsigned short *sp) +{ + uint32_t ui32; + + switch (xdrs->x_op) + { + case XDR_ENCODE: + ui32 = (uint32_t) *sp; + return xdr_putlong (xdrs, (int32_t *)&ui32); + + case XDR_DECODE: + if (!xdr_getlong (xdrs, (int32_t *)&ui32)) + { + return 0; + } + *sp = (unsigned short) ui32; + return 1; + + case XDR_FREE: + return 1; + } + return 0; +} + +static int +xdr_char (XDR *xdrs, char *cp) +{ + int i; + + i = (*cp); + if (!xdr_int (xdrs, &i)) + { + return 0; + } + *cp = i; + return 1; +} + +static int +xdr_u_char (XDR *xdrs, unsigned char *cp) +{ + unsigned int u; + + u = (*cp); + if (!xdr_u_int (xdrs, &u)) + { + return 0; + } + *cp = u; + return 1; +} + +/* + * 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. + */ +static int +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 1; + + /* + * 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 0; + } + if (rndup == 0) + return 1; + return xdr_getbytes (xdrs, (char *)crud, rndup); + + case XDR_ENCODE: + if (!xdr_putbytes (xdrs, cp, cnt)) + { + return 0; + } + if (rndup == 0) + return 1; + return xdr_putbytes (xdrs, xdr_zero, rndup); + + case XDR_FREE: + return 1; + } +#undef BYTES_PER_XDR_UNIT + return 0; +} + + +/* + * XDR null terminated ASCII strings + */ +static int +xdr_string (XDR *xdrs, char **cpp, unsigned int maxsize) +{ + char *sp = *cpp; /* sp is the actual string pointer */ + unsigned int size; + unsigned int nodesize; + + /* + * first deal with the length since xdr strings are counted-strings + */ + switch (xdrs->x_op) + { + case XDR_FREE: + if (sp == NULL) + { + return 1; /* already free */ + } + /* fall through... */ + case XDR_ENCODE: + if (sp == NULL) + return 0; + size = strlen (sp); + break; + case XDR_DECODE: + break; + } + if (!xdr_u_int (xdrs, &size)) + { + return 0; + } + if (size > maxsize) + { + return 0; + } + nodesize = size + 1; + + /* + * now deal with the actual bytes + */ + switch (xdrs->x_op) + { + case XDR_DECODE: + if (nodesize == 0) + { + return 1; + } + if (sp == NULL) + *cpp = sp = (char *) malloc (nodesize); + if (sp == NULL) + { + (void) fputs ("xdr_string: out of memory\n", stderr); + return 0; + } + sp[size] = 0; + /* fall into ... */ + + case XDR_ENCODE: + return xdr_opaque (xdrs, sp, size); + + case XDR_FREE: + free (sp); + *cpp = NULL; + return 1; + } + return 0; +} + + + +/* Floating-point stuff */ + +static int +xdr_float(XDR *xdrs, float *fp) +{ + switch (xdrs->x_op) { + + case XDR_ENCODE: + if (sizeof(float) == sizeof(int32_t)) + return (xdr_putlong(xdrs, (int32_t *)fp)); + else if (sizeof(float) == sizeof(int)) { + int32_t tmp = *(int *)fp; + return (xdr_putlong(xdrs, &tmp)); + } + break; + + case XDR_DECODE: + if (sizeof(float) == sizeof(int32_t)) + return (xdr_getlong(xdrs, (int32_t *)fp)); + else if (sizeof(float) == sizeof(int)) { + int32_t tmp; + if (xdr_getlong(xdrs, &tmp)) { + *(int *)fp = tmp; + return (1); + } + } + break; + + case XDR_FREE: + return (1); + } + return (0); +} + + +static int +xdr_double(XDR *xdrs, double *dp) +{ + /* Gromacs detects floating-point stuff at compile time, which is faster */ +#ifdef GROMACS +# ifndef FLOAT_FORMAT_IEEE754 +# error non-IEEE floating point system, or you defined GROMACS yourself... +# endif + int LSW; +# ifdef IEEE754_BIG_ENDIAN_WORD_ORDER + int LSW=1; +# else + int LSW=0; +# endif /* Big endian word order */ +#else + /* Outside Gromacs we rely on dynamic detection of FP order. */ + int LSW; /* Least significant fp word */ + + double x=0.987654321; /* Just a number */ + unsigned char ix = *((char *)&x); + + /* 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 + */ + 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 */ + fprintf(stderr,"Cannot detect floating-point word order.\n" + "Do you have a non-IEEE system?\n" + "Use system XDR libraries or fix xdr_double().\n"); + abort(); + } +#endif /* end of dynamic detection of fp word order */ + + switch (xdrs->x_op) { + + case XDR_ENCODE: + if (2*sizeof(int32_t) == sizeof(double)) { + int32_t *lp = (int32_t *)dp; + return (xdr_putlong(xdrs, lp+!LSW) && + xdr_putlong(xdrs, lp+LSW)); + } else if (2*sizeof(int) == sizeof(double)) { + int *ip = (int *)dp; + int32_t tmp[2]; + tmp[0] = ip[!LSW]; + tmp[1] = ip[LSW]; + return (xdr_putlong(xdrs, tmp) && + xdr_putlong(xdrs, tmp+1)); + } + break; + + case XDR_DECODE: + if (2*sizeof(int32_t) == sizeof(double)) { + int32_t *lp = (int32_t *)dp; + return (xdr_getlong(xdrs, lp+!LSW) && + xdr_getlong(xdrs, lp+LSW)); + } else if (2*sizeof(int) == sizeof(double)) { + int *ip = (int *)dp; + int32_t tmp[2]; + if (xdr_getlong(xdrs, tmp+!LSW) && + xdr_getlong(xdrs, tmp+LSW)) { + ip[0] = tmp[0]; + ip[1] = tmp[1]; + return (1); + } + } + break; + + case XDR_FREE: + return (1); + } + return (0); +} + + +static int xdrstdio_getlong (XDR *, int32_t *); +static int xdrstdio_putlong (XDR *, int32_t *); +static int xdrstdio_getbytes (XDR *, char *, unsigned int); +static int xdrstdio_putbytes (XDR *, char *, unsigned int); +static unsigned int xdrstdio_getpos (XDR *); +static int xdrstdio_setpos (XDR *, unsigned int); +static void xdrstdio_destroy (XDR *); + +/* + * Ops vector for stdio type XDR + */ +static const struct xdr_ops xdrstdio_ops = + { + xdrstdio_getlong, /* deserialize a long int */ + xdrstdio_putlong, /* serialize a long int */ + 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_destroy, /* destroy stream */ + }; + +/* + * Initialize a stdio xdr stream. + * Sets the xdr stream handle xdrs for use on the stream file. + * Operation flag is set to op. + */ +static void +xdrstdio_create (XDR *xdrs, FILE *file, enum xdr_op op) +{ + xdrs->x_op = op; + + xdrs->x_ops = (struct xdr_ops *) &xdrstdio_ops; + xdrs->x_private = (char *) file; +} + +/* + * 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 int +xdrstdio_getlong (XDR *xdrs, int32_t *lp) +{ + int32_t mycopy; + + if (fread ((char *) & mycopy, 4, 1, (FILE *) xdrs->x_private) != 1) + return 0; + *lp = (int32_t) xdr_ntohl (mycopy); + return 1; +} + +static int +xdrstdio_putlong (XDR *xdrs, int32_t *lp) +{ + int32_t mycopy = xdr_htonl (*lp); + lp = &mycopy; + if (fwrite ((char *) lp, 4, 1, (FILE *) xdrs->x_private) != 1) + return 0; + return 1; +} + +static int +xdrstdio_getbytes (XDR *xdrs, char *addr, unsigned int len) +{ + if ((len != 0) && (fread (addr, (int) len, 1, + (FILE *) xdrs->x_private) != 1)) + return 0; + return 1; +} + +static int +xdrstdio_putbytes (XDR *xdrs, char *addr, unsigned int len) +{ + if ((len != 0) && (fwrite (addr, (int) len, 1, + (FILE *) xdrs->x_private) != 1)) + return 0; + return 1; +} + +/* 32 bit fileseek operations */ +static unsigned int +xdrstdio_getpos (XDR *xdrs) +{ + return (unsigned int) ftell ((FILE *) xdrs->x_private); +} + +static int +xdrstdio_setpos (XDR *xdrs, unsigned int pos) +{ + return fseek ((FILE *) xdrs->x_private, pos, 0) < 0 ? 0 : 1; +} + + + +#endif /* HAVE_RPC_XDR_H not defined */ diff --git a/external/xdrfile/src/xdrfile_trr.c b/external/xdrfile/src/xdrfile_trr.c new file mode 100644 index 0000000..71e27e3 --- /dev/null +++ b/external/xdrfile/src/xdrfile_trr.c @@ -0,0 +1,533 @@ + /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "xdrfile.h" +#include "xdrfile_trr.h" + +#define BUFSIZE 128 +#define GROMACS_MAGIC 1993 + +typedef struct /* This struct describes the order and the */ +/* sizes of the structs in a trjfile, sizes are given in bytes. */ +{ + mybool bDouble; /* Double precision? */ + int ir_size; /* Backward compatibility */ + int e_size; /* Backward compatibility */ + int box_size; /* Non zero if a box is present */ + int vir_size; /* Backward compatibility */ + int pres_size; /* Backward compatibility */ + int top_size; /* Backward compatibility */ + int sym_size; /* Backward compatibility */ + int x_size; /* Non zero if coordinates are present */ + int v_size; /* Non zero if velocities are present */ + int f_size; /* Non zero if forces are present */ + + int natoms; /* The total number of atoms */ + int step; /* Current step number */ + int nre; /* Backward compatibility */ + float tf; /* Current time */ + float lambdaf; /* Current value of lambda */ + double td; /* Current time */ + double lambdad; /* Current value of lambda */ +} t_trnheader; + +static int nFloatSize(t_trnheader *sh,int *nflsz) +{ + int nflsize=0; + + if (sh->box_size) + nflsize = sh->box_size/(DIM*DIM); + else if (sh->x_size) + nflsize = sh->x_size/(sh->natoms*DIM); + else if (sh->v_size) + nflsize = sh->v_size/(sh->natoms*DIM); + else if (sh->f_size) + nflsize = sh->f_size/(sh->natoms*DIM); + else + return exdrHEADER; + + if (((nflsize != sizeof(float)) && (nflsize != sizeof(double)))) + return exdrHEADER; + + *nflsz = nflsize; + + return exdrOK; +} + +extern int do_trnheader(XDRFILE *xd,mybool bRead,t_trnheader *sh) +{ + int magic=GROMACS_MAGIC; + int nflsz,slen,result; + char *version = "GMX_trn_file"; + char buf[BUFSIZE]; + + if (xdrfile_read_int(&magic,1,xd) != 1) { + /* modification by RTM to return the right EOF code + this is what's happening in the XTC code */ + if (bRead) + return exdrENDOFFILE; + else + return exdrINT; + } + if (magic != GROMACS_MAGIC) + return exdrMAGIC; + if (bRead) + { + if (xdrfile_read_int(&slen,1,xd) != 1) + return exdrINT; + if (slen != strlen(version)+1) + return exdrSTRING; + if (xdrfile_read_string(buf,BUFSIZE,xd) <= 0) + return exdrSTRING; + } + else + { + slen = strlen(version)+1; + if (xdrfile_read_int(&slen,1,xd) != 1) + return exdrINT; + if (xdrfile_write_string(version,xd) != (strlen(version)+1) ) + return exdrSTRING; + } + if (xdrfile_read_int(&sh->ir_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->e_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->box_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->vir_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->pres_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->top_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->sym_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->x_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->v_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->f_size,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->natoms,1,xd) != 1) + return exdrINT; + + if ((result = nFloatSize(sh,&nflsz)) != exdrOK) + return result; + sh->bDouble = (nflsz == sizeof(double)); + + if (xdrfile_read_int(&sh->step,1,xd) != 1) + return exdrINT; + if (xdrfile_read_int(&sh->nre,1,xd) != 1) + return exdrINT; + if (sh->bDouble) + { + if (xdrfile_read_double(&sh->td,1,xd) != 1) + return exdrDOUBLE; + sh->tf = sh->td; + if (xdrfile_read_double(&sh->lambdad,1,xd) != 1) + return exdrDOUBLE; + sh->lambdaf = sh->lambdad; + } + else + { + if (xdrfile_read_float(&sh->tf,1,xd) != 1) + return exdrFLOAT; + sh->td = sh->tf; + if (xdrfile_read_float(&sh->lambdaf,1,xd) != 1) + return exdrFLOAT; + sh->lambdad = sh->lambdaf; + } + + return exdrOK; +} + +static int do_htrn(XDRFILE *xd,mybool bRead,t_trnheader *sh, + matrix box,rvec *x,rvec *v,rvec *f) +{ + double pvd[DIM*DIM]; + double *dx=NULL; + float pvf[DIM*DIM]; + float *fx=NULL; + int i,j; + + if (sh->bDouble) + { + if (sh->box_size != 0) + { + if (!bRead) + { + for(i=0; (ivir_size != 0) + { + if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM) + return exdrDOUBLE; + } + + if (sh->pres_size!= 0) + { + if (xdrfile_read_double(pvd,DIM*DIM,xd) != DIM*DIM) + return exdrDOUBLE; + } + + if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) { + dx = (double *)calloc(sh->natoms*DIM,sizeof(dx[0])); + if (NULL == dx) + return exdrNOMEM; + } + if (sh->x_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + if (bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jv_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + for(i=0; (inatoms); i++) + for(j=0; (jf_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + for(i=0; (inatoms); i++) + { + for(j=0; (jx_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) { + free(dx); + } + } + else + /* Float */ + { + if (sh->box_size != 0) + { + if (!bRead) + { + for(i=0; (ivir_size != 0) + { + if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM) + return exdrFLOAT; + } + + if (sh->pres_size!= 0) + { + if (xdrfile_read_float(pvf,DIM*DIM,xd) != DIM*DIM) + return exdrFLOAT; + } + + if ((sh->x_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) { + fx = (float *)calloc(sh->natoms*DIM,sizeof(fx[0])); + if (NULL == fx) + return exdrNOMEM; + } + if (sh->x_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + if (bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jv_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + for(i=0; (inatoms); i++) + for(j=0; (jf_size != 0) + { + if (!bRead) + { + for(i=0; (inatoms); i++) + for(j=0; (jnatoms*DIM,xd) == sh->natoms*DIM) + { + for(i=0; (inatoms); i++) + for(j=0; (jx_size != 0) || (sh->v_size != 0) || (sh->f_size != 0)) { + free(fx); + } + } + return exdrOK; +} + +static int do_trn(XDRFILE *xd,mybool bRead,int *step,float *t,float *lambda, + matrix box,int *natoms,rvec *x,rvec *v,rvec *f) +{ + t_trnheader *sh; + int result; + + sh = (t_trnheader *)calloc(1,sizeof(*sh)); + + if (!bRead) { + sh->box_size = (NULL != box) ? sizeof(matrix):0; + sh->x_size = ((NULL != x) ? (*natoms*sizeof(x[0])):0); + sh->v_size = ((NULL != v) ? (*natoms*sizeof(v[0])):0); + sh->f_size = ((NULL != f) ? (*natoms*sizeof(f[0])):0); + sh->natoms = *natoms; + sh->step = *step; + sh->nre = 0; + sh->td = *t; + sh->lambdad = *lambda; + sh->tf = *t; + sh->lambdaf = *lambda; + } + if ((result = do_trnheader(xd,bRead,sh)) != exdrOK) + return result; + + if (bRead) { + *natoms = sh->natoms; + *step = sh->step; + *t = sh->td; + *lambda = sh->lambdad; + } + if ((result = do_htrn(xd,bRead,sh,box,x,v,f)) != exdrOK) + return result; + + free(sh); + + return exdrOK; +} + +/************************************************************ + * + * The following routines are the exported ones + * + ************************************************************/ + +int read_trr_natoms(char *fn,int *natoms) +{ + XDRFILE *xd; + t_trnheader sh; + int result; + + xd = xdrfile_open(fn,"r"); + if (NULL == xd) + return exdrFILENOTFOUND; + if ((result = do_trnheader(xd,1,&sh)) != exdrOK) { + xdrfile_close(xd); + return result; + } + xdrfile_close(xd); + *natoms = sh.natoms; + + return exdrOK; +} + +int read_trr_nframes(char *fn, unsigned long *nframes) { + XDRFILE *xd; + int result, step; + float time, lambda; + int natoms; + matrix box; + rvec *x; + *nframes = 0; + + read_trr_natoms(fn, &natoms); + x = malloc(natoms * sizeof(*x)); + + xd = xdrfile_open(fn, "r"); + if (NULL == xd) + return exdrFILENOTFOUND; + + do { + result = read_trr(xd, natoms, &step, &time, &lambda, + box, x, NULL, NULL); + if (exdrENDOFFILE != result) { + (*nframes)++; + } + } while (result == exdrOK); + + xdrfile_close(xd); + free(x); + return exdrOK; +} + +int write_trr(XDRFILE *xd,int natoms,int step,float t,float lambda, + matrix box,rvec *x,rvec *v,rvec *f) +{ + return do_trn(xd,0,&step,&t,&lambda,box,&natoms,x,v,f); +} + +int read_trr(XDRFILE *xd,int natoms,int *step,float *t,float *lambda, + matrix box,rvec *x,rvec *v,rvec *f) +{ + return do_trn(xd,1,step,t,lambda,box,&natoms,x,v,f); +} + diff --git a/external/xdrfile/src/xdrfile_xtc.c b/external/xdrfile/src/xdrfile_xtc.c new file mode 100644 index 0000000..c0c5bdb --- /dev/null +++ b/external/xdrfile/src/xdrfile_xtc.c @@ -0,0 +1,164 @@ +/* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- + * + * $Id$ + * + * Copyright (c) 2009-2014, Erik Lindahl & David van der Spoel + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include "xdrfile.h" +#include "xdrfile_xtc.h" + +#define MAGIC 1995 + +enum { FALSE, TRUE }; + +static int xtc_header(XDRFILE *xd,int *natoms,int *step,float *time,mybool bRead) +{ + int result,magic,n=1; + + /* Note: read is same as write. He he he */ + magic = MAGIC; + if ((result = xdrfile_write_int(&magic,n,xd)) != n) + { + if (bRead) + return exdrENDOFFILE; + else + return exdrINT; + } + if (magic != MAGIC) + return exdrMAGIC; + if ((result = xdrfile_write_int(natoms,n,xd)) != n) + return exdrINT; + if ((result = xdrfile_write_int(step,n,xd)) != n) + return exdrINT; + if ((result = xdrfile_write_float(time,n,xd)) != n) + return exdrFLOAT; + + return exdrOK; +} + +static int xtc_coord(XDRFILE *xd,int *natoms,matrix box,rvec *x,float *prec, + mybool bRead) +{ + int i,j,result; + + /* box */ + result = xdrfile_read_float(box[0],DIM*DIM,xd); + if (DIM*DIM != result) + return exdrFLOAT; + else + { + if (bRead) + { + result = xdrfile_decompress_coord_float(x[0],natoms,prec,xd); + if (result != *natoms) + return exdr3DX; + } + else + { + result = xdrfile_compress_coord_float(x[0],*natoms,*prec,xd); + if (result != *natoms) + return exdr3DX; + } + } + return exdrOK; +} + +int read_xtc_natoms(char *fn,int *natoms) +{ + XDRFILE *xd; + int step,result; + float time; + + xd = xdrfile_open(fn,"r"); + if (NULL == xd) + return exdrFILENOTFOUND; + result = xtc_header(xd,natoms,&step,&time,TRUE); + xdrfile_close(xd); + + return result; +} + +int read_xtc_nframes(char* fn, unsigned long *nframes) { + XDRFILE *xd; + int result, step; + float time; + int natoms; + matrix box; + rvec *x; + float prec; + *nframes = 0; + + read_xtc_natoms(fn, &natoms); + x = malloc(natoms * sizeof(*x)); + + xd = xdrfile_open(fn, "r"); + if (NULL == xd) + return exdrFILENOTFOUND; + + do { + result = read_xtc(xd, natoms, &step, &time, box, x, &prec); + if (exdrENDOFFILE != result) { + (*nframes)++; + } + } while (result == exdrOK); + + xdrfile_close(xd); + free(x); + return exdrOK; +} + +int read_xtc(XDRFILE *xd, + int natoms,int *step,float *time, + matrix box,rvec *x,float *prec) +/* Read subsequent frames */ +{ + int result; + + if ((result = xtc_header(xd,&natoms,step,time,TRUE)) != exdrOK) + return result; + + if ((result = xtc_coord(xd,&natoms,box,x,prec,1)) != exdrOK) + return result; + + return exdrOK; +} + +int write_xtc(XDRFILE *xd, + int natoms,int step,float time, + matrix box,rvec *x,float prec) +/* Write a frame to xtc file */ +{ + int result; + + if ((result = xtc_header(xd,&natoms,&step,&time,FALSE)) != exdrOK) + return result; + + if ((result = xtc_coord(xd,&natoms,box,x,&prec,0)) != exdrOK) + return result; + + return exdrOK; +} diff --git a/src/c_abi/mod.rs b/src/c_abi/mod.rs new file mode 100644 index 0000000..3a71f91 --- /dev/null +++ b/src/c_abi/mod.rs @@ -0,0 +1,5 @@ +#[allow(non_upper_case_globals)] +pub mod xdrfile; +pub mod xdrfile_trr; +pub mod xdrfile_xtc; +pub mod xdr_seek; \ No newline at end of file diff --git a/src/c_abi/xdr_seek.rs b/src/c_abi/xdr_seek.rs new file mode 100644 index 0000000..a9db92e --- /dev/null +++ b/src/c_abi/xdr_seek.rs @@ -0,0 +1,160 @@ +use super::xdrfile::*; + +pub const _STDINT_H: u32 = 1; +pub const _FEATURES_H: u32 = 1; +pub const _DEFAULT_SOURCE: u32 = 1; +pub const __USE_ISOC11: u32 = 1; +pub const __USE_ISOC99: u32 = 1; +pub const __USE_ISOC95: u32 = 1; +pub const __USE_POSIX_IMPLICITLY: u32 = 1; +pub const _POSIX_SOURCE: u32 = 1; +pub const _POSIX_C_SOURCE: u32 = 200809; +pub const __USE_POSIX: u32 = 1; +pub const __USE_POSIX2: u32 = 1; +pub const __USE_POSIX199309: u32 = 1; +pub const __USE_POSIX199506: u32 = 1; +pub const __USE_XOPEN2K: u32 = 1; +pub const __USE_XOPEN2K8: u32 = 1; +pub const _ATFILE_SOURCE: u32 = 1; +pub const __USE_MISC: u32 = 1; +pub const __USE_ATFILE: u32 = 1; +pub const __USE_FORTIFY_LEVEL: u32 = 0; +pub const _STDC_PREDEF_H: u32 = 1; +pub const __STDC_IEC_559__: u32 = 1; +pub const __STDC_IEC_559_COMPLEX__: u32 = 1; +pub const __STDC_ISO_10646__: u32 = 201605; +pub const __STDC_NO_THREADS__: u32 = 1; +pub const __GNU_LIBRARY__: u32 = 6; +pub const __GLIBC__: u32 = 2; +pub const __GLIBC_MINOR__: u32 = 24; +pub const _SYS_CDEFS_H: u32 = 1; +pub const __WORDSIZE: u32 = 64; +pub const __WORDSIZE_TIME64_COMPAT32: u32 = 1; +pub const __SYSCALL_WORDSIZE: u32 = 64; +pub const _BITS_WCHAR_H: u32 = 1; +pub const INT8_MIN: i32 = -128; +pub const INT16_MIN: i32 = -32768; +pub const INT32_MIN: i32 = -2147483648; +pub const INT8_MAX: u32 = 127; +pub const INT16_MAX: u32 = 32767; +pub const INT32_MAX: u32 = 2147483647; +pub const UINT8_MAX: u32 = 255; +pub const UINT16_MAX: u32 = 65535; +pub const UINT32_MAX: u32 = 4294967295; +pub const INT_LEAST8_MIN: i32 = -128; +pub const INT_LEAST16_MIN: i32 = -32768; +pub const INT_LEAST32_MIN: i32 = -2147483648; +pub const INT_LEAST8_MAX: u32 = 127; +pub const INT_LEAST16_MAX: u32 = 32767; +pub const INT_LEAST32_MAX: u32 = 2147483647; +pub const UINT_LEAST8_MAX: u32 = 255; +pub const UINT_LEAST16_MAX: u32 = 65535; +pub const UINT_LEAST32_MAX: u32 = 4294967295; +pub const INT_FAST8_MIN: i32 = -128; +pub const INT_FAST16_MIN: i64 = -9223372036854775808; +pub const INT_FAST32_MIN: i64 = -9223372036854775808; +pub const INT_FAST8_MAX: u32 = 127; +pub const INT_FAST16_MAX: u64 = 9223372036854775807; +pub const INT_FAST32_MAX: u64 = 9223372036854775807; +pub const UINT_FAST8_MAX: u32 = 255; +pub const UINT_FAST16_MAX: i32 = -1; +pub const UINT_FAST32_MAX: i32 = -1; +pub const INTPTR_MIN: i64 = -9223372036854775808; +pub const INTPTR_MAX: u64 = 9223372036854775807; +pub const UINTPTR_MAX: i32 = -1; +pub const PTRDIFF_MIN: i64 = -9223372036854775808; +pub const PTRDIFF_MAX: u64 = 9223372036854775807; +pub const SIG_ATOMIC_MIN: i32 = -2147483648; +pub const SIG_ATOMIC_MAX: u32 = 2147483647; +pub const SIZE_MAX: i32 = -1; +pub const WINT_MIN: u32 = 0; +pub const WINT_MAX: u32 = 4294967295; +pub const DIM: u32 = 3; +pub type int_least8_t = ::std::os::raw::c_schar; +pub type int_least16_t = ::std::os::raw::c_short; +pub type int_least32_t = ::std::os::raw::c_int; +pub type int_least64_t = ::std::os::raw::c_long; +pub type uint_least8_t = ::std::os::raw::c_uchar; +pub type uint_least16_t = ::std::os::raw::c_ushort; +pub type uint_least32_t = ::std::os::raw::c_uint; +pub type uint_least64_t = ::std::os::raw::c_ulong; +pub type int_fast8_t = ::std::os::raw::c_schar; +pub type int_fast16_t = ::std::os::raw::c_long; +pub type int_fast32_t = ::std::os::raw::c_long; +pub type int_fast64_t = ::std::os::raw::c_long; +pub type uint_fast8_t = ::std::os::raw::c_uchar; +pub type uint_fast16_t = ::std::os::raw::c_ulong; +pub type uint_fast32_t = ::std::os::raw::c_ulong; +pub type uint_fast64_t = ::std::os::raw::c_ulong; +pub type intmax_t = ::std::os::raw::c_long; +pub type uintmax_t = ::std::os::raw::c_ulong; + +extern "C" { + pub fn xdr_tell(xd: *mut XDRFILE) -> i64; +} +extern "C" { + pub fn xdr_seek( + xd: *mut XDRFILE, + pos: i64, + whence: ::std::os::raw::c_int, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn xdr_flush(xd: *mut XDRFILE) -> ::std::os::raw::c_int; +} + + +#[cfg(test)] +mod tests { + + use super::*; + use super::super::xdrfile_xtc::*; + use std::ffi::CString; + + #[test] + fn test_xdr_tell() { + let path = CString::new("tests/1l2y.xtc").unwrap(); + let num_atoms = 304; + let mut time: f32 = 2.0; + let mut step: i32 = 5; + let box_vec = [[0.0; 3]; 3]; + let x_p = std::ptr::null_mut(); + let mut prec: f32 = 0.0; + + unsafe { + let mode = CString::new("r").unwrap(); + let xdr = xdrfile_open(path.as_ptr(), mode.as_ptr()); + assert!(!xdr.is_null()); + + let tell = xdr_tell(xdr); + assert!(tell == 0, "{}", tell); + + read_xtc(xdr, num_atoms, &mut step, &mut time, + box_vec.as_ptr() as *mut Matrix, + x_p, + &mut prec); + + let tell = xdr_tell(xdr); + assert!(tell > 0, "{}", tell); + } + } + + #[test] + fn test_xdr_seek() { + let path = CString::new("tests/1l2y.xtc").unwrap(); + + unsafe { + let mode = CString::new("r").unwrap(); + let xdr = xdrfile_open(path.as_ptr(), mode.as_ptr()); + assert!(!xdr.is_null()); + + let tell = xdr_tell(xdr); + assert!(tell == 0, "{}", tell); + + xdr_seek(xdr, 500, 0); + + let tell = xdr_tell(xdr); + assert!(tell == 500, "{}", tell); + } + } +} \ No newline at end of file diff --git a/src/c_abi/xdrfile.rs b/src/c_abi/xdrfile.rs new file mode 100644 index 0000000..e3721d3 --- /dev/null +++ b/src/c_abi/xdrfile.rs @@ -0,0 +1,520 @@ + +pub const DIM: u32 = 3; + +#[repr(C)] +#[derive(Debug, Copy, Clone)] +pub struct XDRFILE { + _unused: [u8; 0], +} + + +pub type BindgenTy1 = u32; +pub const exdrOK: BindgenTy1 = 0; +pub const exdrHEADER: BindgenTy1 = 1; +pub const exdrSTRING: BindgenTy1 = 2; +pub const exdrDOUBLE: BindgenTy1 = 3; +pub const exdrINT: BindgenTy1 = 4; +pub const exdrFLOAT: BindgenTy1 = 5; +pub const exdrUINT: BindgenTy1 = 6; +pub const exdr3DX: BindgenTy1 = 7; +pub const exdrCLOSE: BindgenTy1 = 8; +pub const exdrMAGIC: BindgenTy1 = 9; +pub const exdrNOMEM: BindgenTy1 = 10; +pub const exdrENDOFFILE: BindgenTy1 = 11; +pub const exdrFILENOTFOUND: BindgenTy1 = 12; +pub const exdrNR: BindgenTy1 = 13; + +extern "C" { + pub static mut exdr_message: [*mut ::std::os::raw::c_char; 13usize]; +} + +pub type Matrix = [[::std::os::raw::c_float; 3usize]; 3usize]; +pub type Rvec = [::std::os::raw::c_float; 3usize]; +pub type Mybool = ::std::os::raw::c_int; + +extern "C" { + #[doc = " \\brief Open a portable binary file, just like fopen()"] + #[doc = ""] + #[doc = " Use this routine much like calls to the standard library function"] + #[doc = " fopen(). The only difference is that the returned pointer should only"] + #[doc = " be used with routines defined in this header."] + #[doc = ""] + #[doc = " \\param path Full or relative path (including name) of the file"] + #[doc = " \\param mode \"r\" for reading, \"w\" for writing, \"a\" for append."] + #[doc = ""] + #[doc = " \\return Pointer to abstract xdr file datatype, or NULL if an error occurs."] + #[doc = ""] + pub fn xdrfile_open( + path: *const ::std::os::raw::c_char, + mode: *const ::std::os::raw::c_char, + ) -> *mut XDRFILE; +} +extern "C" { + #[doc = " \\brief Close a previously opened portable binary file, just like fclose()"] + #[doc = ""] + #[doc = " Use this routine much like calls to the standard library function"] + #[doc = " fopen(). The only difference is that it is used for an XDRFILE handle"] + #[doc = " instead of a FILE handle."] + #[doc = ""] + #[doc = " \\param xfp Pointer to an abstract XDRFILE datatype"] + #[doc = ""] + #[doc = " \\return 0 on success, non-zero on error."] + pub fn xdrfile_close(xfp: *mut XDRFILE) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a char type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of characters to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of characters read"] + pub fn xdrfile_read_char( + ptr: *mut ::std::os::raw::c_char, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a characters type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of characters to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of characters written"] + pub fn xdrfile_write_char( + ptr: *mut ::std::os::raw::c_char, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a unsigned \\a char type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of unsigned characters to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned characters read"] + pub fn xdrfile_read_uchar( + ptr: *mut ::std::os::raw::c_uchar, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a unsigned \\a characters type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of unsigned characters to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned characters written"] + pub fn xdrfile_write_uchar( + ptr: *mut ::std::os::raw::c_uchar, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a short type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of shorts to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of shorts read"] + pub fn xdrfile_read_short( + ptr: *mut ::std::os::raw::c_short, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a short type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of shorts to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of shorts written"] + pub fn xdrfile_write_short( + ptr: *mut ::std::os::raw::c_short, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a unsigned \\a short type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of unsigned shorts to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned shorts read"] + pub fn xdrfile_read_ushort( + ptr: *mut ::std::os::raw::c_ushort, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a unsigned \\a short type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of unsigned shorts to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned shorts written"] + pub fn xdrfile_write_ushort( + ptr: *mut ::std::os::raw::c_ushort, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a integer type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of integers to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of integers read"] + #[doc = ""] + #[doc = " The integer data type is assumed to be less than or equal to 32 bits."] + #[doc = ""] + #[doc = " We do not provide any routines for reading/writing 64-bit integers, since"] + #[doc = " - Not all XDR implementations support it"] + #[doc = " - Not all machines have 64-bit integers"] + #[doc = ""] + #[doc = " Split your 64-bit data into two 32-bit integers for portability!"] + pub fn xdrfile_read_int( + ptr: *mut ::std::os::raw::c_int, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a integer type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of integers to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of integers written"] + #[doc = ""] + #[doc = " The integer data type is assumed to be less than or equal to 32 bits."] + #[doc = ""] + #[doc = " We do not provide any routines for reading/writing 64-bit integers, since"] + #[doc = " - Not all XDR implementations support it"] + #[doc = " - Not all machines have 64-bit integers"] + #[doc = ""] + #[doc = " Split your 64-bit data into two 32-bit integers for portability!"] + pub fn xdrfile_write_int( + ptr: *mut ::std::os::raw::c_int, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a unsigned \\a integers type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of unsigned integers to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned integers read"] + #[doc = ""] + #[doc = " The integer data type is assumed to be less than or equal to 32 bits."] + #[doc = ""] + #[doc = " We do not provide any routines for reading/writing 64-bit integers, since"] + #[doc = " - Not all XDR implementations support it"] + #[doc = " - Not all machines have 64-bit integers"] + #[doc = ""] + #[doc = " Split your 64-bit data into two 32-bit integers for portability!"] + pub fn xdrfile_read_uint( + ptr: *mut ::std::os::raw::c_uint, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a unsigned \\a integer type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of unsigned integers to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of unsigned integers written"] + #[doc = ""] + #[doc = " The integer data type is assumed to be less than or equal to 32 bits."] + #[doc = ""] + #[doc = " We do not provide any routines for reading/writing 64-bit integers, since"] + #[doc = " - Not all XDR implementations support it"] + #[doc = " - Not all machines have 64-bit integers"] + #[doc = ""] + #[doc = " Split your 64-bit data into two 32-bit integers for portability!"] + pub fn xdrfile_write_uint( + ptr: *mut ::std::os::raw::c_uint, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a float type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of floats to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of floats read"] + pub fn xdrfile_read_float( + ptr: *mut ::std::os::raw::c_float, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a float type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of floats to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of floats written"] + pub fn xdrfile_write_float( + ptr: *mut ::std::os::raw::c_float, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read one or more \\a double type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param ndata Number of doubles to read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of doubles read"] + pub fn xdrfile_read_double( + ptr: *mut ::std::os::raw::c_double, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write one or more \\a double type variable(s)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param ndata Number of double to write."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of doubles written"] + pub fn xdrfile_write_double( + ptr: *mut ::std::os::raw::c_double, + ndata: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read a string (array of characters)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param maxlen Maximum length of string. If no end-of-string is encountered,"] + #[doc = " one byte less than this is read and end-of-string appended."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of characters read, including end-of-string"] + pub fn xdrfile_read_string( + ptr: *mut ::std::os::raw::c_char, + maxlen: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write a string (array of characters)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of characters written, including end-of-string"] + pub fn xdrfile_write_string( + ptr: *mut ::std::os::raw::c_char, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Read raw bytes from file (unknown datatype)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be written"] + #[doc = " \\param nbytes Number of bytes to read. No conversion whatsoever is done."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of bytes read from file"] + pub fn xdrfile_read_opaque( + ptr: *mut ::std::os::raw::c_char, + nbytes: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Write raw bytes to file (unknown datatype)"] + #[doc = ""] + #[doc = " \\param ptr Pointer to memory where data should be read"] + #[doc = " \\param nbytes Number of bytes to write. No conversion whatsoever is done."] + #[doc = " \\param xfp Handle to portable binary file, created with xdrfile_open()"] + #[doc = ""] + #[doc = " \\return Number of bytes written to file"] + pub fn xdrfile_write_opaque( + ptr: *mut ::std::os::raw::c_char, + nbytes: ::std::os::raw::c_int, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Compress coordiates in a float array to XDR file"] + #[doc = ""] + #[doc = " This routine will perform \\a lossy compression on the three-dimensional"] + #[doc = " coordinate data data specified and store it in the XDR file."] + #[doc = ""] + #[doc = " The lossy part of the compression consists of multiplying each"] + #[doc = " coordinate with the precision argument and then rounding to integers."] + #[doc = " We suggest a default value of 1000.0, which means you are guaranteed"] + #[doc = " three decimals of accuracy. The only limitation is that scaled coordinates"] + #[doc = " must still fit in an integer variable, so if the precision is 1000.0 the"] + #[doc = " coordinate magnitudes must be less than +-2e6."] + #[doc = ""] + #[doc = " \\param ptr Pointer to coordinates to compress (length 3*ncoord)"] + #[doc = " \\param ncoord Number of coordinate triplets in data"] + #[doc = " \\param precision Scaling factor for lossy compression. If it is <=0,"] + #[doc = " the default value of 1000.0 is used."] + #[doc = " \\param xfp Handle to portably binary file"] + #[doc = ""] + #[doc = " \\return Number of coordinate triplets written."] + #[doc = " IMPORTANT: Check that this is equal to ncoord - if it is"] + #[doc = " negative, an error occured. This should not happen with"] + #[doc = "\t \t normal data, but if your coordinates are NaN or very"] + #[doc = " large (>1e6) it is not possible to use the compression."] + #[doc = ""] + #[doc = " \\warning The compression algorithm is not part of the XDR standard,"] + #[doc = " and very complicated, so you will need this xdrfile module"] + #[doc = " to read it later."] + pub fn xdrfile_compress_coord_float( + ptr: *mut ::std::os::raw::c_float, + ncoord: ::std::os::raw::c_int, + precision: ::std::os::raw::c_float, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Decompress coordiates from XDR file to array of floats"] + #[doc = ""] + #[doc = " This routine will decompress three-dimensional coordinate data previously"] + #[doc = " stored in an XDR file and store it in the specified array of floats."] + #[doc = ""] + #[doc = " The precision used during the earlier compression is read from the file"] + #[doc = " and returned - you cannot adjust the accuracy at this stage."] + #[doc = ""] + #[doc = " \\param ptr Pointer to coordinates to compress (length>= 3*ncoord)"] + #[doc = " \\param ncoord Max number of coordinate triplets to read on input, actual"] + #[doc = " number of coordinate triplets read on return. If this"] + #[doc = " is smaller than the number of coordinates in the frame an"] + #[doc = " error will occur."] + #[doc = " \\param precision The precision used in the previous compression will be"] + #[doc = " written to this variable on return."] + #[doc = " \\param xfp Handle to portably binary file"] + #[doc = ""] + #[doc = " \\return Number of coordinate triplets read. If this is negative,"] + #[doc = " an error occured."] + #[doc = ""] + #[doc = " \\warning Since we cannot count on being able to set/get the"] + #[doc = " position of large files (>2Gb), it is not possible to"] + #[doc = " recover from errors by re-reading the frame if the"] + #[doc = " storage area you provided was too small. To avoid this"] + #[doc = " from happening, we recommend that you store the number of"] + #[doc = " coordinates triplet as an integer either in a header or"] + #[doc = " just before the compressed coordinate data, so you can"] + #[doc = " read it first and allocated enough memory."] + pub fn xdrfile_decompress_coord_float( + ptr: *mut ::std::os::raw::c_float, + ncoord: *mut ::std::os::raw::c_int, + precision: *mut ::std::os::raw::c_float, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Compress coordiates in a double array to XDR file"] + #[doc = ""] + #[doc = " This routine will perform \\a lossy compression on the three-dimensional"] + #[doc = " coordinate data data specified and store it in the XDR file. Double will"] + #[doc = " NOT give you any extra precision since the coordinates are compressed. This"] + #[doc = " routine just avoids allocating a temporary array of floats."] + #[doc = ""] + #[doc = " The lossy part of the compression consists of multiplying each"] + #[doc = " coordinate with the precision argument and then rounding to integers."] + #[doc = " We suggest a default value of 1000.0, which means you are guaranteed"] + #[doc = " three decimals of accuracy. The only limitation is that scaled coordinates"] + #[doc = " must still fit in an integer variable, so if the precision is 1000.0 the"] + #[doc = " coordinate magnitudes must be less than +-2e6."] + #[doc = ""] + #[doc = " \\param ptr Pointer to coordinates to compress (length 3*ncoord)"] + #[doc = " \\param ncoord Number of coordinate triplets in data"] + #[doc = " \\param precision Scaling factor for lossy compression. If it is <=0, the"] + #[doc = " default value of 1000.0 is used."] + #[doc = " \\param xfp Handle to portably binary file"] + #[doc = ""] + #[doc = " \\return Number of coordinate triplets written."] + #[doc = " IMPORTANT: Check that this is equal to ncoord - if it is"] + #[doc = " negative, an error occured. This should not happen with"] + #[doc = " normal data, but if your coordinates are NaN or very"] + #[doc = " large (>1e6) it is not possible to use the compression."] + #[doc = ""] + #[doc = " \\warning The compression algorithm is not part of the XDR standard,"] + #[doc = " and very complicated, so you will need this xdrfile module"] + #[doc = " to read it later."] + pub fn xdrfile_compress_coord_double( + ptr: *mut ::std::os::raw::c_double, + ncoord: ::std::os::raw::c_int, + precision: ::std::os::raw::c_double, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} +extern "C" { + #[doc = " \\brief Decompress coordiates from XDR file to array of doubles"] + #[doc = ""] + #[doc = " This routine will decompress three-dimensional coordinate data previously"] + #[doc = " stored in an XDR file and store it in the specified array of doubles."] + #[doc = " Double will NOT give you any extra precision since the coordinates are"] + #[doc = " compressed. This routine just avoids allocating a temporary array of floats."] + #[doc = ""] + #[doc = " The precision used during the earlier compression is read from the file"] + #[doc = " and returned - you cannot adjust the accuracy at this stage."] + #[doc = ""] + #[doc = " \\param ptr Pointer to coordinates to compress (length>= 3*ncoord)"] + #[doc = " \\param ncoord Max number of coordinate triplets to read on input, actual"] + #[doc = " number of coordinate triplets read on return. If this"] + #[doc = " is smaller than the number of coordinates in the frame an"] + #[doc = " error will occur."] + #[doc = " \\param precision The precision used in the previous compression will be"] + #[doc = " written to this variable on return."] + #[doc = " \\param xfp Handle to portably binary file"] + #[doc = ""] + #[doc = " \\return Number of coordinate triplets read. If this is negative,"] + #[doc = " an error occured."] + #[doc = ""] + #[doc = " \\warning Since we cannot count on being able to set/get the"] + #[doc = " position of large files (>2Gb), it is not possible to"] + #[doc = " recover from errors by re-reading the frame if the"] + #[doc = " storage area you provided was too small. To avoid this"] + #[doc = " from happening, we recommend that you store the number of"] + #[doc = " coordinates triplet as an integer either in a header or"] + #[doc = " just before the compressed coordinate data, so you can"] + #[doc = " read it first and allocated enough memory."] + pub fn xdrfile_decompress_coord_double( + ptr: *mut ::std::os::raw::c_double, + ncoord: *mut ::std::os::raw::c_int, + precision: *mut ::std::os::raw::c_double, + xfp: *mut XDRFILE, + ) -> ::std::os::raw::c_int; +} diff --git a/src/c_abi/xdrfile_trr.rs b/src/c_abi/xdrfile_trr.rs new file mode 100644 index 0000000..d80effb --- /dev/null +++ b/src/c_abi/xdrfile_trr.rs @@ -0,0 +1,133 @@ +use super::xdrfile::*; + + +extern "C" { + pub fn read_trr_natoms( + fn_: *const ::std::os::raw::c_char, + natoms: *const ::std::os::raw::c_int, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn read_trr_nframes( + fn_: *const ::std::os::raw::c_char, + nframes: *const ::std::os::raw::c_ulong, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn read_trr( + xd: *mut XDRFILE, + natoms: ::std::os::raw::c_int, + step: *mut ::std::os::raw::c_int, + t: *mut ::std::os::raw::c_float, + lambda: *mut ::std::os::raw::c_float, + box_vec: *mut Matrix, + x: *mut Rvec, + v: *mut Rvec, + f: *mut Rvec, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn write_trr( + xd: *mut XDRFILE, + natoms: ::std::os::raw::c_int, + step: ::std::os::raw::c_int, + t: ::std::os::raw::c_float, + lambda: ::std::os::raw::c_float, + box_vec: *const Matrix, + x: *const Rvec, + v: *const Rvec, + f: *const Rvec, + ) -> ::std::os::raw::c_int; +} + +#[cfg(test)] +mod tests { + + use super::*; + use std::ffi::CString; + use tempfile::NamedTempFile; + + #[test] + fn test_read_trr_natoms() { + let path = CString::new("tests/1l2y.trr").unwrap(); + + let mut natoms = 0; + unsafe { + read_trr_natoms(path.as_ptr() as *const i8, &mut natoms); + } + assert!(natoms == 304); + } + + #[test] + fn test_read_trr_nframes() { + let path = CString::new("tests/1l2y.trr").unwrap(); + let mut nframes: u64 = 0; + + unsafe { + let code = read_trr_nframes(path.as_ptr() as *const i8, &mut nframes); + assert!(code as u32 == exdrOK); + } + assert!(nframes == 38, "{:?}", nframes); + } + + #[test] + fn test_read_write_trr() { + let tempfile = NamedTempFile::new().unwrap(); + let tmp_path = CString::new(tempfile.path().to_str().unwrap()).unwrap(); + + // write atoms to tempfile + let natoms: i32 = 2; + let time: f32 = 2.0; + let lambda: f32 = 1.0; + let step: i32 = 5; + + let box_vec: Matrix = [[1.0, 2.0, 3.0], [2.0, 1.0, 3.0], [3.0, 2.0, 1.0]]; + let x: Vec = vec![[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]; + let v: Vec = vec![[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]; + let f: Vec = vec![[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]; + + unsafe { + let mode = CString::new("w").unwrap(); + let xdr = xdrfile_open(tmp_path.as_ptr(), mode.as_ptr()); + let write_code = write_trr(xdr, natoms, step, time, lambda, + box_vec.as_ptr() as *mut Matrix, + x.as_ptr() as *mut Rvec, + v.as_ptr() as *mut Rvec, + f.as_ptr() as *mut Rvec); + assert!(write_code as u32 == exdrOK); + xdrfile_close(xdr); + } + + // read atoms from tempfile + let mut time2: f32 = 0.0; + let mut lambda2: f32 = 0.0; + let mut step2: i32 = 0; + + let box_vec2: Matrix = [[0.0, 0.0, 0.0]; 3]; + let x2: Vec = vec![[0.0, 0.0, 0.0]; 2]; + let v2: Vec = vec![[0.0, 0.0, 0.0]; 2]; + let f2: Vec = vec![[0.0, 0.0, 0.0]; 2]; + + unsafe { + let mode = CString::new("r").unwrap(); + let xdr = xdrfile_open(tmp_path.as_ptr(), mode.as_ptr()); + let read_code = read_trr(xdr, natoms, &mut step2, &mut time2, + &mut lambda2, box_vec2.as_ptr() as *mut Matrix, + x2.as_ptr() as *mut Rvec, + v2.as_ptr() as *mut Rvec, + f2.as_ptr() as *mut Rvec); + assert!(read_code as u32 == exdrOK); + xdrfile_close(xdr); + } + + // make sure everything is still the same + assert!(step2 == step); + assert!(time2 == time); + assert!(lambda == lambda2); + assert!(box_vec2 == box_vec); + assert!(x2 == x); + assert!(v2 == v); + assert!(f2 == f); + } + +} \ No newline at end of file diff --git a/src/c_abi/xdrfile_xtc.rs b/src/c_abi/xdrfile_xtc.rs new file mode 100644 index 0000000..e074d51 --- /dev/null +++ b/src/c_abi/xdrfile_xtc.rs @@ -0,0 +1,114 @@ +use super::xdrfile::*; + + +extern "C" { + pub fn read_xtc_natoms( + fn_: *const ::std::os::raw::c_char, + natoms: *const ::std::os::raw::c_int, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn read_xtc_nframes( + fn_: *const ::std::os::raw::c_char, + nframes: *const ::std::os::raw::c_ulong, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn read_xtc( + xd: *mut XDRFILE, + natoms: ::std::os::raw::c_int, + step: *mut ::std::os::raw::c_int, + time: *mut ::std::os::raw::c_float, + box_vec: *mut Matrix, + x: *mut Rvec, + prec: *mut ::std::os::raw::c_float, + ) -> ::std::os::raw::c_int; +} +extern "C" { + pub fn write_xtc( + xd: *mut XDRFILE, + natoms: ::std::os::raw::c_int, + step: ::std::os::raw::c_int, + time: ::std::os::raw::c_float, + box_vec: *mut Matrix, + x: *mut Rvec, + prec: ::std::os::raw::c_float, + ) -> ::std::os::raw::c_int; +} + +#[cfg(test)] +mod tests { + + use super::*; + use std::ffi::CString; + use tempfile::NamedTempFile; + + #[test] + fn test_read_xtc_natoms() { + let path = CString::new("tests/1l2y.xtc").unwrap(); + + let mut natoms = 0; + unsafe { + read_xtc_natoms(path.as_ptr() as *mut i8, &mut natoms); + } + assert!(natoms == 304); + } + + #[test] + fn test_read_xtc_nframes() { + let path = CString::new("tests/1l2y.xtc").unwrap(); + let mut nframes: u64 = 0; + + unsafe { + let code = read_xtc_nframes(path.as_ptr() as *const i8, &mut nframes); + assert!(code as u32 == exdrOK); + } + assert!(nframes == 38, "{:?}", nframes); + } + + #[test] + fn test_read_write_xtc() { + let tempfile = NamedTempFile::new().unwrap(); + let tmp_path = CString::new(tempfile.path().to_str().unwrap()).unwrap(); + + // write atoms to tempfile + let natoms: i32 = 2; + let time: f32 = 2.0; + let step: i32 = 5; + let box_vec: Matrix = [[1.0, 2.0, 3.0], [2.0, 1.0, 3.0], [3.0, 2.0, 1.0]]; + let x: Vec = vec![[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]; + + unsafe { + let mode = CString::new("w").unwrap(); + let xdr = xdrfile_open(tmp_path.as_ptr(), mode.as_ptr()); + let write_code = write_xtc(xdr, natoms, step, time, + box_vec.as_ptr() as *mut Matrix, x.as_ptr() as *mut Rvec, + 1000.0); + assert!(write_code as u32 == exdrOK); + xdrfile_close(xdr); + } + + // read atoms from tempfile + let mut time2: f32 = 0.0; + let mut step2: i32 = 0; + let box_vec2: Matrix = [[0.0, 0.0, 0.0]; 3]; + let x2: Vec = vec![[0.0, 0.0, 0.0]; 2]; + let mut prec: f32 = 0.0; + + unsafe { + let mode = CString::new("r").unwrap(); + let xdr = xdrfile_open(tmp_path.as_ptr(), mode.as_ptr()); + let read_code = read_xtc(xdr, natoms, &mut step2, &mut time2, + box_vec2.as_ptr() as *mut Matrix, x2.as_ptr() as *mut Rvec, + &mut prec); + assert!(read_code as u32 == exdrOK); + xdrfile_close(xdr); + } + + // make sure everything is still the same + assert!(step2 == step); + assert!(time2 == time); + assert!(box_vec2 == box_vec); + assert!(x2 == x); + } +} \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..268972c --- /dev/null +++ b/src/lib.rs @@ -0,0 +1 @@ +pub mod c_abi; diff --git a/tests/1l2y.trr b/tests/1l2y.trr new file mode 100644 index 0000000..89dc66b Binary files /dev/null and b/tests/1l2y.trr differ diff --git a/tests/1l2y.xtc b/tests/1l2y.xtc new file mode 100644 index 0000000..8855991 Binary files /dev/null and b/tests/1l2y.xtc differ