--- /dev/null
- if((APPLE OR CYGWIN OR ${CMAKE_SYSTEM_NAME} MATCHES "Linux|.*BSD") AND NOT GMX_BUILD_MDRUN_ONLY)
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+# Manage the Gromacs shared library setup.
+
+########################################################################
+# Shared/static library settings
+########################################################################
+# Determine the defaults (this block has no effect if the variables have
+# already been set)
++if((APPLE OR CYGWIN OR ${CMAKE_SYSTEM_NAME} MATCHES "Linux|.*BSD|GNU") AND NOT GMX_BUILD_MDRUN_ONLY)
+ # Maybe Solaris should be here? Patch this if you know!
+ SET(SHARED_LIBS_DEFAULT ON)
+elseif(WIN32 OR ${CMAKE_SYSTEM_NAME} MATCHES "BlueGene")
+ # Support for shared libs on native Windows is a bit new. Its
+ # default might change later if/when we sort things out. Also,
+ # Cray should go here. What variable value can detect it?
+ SET(SHARED_LIBS_DEFAULT OFF)
+else()
+ if (NOT DEFINED BUILD_SHARED_LIBS)
+ message(STATUS "Defaulting to building static libraries")
+ endif()
+ SET(SHARED_LIBS_DEFAULT OFF)
+endif()
+if (GMX_PREFER_STATIC_LIBS)
+ if (NOT DEFINED BUILD_SHARED_LIBS AND SHARED_LIBS_DEFAULT)
+ message("Searching for static libraries requested, so the GROMACS libraries will also be static (BUILD_SHARED_LIBS=OFF)")
+ endif()
+ set(SHARED_LIBS_DEFAULT OFF)
+endif()
+set(GMX_PREFER_STATIC_LIBS_DEFAULT OFF)
+if (WIN32 AND NOT CYGWIN AND NOT BUILD_SHARED_LIBS)
+ set(GMX_PREFER_STATIC_LIBS_DEFAULT ON)
+endif()
+
+# Declare the user-visible options
+option(BUILD_SHARED_LIBS "Enable shared libraries (can be problematic e.g. with MPI, or on some HPC systems)" ${SHARED_LIBS_DEFAULT})
+if(BUILD_SHARED_LIBS AND GMX_BUILD_MDRUN_ONLY)
+ message(WARNING "Both BUILD_SHARED_LIBS and GMX_BUILD_MDRUN_ONLY are set. Generally, an mdrun-only build should prefer to use static libraries, which is the default if you make a fresh build tree. You may be re-using an old build tree, and so may wish to set BUILD_SHARED_LIBS=off yourself.")
+endif()
+
+if (UNIX)
+ set(GMX_PREFER_STATIC_LIBS_DESCRIPTION
+ "When finding libraries prefer static archives (it will only work if static versions of external dependencies are available and found)")
+elseif (WIN32 AND NOT CYGWIN)
+ set(GMX_PREFER_STATIC_LIBS_DESCRIPTION
+ "When finding libraries prefer static system libraries (MT instead of MD)")
+endif()
+option(GMX_PREFER_STATIC_LIBS "${GMX_PREFER_STATIC_LIBS_DESCRIPTION}"
+ ${GMX_PREFER_STATIC_LIBS_DEFAULT})
+mark_as_advanced(GMX_PREFER_STATIC_LIBS)
+
+# Act on the set values
+if (UNIX AND GMX_PREFER_STATIC_LIBS)
+ if(BUILD_SHARED_LIBS)
+ # Warn the user about the combination. But don't overwrite the request.
+ message(WARNING "Searching for static libraries requested, and building shared Gromacs libraries requested. This might cause problems linking later.")
+ endif()
+ # On Linux .a is the static library suffix, on Mac OS X .lib can also
+ # be used, so we'll add both to the preference list.
+ SET(CMAKE_FIND_LIBRARY_SUFFIXES ".lib;.a" ${CMAKE_FIND_LIBRARY_SUFFIXES})
+endif()
+
+# ==========
+# Only things for managing shared libraries and build types on Windows follow
+
+# Change the real CMake variables so we prefer static linking. This
+# should be a function so we can have proper local variables while
+# avoiding duplicating code.
+function(gmx_manage_prefer_static_libs_flags build_type)
+ if("${build_type}" STREQUAL "")
+ set(punctuation "") # for general compiler flags (e.g.) CMAKE_CXX_FLAGS
+ else()
+ set(punctuation "_") # for build-type-specific compiler flags (e.g.) CMAKE_CXX_FLAGS_RELEASE
+ endif()
+
+ # Change the real CMake variables for the given build type in each
+ # language, in the parent scope.
+ foreach(language C CXX)
+ string(REPLACE /MD /MT CMAKE_${language}_FLAGS${punctuation}${build_type} ${CMAKE_${language}_FLAGS${punctuation}${build_type}} PARENT_SCOPE)
+ endforeach()
+endfunction()
+
+IF( WIN32 AND NOT CYGWIN)
+ if (NOT BUILD_SHARED_LIBS)
+ if(NOT GMX_PREFER_STATIC_LIBS)
+ message(WARNING "Shared system libraries requested, and static Gromacs libraries requested.")
+ endif()
+ else()
+ message(FATAL_ERROR "BUILD_SHARED_LIBS=ON not yet working for Windows in the master branch")
+ if(GMX_PREFER_STATIC_LIBS)
+ #this combination segfaults (illegal passing of file handles)
+ message(FATAL_ERROR "Static system libraries requested, and shared Gromacs libraries requested.")
+ endif()
+ add_definitions(-DUSE_VISIBILITY -DTMPI_USE_VISIBILITY)
+ set(PKG_CFLAGS "$PKG_CFLAGS -DUSE_VISIBILITY -DTMPI_USE_VISIBILITY")
+ endif()
+
+ IF (GMX_PREFER_STATIC_LIBS)
+ foreach(build_type "" ${build_types_with_explicit_flags})
+ gmx_manage_prefer_static_libs_flags("${build_type}")
+ endforeach()
+ ENDIF()
+ IF( CMAKE_C_COMPILER_ID MATCHES "Intel" )
+ if(BUILD_SHARED_LIBS) #not sure why incremental building with shared libs doesn't work
+ STRING(REPLACE "/INCREMENTAL:YES" "" CMAKE_SHARED_LINKER_FLAGS ${CMAKE_SHARED_LINKER_FLAGS})
+ endif()
+ ENDIF()
+ENDIF()
--- /dev/null
-
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "xdrf.h"
+#include "xdr_datatype.h"
+#include "futil.h"
+
+/* This is just for clarity - it can never be anything but 4! */
+#define XDR_INT_SIZE 4
+
+/* same order as the definition of xdr_datatype */
+const char *xdr_datatype_names[] =
+{
+ "int",
+ "float",
+ "double",
+ "large int",
+ "char",
+ "string"
+};
+
+
+/*___________________________________________________________________________
+ |
+ | what follows are the C routine to read/write compressed coordinates together
+ | with some routines to assist in this task (those are marked
+ | static and cannot be called from user programs)
+ */
+#define MAXABS INT_MAX-2
+
+#ifndef MIN
+#define MIN(x, y) ((x) < (y) ? (x) : (y))
+#endif
+#ifndef MAX
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+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))
+
+
+/*____________________________________________________________________________
+ |
+ | sendbits - encode num into buf using the specified number of bits
+ |
+ | This routines appends the value of num to the bits already present in
+ | the array buf. You need to give it the number of bits to use and you
+ | better make sure that this number of bits is enough to hold the value
+ | Also num must be positive.
+ |
+ */
+
+static void sendbits(int buf[], int num_of_bits, int num)
+{
+
+ unsigned int cnt, lastbyte;
+ int lastbits;
+ unsigned char * cbuf;
+
+ cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
+ cnt = (unsigned int) buf[0];
+ lastbits = buf[1];
+ lastbyte = (unsigned int) buf[2];
+ while (num_of_bits >= 8)
+ {
+ lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
+ cbuf[cnt++] = lastbyte >> lastbits;
+ num_of_bits -= 8;
+ }
+ if (num_of_bits > 0)
+ {
+ lastbyte = (lastbyte << num_of_bits) | num;
+ lastbits += num_of_bits;
+ if (lastbits >= 8)
+ {
+ lastbits -= 8;
+ cbuf[cnt++] = lastbyte >> lastbits;
+ }
+ }
+ buf[0] = cnt;
+ buf[1] = lastbits;
+ buf[2] = lastbyte;
+ if (lastbits > 0)
+ {
+ cbuf[cnt] = lastbyte << (8 - lastbits);
+ }
+}
+
+/*_________________________________________________________________________
+ |
+ | sizeofint - calculate bitsize of an integer
+ |
+ | return the number of bits needed to store an integer with given max size
+ |
+ */
+
+static int sizeofint(const int size)
+{
+ int num = 1;
+ int num_of_bits = 0;
+
+ while (size >= num && num_of_bits < 32)
+ {
+ num_of_bits++;
+ num <<= 1;
+ }
+ return num_of_bits;
+}
+
+/*___________________________________________________________________________
+ |
+ | sizeofints - calculate 'bitsize' of compressed ints
+ |
+ | given the number of small unsigned integers and the maximum value
+ | return the number of bits needed to read or write them with the
+ | routines receiveints and sendints. You need this parameter when
+ | calling these routines. Note that for many calls I can use
+ | the variable 'smallidx' which is exactly the number of bits, and
+ | So I don't need to call 'sizeofints for those calls.
+ */
+
+static int sizeofints( const int num_of_ints, unsigned int sizes[])
+{
+ int i, num;
+ int bytes[32];
+ unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
+ num_of_bytes = 1;
+ bytes[0] = 1;
+ num_of_bits = 0;
+ for (i = 0; i < num_of_ints; i++)
+ {
+ tmp = 0;
+ for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
+ {
+ tmp = bytes[bytecnt] * sizes[i] + tmp;
+ bytes[bytecnt] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ while (tmp != 0)
+ {
+ bytes[bytecnt++] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ num_of_bytes = bytecnt;
+ }
+ num = 1;
+ num_of_bytes--;
+ while (bytes[num_of_bytes] >= num)
+ {
+ num_of_bits++;
+ num *= 2;
+ }
+ return num_of_bits + num_of_bytes * 8;
+
+}
+
+/*____________________________________________________________________________
+ |
+ | sendints - send a small set of small integers in compressed format
+ |
+ | this routine is used internally by xdr3dfcoord, to send a set of
+ | small integers to the buffer.
+ | Multiplication with fixed (specified maximum ) sizes is used to get
+ | to one big, multibyte integer. Allthough the routine could be
+ | modified to handle sizes bigger than 16777216, or more than just
+ | a few integers, this is not done, because the gain in compression
+ | isn't worth the effort. Note that overflowing the multiplication
+ | or the byte buffer (32 bytes) is unchecked and causes bad results.
+ |
+ */
+
+static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
+ unsigned int sizes[], unsigned int nums[])
+{
+
+ int i, num_of_bytes, bytecnt;
+ unsigned int bytes[32], tmp;
+
+ tmp = nums[0];
+ num_of_bytes = 0;
+ do
+ {
+ bytes[num_of_bytes++] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ while (tmp != 0);
+
+ for (i = 1; i < num_of_ints; i++)
+ {
+ if (nums[i] >= sizes[i])
+ {
+ fprintf(stderr, "major breakdown in sendints num %u doesn't "
+ "match size %u\n", nums[i], sizes[i]);
+ exit(1);
+ }
+ /* use one step multiply */
+ tmp = nums[i];
+ for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
+ {
+ tmp = bytes[bytecnt] * sizes[i] + tmp;
+ bytes[bytecnt] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ while (tmp != 0)
+ {
+ bytes[bytecnt++] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ num_of_bytes = bytecnt;
+ }
+ if (num_of_bits >= num_of_bytes * 8)
+ {
+ for (i = 0; i < num_of_bytes; i++)
+ {
+ sendbits(buf, 8, bytes[i]);
+ }
+ sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
+ }
+ else
+ {
+ for (i = 0; i < num_of_bytes-1; i++)
+ {
+ sendbits(buf, 8, bytes[i]);
+ }
+ sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
+ }
+}
+
+
+/*___________________________________________________________________________
+ |
+ | receivebits - decode number from buf using specified number of bits
+ |
+ | extract the number of bits from the array buf and construct an integer
+ | from it. Return that value.
+ |
+ */
+
+static int receivebits(int buf[], int num_of_bits)
+{
+
+ int cnt, num, lastbits;
+ unsigned int lastbyte;
+ unsigned char * cbuf;
+ int mask = (1 << num_of_bits) -1;
+
+ cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
+ cnt = buf[0];
+ lastbits = (unsigned int) buf[1];
+ lastbyte = (unsigned int) buf[2];
+
+ num = 0;
+ while (num_of_bits >= 8)
+ {
+ lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
+ num |= (lastbyte >> lastbits) << (num_of_bits - 8);
+ num_of_bits -= 8;
+ }
+ if (num_of_bits > 0)
+ {
+ if (lastbits < num_of_bits)
+ {
+ lastbits += 8;
+ lastbyte = (lastbyte << 8) | cbuf[cnt++];
+ }
+ lastbits -= num_of_bits;
+ num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
+ }
+ num &= mask;
+ buf[0] = cnt;
+ buf[1] = lastbits;
+ buf[2] = lastbyte;
+ return num;
+}
+
+/*____________________________________________________________________________
+ |
+ | receiveints - decode 'small' integers from the buf array
+ |
+ | this routine is the inverse from sendints() and decodes the small integers
+ | written to buf by calculating the remainder and doing divisions with
+ | the given sizes[]. You need to specify the total number of bits to be
+ | used from buf in num_of_bits.
+ |
+ */
+
+static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
+ unsigned int sizes[], int nums[])
+{
+ int bytes[32];
+ int i, j, num_of_bytes, p, num;
+
+ bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
+ num_of_bytes = 0;
+ while (num_of_bits > 8)
+ {
+ bytes[num_of_bytes++] = receivebits(buf, 8);
+ num_of_bits -= 8;
+ }
+ if (num_of_bits > 0)
+ {
+ bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
+ }
+ for (i = num_of_ints-1; i > 0; i--)
+ {
+ num = 0;
+ for (j = num_of_bytes-1; j >= 0; j--)
+ {
+ num = (num << 8) | bytes[j];
+ p = num / sizes[i];
+ bytes[j] = p;
+ num = num - p * sizes[i];
+ }
+ nums[i] = num;
+ }
+ nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
+}
+
+/*____________________________________________________________________________
+ |
+ | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
+ |
+ | this routine reads or writes (depending on how you opened the file with
+ | xdropen() ) a large number of 3d coordinates (stored in *fp).
+ | The number of coordinates triplets to write is given by *size. On
+ | read this number may be zero, in which case it reads as many as were written
+ | or it may specify the number if triplets to read (which should match the
+ | number written).
+ | Compression is achieved by first converting all floating numbers to integer
+ | using multiplication by *precision and rounding to the nearest integer.
+ | Then the minimum and maximum value are calculated to determine the range.
+ | The limited range of integers so found, is used to compress the coordinates.
+ | In addition the differences between succesive coordinates is calculated.
+ | If the difference happens to be 'small' then only the difference is saved,
+ | compressing the data even more. The notion of 'small' is changed dynamically
+ | and is enlarged or reduced whenever needed or possible.
+ | Extra compression is achieved in the case of GROMOS and coordinates of
+ | water molecules. GROMOS first writes out the Oxygen position, followed by
+ | the two hydrogens. In order to make the differences smaller (and thereby
+ | compression the data better) the order is changed into first one hydrogen
+ | then the oxygen, followed by the other hydrogen. This is rather special, but
+ | it shouldn't harm in the general case.
+ |
+ */
+
+int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
+{
+ int *ip = NULL;
+ int *buf = NULL;
+ gmx_bool bRead;
+
+ /* preallocate a small buffer and ip on the stack - if we need more
+ we can always malloc(). This is faster for small values of size: */
+ unsigned prealloc_size = 3*16;
+ int prealloc_ip[3*16], prealloc_buf[3*20];
+ int we_should_free = 0;
+
+ int minint[3], maxint[3], mindiff, *lip, diff;
+ int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
+ int minidx, maxidx;
+ unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
+ int flag, k;
+ int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
+ float *lfp, lf;
+ int tmp, *thiscoord, prevcoord[3];
+ unsigned int tmpcoord[30];
+
+ int bufsize, xdrid, lsize;
+ unsigned int bitsize;
+ float inv_precision;
+ int errval = 1;
+ int rc;
+
+ bRead = (xdrs->x_op == XDR_DECODE);
+ bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
+ prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
+
+ if (!bRead)
+ {
+ /* xdrs is open for writing */
+
+ if (xdr_int(xdrs, size) == 0)
+ {
+ return 0;
+ }
+ size3 = *size * 3;
+ /* when the number of coordinates is small, don't try to compress; just
+ * write them as floats using xdr_vector
+ */
+ if (*size <= 9)
+ {
+ return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
+ (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
+ }
+
+ if (xdr_float(xdrs, precision) == 0)
+ {
+ return 0;
+ }
+
+ if (size3 <= prealloc_size)
+ {
+ ip = prealloc_ip;
+ buf = prealloc_buf;
+ }
+ else
+ {
+ we_should_free = 1;
+ bufsize = size3 * 1.2;
+ ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
+ buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
+ if (ip == NULL || buf == NULL)
+ {
+ fprintf(stderr, "malloc failed\n");
+ exit(1);
+ }
+ }
+ /* buf[0-2] are special and do not contain actual data */
+ buf[0] = buf[1] = buf[2] = 0;
+ minint[0] = minint[1] = minint[2] = INT_MAX;
+ maxint[0] = maxint[1] = maxint[2] = INT_MIN;
+ prevrun = -1;
+ lfp = fp;
+ lip = ip;
+ mindiff = INT_MAX;
+ oldlint1 = oldlint2 = oldlint3 = 0;
+ while (lfp < fp + size3)
+ {
+ /* find nearest integer */
+ if (*lfp >= 0.0)
+ {
+ lf = *lfp * *precision + 0.5;
+ }
+ else
+ {
+ lf = *lfp * *precision - 0.5;
+ }
+ if (fabs(lf) > MAXABS)
+ {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint1 = lf;
+ if (lint1 < minint[0])
+ {
+ minint[0] = lint1;
+ }
+ if (lint1 > maxint[0])
+ {
+ maxint[0] = lint1;
+ }
+ *lip++ = lint1;
+ lfp++;
+ if (*lfp >= 0.0)
+ {
+ lf = *lfp * *precision + 0.5;
+ }
+ else
+ {
+ lf = *lfp * *precision - 0.5;
+ }
+ if (fabs(lf) > MAXABS)
+ {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint2 = lf;
+ if (lint2 < minint[1])
+ {
+ minint[1] = lint2;
+ }
+ if (lint2 > maxint[1])
+ {
+ maxint[1] = lint2;
+ }
+ *lip++ = lint2;
+ lfp++;
+ if (*lfp >= 0.0)
+ {
+ lf = *lfp * *precision + 0.5;
+ }
+ else
+ {
+ lf = *lfp * *precision - 0.5;
+ }
+ if (fabs(lf) > MAXABS)
+ {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint3 = lf;
+ if (lint3 < minint[2])
+ {
+ minint[2] = lint3;
+ }
+ if (lint3 > maxint[2])
+ {
+ maxint[2] = lint3;
+ }
+ *lip++ = lint3;
+ lfp++;
+ diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
+ if (diff < mindiff && lfp > fp + 3)
+ {
+ mindiff = diff;
+ }
+ oldlint1 = lint1;
+ oldlint2 = lint2;
+ oldlint3 = lint3;
+ }
+ if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
+ (xdr_int(xdrs, &(minint[1])) == 0) ||
+ (xdr_int(xdrs, &(minint[2])) == 0) ||
+ (xdr_int(xdrs, &(maxint[0])) == 0) ||
+ (xdr_int(xdrs, &(maxint[1])) == 0) ||
+ (xdr_int(xdrs, &(maxint[2])) == 0))
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+ if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
+ (float)maxint[1] - (float)minint[1] >= MAXABS ||
+ (float)maxint[2] - (float)minint[2] >= MAXABS)
+ {
+ /* turning value in unsigned by subtracting minint
+ * would cause overflow
+ */
+ errval = 0;
+ }
+ sizeint[0] = maxint[0] - minint[0]+1;
+ sizeint[1] = maxint[1] - minint[1]+1;
+ sizeint[2] = maxint[2] - minint[2]+1;
+
+ /* check if one of the sizes is to big to be multiplied */
+ if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
+ {
+ bitsizeint[0] = sizeofint(sizeint[0]);
+ bitsizeint[1] = sizeofint(sizeint[1]);
+ bitsizeint[2] = sizeofint(sizeint[2]);
+ bitsize = 0; /* flag the use of large sizes */
+ }
+ else
+ {
+ bitsize = sizeofints(3, sizeint);
+ }
+ lip = ip;
+ luip = (unsigned int *) ip;
+ smallidx = FIRSTIDX;
+ while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
+ {
+ smallidx++;
+ }
+ if (xdr_int(xdrs, &smallidx) == 0)
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+ maxidx = MIN(LASTIDX, smallidx + 8);
+ minidx = maxidx - 8; /* often this equal smallidx */
+ smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 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)
+ {
+ sendbits(buf, bitsizeint[0], tmpcoord[0]);
+ sendbits(buf, bitsizeint[1], tmpcoord[1]);
+ sendbits(buf, bitsizeint[2], tmpcoord[2]);
+ }
+ else
+ {
+ sendints(buf, 3, bitsize, sizeint, tmpcoord);
+ }
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+ thiscoord = thiscoord + 3;
+ i++;
+
+ run = 0;
+ if (is_small == 0 && is_smaller == -1)
+ {
+ is_smaller = 0;
+ }
+ while (is_small && run < 8*3)
+ {
+ if (is_smaller == -1 && (
+ SQR(thiscoord[0] - prevcoord[0]) +
+ SQR(thiscoord[1] - prevcoord[1]) +
+ SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller))
+ {
+ is_smaller = 0;
+ }
+
+ tmpcoord[run++] = thiscoord[0] - prevcoord[0] + 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;
+ sendbits(buf, 1, 1); /* flag the change in run-length */
+ sendbits(buf, 5, run+is_smaller+1);
+ }
+ else
+ {
+ sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
+ }
+ for (k = 0; k < run; k += 3)
+ {
+ sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
+ }
+ if (is_smaller != 0)
+ {
+ smallidx += is_smaller;
+ if (is_smaller < 0)
+ {
+ smallnum = smaller;
+ smaller = magicints[smallidx-1] / 2;
+ }
+ else
+ {
+ smaller = smallnum;
+ smallnum = magicints[smallidx] / 2;
+ }
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
+ }
+ }
+ if (buf[1] != 0)
+ {
+ buf[0]++;
+ }
+ /* buf[0] holds the length in bytes */
+ if (xdr_int(xdrs, &(buf[0])) == 0)
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+
+ rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return rc;
+
+ }
+ else
+ {
+
+ /* xdrs is open for reading */
+
+ if (xdr_int(xdrs, &lsize) == 0)
+ {
+ return 0;
+ }
+ if (*size != 0 && lsize != *size)
+ {
+ fprintf(stderr, "wrong number of coordinates in xdr3dfcoord; "
+ "%d arg vs %d in file", *size, lsize);
+ }
+ *size = lsize;
+ size3 = *size * 3;
+ if (*size <= 9)
+ {
+ *precision = -1;
+ return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
+ (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
+ }
+ if (xdr_float(xdrs, precision) == 0)
+ {
+ return 0;
+ }
+
+ if (size3 <= prealloc_size)
+ {
+ ip = prealloc_ip;
+ buf = prealloc_buf;
+ }
+ else
+ {
+ we_should_free = 1;
+ bufsize = size3 * 1.2;
+ ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
+ buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
+ if (ip == NULL || buf == NULL)
+ {
+ fprintf(stderr, "malloc failed\n");
+ exit(1);
+ }
+ }
+
+ buf[0] = buf[1] = buf[2] = 0;
+
+ if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
+ (xdr_int(xdrs, &(minint[1])) == 0) ||
+ (xdr_int(xdrs, &(minint[2])) == 0) ||
+ (xdr_int(xdrs, &(maxint[0])) == 0) ||
+ (xdr_int(xdrs, &(maxint[1])) == 0) ||
+ (xdr_int(xdrs, &(maxint[2])) == 0))
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 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);
+ }
+
+ if (xdr_int(xdrs, &smallidx) == 0)
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+ maxidx = MIN(LASTIDX, smallidx + 8);
+ minidx = maxidx - 8; /* often this equal smallidx */
+ smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
+ smallnum = magicints[smallidx] / 2;
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
+ larger = magicints[maxidx];
+
+ /* buf[0] holds the length in bytes */
+
+ if (xdr_int(xdrs, &(buf[0])) == 0)
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+
+ if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
+ {
+ if (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 0;
+ }
+
+
+
+ buf[0] = buf[1] = buf[2] = 0;
+
+ lfp = fp;
+ inv_precision = 1.0 / *precision;
+ run = 0;
+ i = 0;
+ lip = ip;
+ while (i < lsize)
+ {
+ thiscoord = (int *)(lip) + i * 3;
+
+ if (bitsize == 0)
+ {
+ thiscoord[0] = receivebits(buf, bitsizeint[0]);
+ thiscoord[1] = receivebits(buf, bitsizeint[1]);
+ thiscoord[2] = receivebits(buf, bitsizeint[2]);
+ }
+ else
+ {
+ receiveints(buf, 3, bitsize, sizeint, thiscoord);
+ }
+
+ i++;
+ thiscoord[0] += minint[0];
+ thiscoord[1] += minint[1];
+ thiscoord[2] += minint[2];
+
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+
+
+ flag = receivebits(buf, 1);
+ is_smaller = 0;
+ if (flag == 1)
+ {
+ run = receivebits(buf, 5);
+ is_smaller = run % 3;
+ run -= is_smaller;
+ is_smaller--;
+ }
+ if (run > 0)
+ {
+ thiscoord += 3;
+ for (k = 0; k < run; k += 3)
+ {
+ receiveints(buf, 3, smallidx, sizesmall, thiscoord);
+ i++;
+ thiscoord[0] += prevcoord[0] - 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 (we_should_free)
+ {
+ free(ip);
+ free(buf);
+ }
+ return 1;
+}
+
+
+
+/******************************************************************
+
+ XTC files have a relatively simple structure.
+ They have a header of 16 bytes and the rest are
+ the compressed coordinates of the files. Due to the
+ compression 00 is not present in the coordinates.
+ The first 4 bytes of the header are the magic number
+ 1995 (0x000007CB). If we find this number we are guaranteed
+ to be in the header, due to the presence of so many zeros.
+ The second 4 bytes are the number of atoms in the frame, and is
+ assumed to be constant. The third 4 bytes are the frame number.
+ The last 4 bytes are a floating point representation of the time.
+
+ ********************************************************************/
+
+/* Must match definition in xtcio.c */
+#ifndef XTC_MAGIC
+#define XTC_MAGIC 1995
+#endif
+
+static const int header_size = 16;
+
+/* Check if we are at the header start.
+ At the same time it will also read 1 int */
+static int xtc_at_header_start(FILE *fp, XDR *xdrs,
+ int natoms, int * timestep, float * time)
+{
+ int i_inp[3];
+ float f_inp[10];
+ int i;
+ gmx_off_t off;
+
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+ /* read magic natoms and timestep */
+ for (i = 0; i < 3; i++)
+ {
+ if (!xdr_int(xdrs, &(i_inp[i])))
+ {
+ gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
+ return -1;
+ }
+ }
+ /* quick return */
+ if (i_inp[0] != XTC_MAGIC)
+ {
+ if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
+ {
+ return -1;
+ }
+ return 0;
+ }
+ /* read time and box */
+ for (i = 0; i < 10; i++)
+ {
+ if (!xdr_float(xdrs, &(f_inp[i])))
+ {
+ gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET);
+ return -1;
+ }
+ }
+ /* Make a rigourous check to see if we are in the beggining of a header
+ Hopefully there are no ambiguous cases */
+ /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
+ right triangle and that the first element must be nonzero unless the entire matrix is zero
+ */
+ if (i_inp[1] == natoms &&
+ ((f_inp[1] != 0 && f_inp[6] == 0) ||
+ (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
+ {
+ if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
+ {
+ return -1;
+ }
+ *time = f_inp[0];
+ *timestep = i_inp[2];
+ return 1;
+ }
+ if (gmx_fseek(fp, off+XDR_INT_SIZE, SEEK_SET))
+ {
+ return -1;
+ }
+ return 0;
+}
+
+static int
+xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
+{
+ gmx_off_t off;
+ int step;
+ float time;
+ int ret;
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+
+ /* read one int just to make sure we dont read this frame but the next */
+ xdr_int(xdrs, &step);
+ while (1)
+ {
+ ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
+ if (ret == 1)
+ {
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ return -1;
+ }
+ return step;
+ }
+ else if (ret == -1)
+ {
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ return -1;
+ }
+ }
+ }
+ return -1;
+}
+
+
+static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
+ gmx_bool * bOK)
+{
+ gmx_off_t off;
+ float time;
+ int step;
+ int ret;
+ *bOK = 0;
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+ /* read one int just to make sure we dont read this frame but the next */
+ xdr_int(xdrs, &step);
+ while (1)
+ {
+ ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
+ if (ret == 1)
+ {
+ *bOK = 1;
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ *bOK = 0;
+ return -1;
+ }
+ return time;
+ }
+ else if (ret == -1)
+ {
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ return -1;
+ }
+ return -1;
+ }
+ }
+ return -1;
+}
+
+
+static float
+xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
+{
+ gmx_off_t off;
+ int step;
+ float time;
+ int ret;
+ *bOK = 0;
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+
+ while (1)
+ {
+ ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
+ if (ret == 1)
+ {
+ *bOK = 1;
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ *bOK = 0;
+ return -1;
+ }
+ return time;
+ }
+ else if (ret == -1)
+ {
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ return -1;
+ }
+ return -1;
+ }
+ else if (ret == 0)
+ {
+ /*Go back.*/
+ if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
+ {
+ return -1;
+ }
+ }
+ }
+ return -1;
+}
+
+/* Currently not used, just for completeness */
+static int
+xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
+{
+ gmx_off_t off;
+ int ret;
+ int step;
+ float time;
+ *bOK = 0;
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+
+
+ while (1)
+ {
+ ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
+ if (ret == 1)
+ {
+ *bOK = 1;
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ *bOK = 0;
+ return -1;
+ }
+ return step;
+ }
+ else if (ret == -1)
+ {
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ return -1;
+ }
+ return -1;
+
+ }
+ else if (ret == 0)
+ {
+ /*Go back.*/
+ if (gmx_fseek(fp, -2*XDR_INT_SIZE, SEEK_CUR))
+ {
+ return -1;
+ }
+ }
+ }
+ return -1;
+}
+
+
+
+static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
+{
+ int inp;
+ gmx_off_t res;
+ int ret;
+ int step;
+ float time;
+ /* read one int just to make sure we dont read this frame but the next */
+ xdr_int(xdrs, &step);
+ while (1)
+ {
+ ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
+ if (ret == 1)
+ {
+ if ((res = gmx_ftell(fp)) >= 0)
+ {
+ return res - XDR_INT_SIZE;
+ }
+ else
+ {
+ return res;
+ }
+ }
+ else if (ret == -1)
+ {
+ return -1;
+ }
+ }
+ return -1;
+}
+
+
+static
+float
+xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
+{
+ float res;
+ float tinit;
+ gmx_off_t off;
+
+ *bOK = 0;
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+
+ tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
+
+ if (!(*bOK))
+ {
+ return -1;
+ }
+
+ res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
+
+ if (!(*bOK))
+ {
+ return -1;
+ }
+
+ res -= tinit;
+ if (0 != gmx_fseek(fp, off, SEEK_SET))
+ {
+ *bOK = 0;
+ return -1;
+ }
+ return res;
+}
+
- if (fr != frame && abs(low-high) > header_size)
+int
+xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
+{
+ gmx_off_t low = 0;
+ gmx_off_t high, pos;
+
+
+ /* round to 4 bytes */
+ int fr;
+ gmx_off_t offset;
+ if (gmx_fseek(fp, 0, SEEK_END))
+ {
+ return -1;
+ }
+
+ if ((high = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+
+ /* round to 4 bytes */
+ high /= XDR_INT_SIZE;
+ high *= XDR_INT_SIZE;
+ offset = ((high/2)/XDR_INT_SIZE)*XDR_INT_SIZE;
+
+ if (gmx_fseek(fp, offset, SEEK_SET))
+ {
+ return -1;
+ }
+
+ while (1)
+ {
+ fr = xtc_get_next_frame_number(fp, xdrs, natoms);
+ if (fr < 0)
+ {
+ return -1;
+ }
- if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
- - time) >= dt && dt_sign >= 0)
- || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
- > header_size))
++ if (fr != frame && llabs(low-high) > header_size)
+ {
+ if (fr < frame)
+ {
+ low = offset;
+ }
+ else
+ {
+ high = offset;
+ }
+ /* round to 4 bytes */
+ offset = (((high+low)/2)/4)*4;
+
+ if (gmx_fseek(fp, offset, SEEK_SET))
+ {
+ return -1;
+ }
+ }
+ else
+ {
+ break;
+ }
+ }
+ if (offset <= header_size)
+ {
+ offset = low;
+ }
+
+ if (gmx_fseek(fp, offset, SEEK_SET))
+ {
+ return -1;
+ }
+
+ if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
+ {
+ /* we probably hit an end of file */
+ return -1;
+ }
+
+ if (gmx_fseek(fp, pos, SEEK_SET))
+ {
+ return -1;
+ }
+
+ return 0;
+}
+
+
+
+int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
+{
+ float t;
+ float dt;
+ gmx_bool bOK = FALSE;
+ gmx_off_t low = 0;
+ gmx_off_t high, offset, pos;
+ int res;
+ int dt_sign = 0;
+
+ if (bSeekForwardOnly)
+ {
+ low = gmx_ftell(fp);
+ }
+ if (gmx_fseek(fp, 0, SEEK_END))
+ {
+ return -1;
+ }
+
+ if ((high = gmx_ftell(fp)) < 0)
+ {
+ return -1;
+ }
+ /* round to int */
+ high /= XDR_INT_SIZE;
+ high *= XDR_INT_SIZE;
+ offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
+
+ if (gmx_fseek(fp, offset, SEEK_SET))
+ {
+ return -1;
+ }
+
+
+ /*
+ * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
+ dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
+
+ if (!bOK)
+ {
+ return -1;
+ }
+ */
+
+ while (1)
+ {
+ dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
+ if (!bOK)
+ {
+ return -1;
+ }
+ else
+ {
+ if (dt > 0)
+ {
+ if (dt_sign == -1)
+ {
+ /* Found a place in the trajectory that has positive time step while
+ other has negative time step */
+ return -2;
+ }
+ dt_sign = 1;
+ }
+ else if (dt < 0)
+ {
+ if (dt_sign == 1)
+ {
+ /* Found a place in the trajectory that has positive time step while
+ other has negative time step */
+ return -2;
+ }
+ dt_sign = -1;
+ }
+ }
+ t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
+ if (!bOK)
+ {
+ return -1;
+ }
+
+ /* If we are before the target time and the time step is positive or 0, or we have
+ after the target time and the time step is negative, or the difference between
+ the current time and the target time is bigger than dt and above all the distance between high
+ and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
+ if we reached the solution */
- if (abs(low - high) <= header_size)
++ if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) ||
++ ((t - time) >= dt && dt_sign >= 0) || ((time - t) >= -dt && dt_sign < 0)) &&
++ (llabs(low - high) > header_size))
+ {
+ if (dt >= 0 && dt_sign != -1)
+ {
+ if (t < time)
+ {
+ low = offset;
+ }
+ else
+ {
+ high = offset;
+ }
+ }
+ else if (dt <= 0 && dt_sign == -1)
+ {
+ if (t >= time)
+ {
+ low = offset;
+ }
+ else
+ {
+ high = offset;
+ }
+ }
+ else
+ {
+ /* We should never reach here */
+ return -1;
+ }
+ /* round to 4 bytes and subtract header*/
+ offset = (((high + low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
+ if (gmx_fseek(fp, offset, SEEK_SET))
+ {
+ return -1;
+ }
+ }
+ else
+ {
++ if (llabs(low - high) <= header_size)
+ {
+ break;
+ }
+ /* re-estimate dt */
+ if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
+ {
+ if (bOK)
+ {
+ dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
+ }
+ }
+ if (t >= time && t - time < dt)
+ {
+ break;
+ }
+ }
+ }
+
+ if (offset <= header_size)
+ {
+ offset = low;
+ }
+
+ gmx_fseek(fp, offset, SEEK_SET);
+
+ if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
+ {
+ return -1;
+ }
+
+ if (gmx_fseek(fp, pos, SEEK_SET))
+ {
+ return -1;
+ }
+ return 0;
+}
+
+float
+xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
+{
+ float time;
+ gmx_off_t off;
+ int res;
+ *bOK = 1;
+ off = gmx_ftell(fp);
+ if (off < 0)
+ {
+ *bOK = 0;
+ return -1;
+ }
+
+ if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END)) != 0)
+ {
+ *bOK = 0;
+ return -1;
+ }
+
+ time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
+ if (!(*bOK))
+ {
+ return -1;
+ }
+
+ if ( (res = gmx_fseek(fp, off, SEEK_SET)) != 0)
+ {
+ *bOK = 0;
+ return -1;
+ }
+ return time;
+}
+
+
+int
+xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
+{
+ int frame;
+ gmx_off_t off;
+ int res;
+ *bOK = 1;
+
+ if ((off = gmx_ftell(fp)) < 0)
+ {
+ *bOK = 0;
+ return -1;
+ }
+
+
+ if (gmx_fseek(fp, -3*XDR_INT_SIZE, SEEK_END))
+ {
+ *bOK = 0;
+ return -1;
+ }
+
+ frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
+ if (!bOK)
+ {
+ return -1;
+ }
+
+
+ if (gmx_fseek(fp, off, SEEK_SET))
+ {
+ *bOK = 0;
+ return -1;
+ }
+
+ return frame;
+}
--- /dev/null
- if (abs(ct[0]-1.0) > 1e-6)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include "typedefs.h"
+#include "gromacs/utility/smalloc.h"
+#include "vec.h"
+#include "geminate.h"
+
+#include "gromacs/legacyheaders/gmx_fatal.h"
+#include "gromacs/utility/gmxomp.h"
+
+static void missing_code_message()
+{
+ fprintf(stderr, "You have requested code to run that is deprecated.\n");
+ fprintf(stderr, "Revert to an older GROMACS version or help in porting the code.\n");
+}
+
+/* The first few sections of this file contain functions that were adopted,
+ * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
+ * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
+ * This is also the case with the function eq10v2().
+ *
+ * The parts menetioned in the previous paragraph were contributed under the BSD license.
+ */
+
+
+/* This first part is from complex.c which I recieved from Omer Markowitch.
+ * - Erik Marklund
+ *
+ * ------------- from complex.c ------------- */
+
+/* Complexation of a paired number (r,i) */
+static gem_complex gem_cmplx(double r, double i)
+{
+ gem_complex value;
+ value.r = r;
+ value.i = i;
+ return value;
+}
+
+/* Complexation of a real number, x */
+static gem_complex gem_c(double x)
+{
+ gem_complex value;
+ value.r = x;
+ value.i = 0;
+ return value;
+}
+
+/* Magnitude of a complex number z */
+static double gem_cx_abs(gem_complex z)
+{
+ return (sqrt(z.r*z.r+z.i*z.i));
+}
+
+/* Addition of two complex numbers z1 and z2 */
+static gem_complex gem_cxadd(gem_complex z1, gem_complex z2)
+{
+ gem_complex value;
+ value.r = z1.r+z2.r;
+ value.i = z1.i+z2.i;
+ return value;
+}
+
+/* Addition of a complex number z1 and a real number r */
+static gem_complex gem_cxradd(gem_complex z, double r)
+{
+ gem_complex value;
+ value.r = z.r + r;
+ value.i = z.i;
+ return value;
+}
+
+/* Subtraction of two complex numbers z1 and z2 */
+static gem_complex gem_cxsub(gem_complex z1, gem_complex z2)
+{
+ gem_complex value;
+ value.r = z1.r-z2.r;
+ value.i = z1.i-z2.i;
+ return value;
+}
+
+/* Multiplication of two complex numbers z1 and z2 */
+static gem_complex gem_cxmul(gem_complex z1, gem_complex z2)
+{
+ gem_complex value;
+ value.r = z1.r*z2.r-z1.i*z2.i;
+ value.i = z1.r*z2.i+z1.i*z2.r;
+ return value;
+}
+
+/* Square of a complex number z */
+static gem_complex gem_cxsq(gem_complex z)
+{
+ gem_complex value;
+ value.r = z.r*z.r-z.i*z.i;
+ value.i = z.r*z.i*2.;
+ return value;
+}
+
+/* multiplication of a complex number z and a real number r */
+static gem_complex gem_cxrmul(gem_complex z, double r)
+{
+ gem_complex value;
+ value.r = z.r*r;
+ value.i = z.i*r;
+ return value;
+}
+
+/* Division of two complex numbers z1 and z2 */
+static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2)
+{
+ gem_complex value;
+ double num;
+ num = z2.r*z2.r+z2.i*z2.i;
+ if (num == 0.)
+ {
+ fprintf(stderr, "ERROR in gem_cxdiv function\n");
+ }
+ value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num;
+ return value;
+}
+
+/* division of a complex z number by a real number x */
+static gem_complex gem_cxrdiv(gem_complex z, double r)
+{
+ gem_complex value;
+ value.r = z.r/r;
+ value.i = z.i/r;
+ return value;
+}
+
+/* division of a real number r by a complex number x */
+static gem_complex gem_rxcdiv(double r, gem_complex z)
+{
+ gem_complex value;
+ double f;
+ f = r/(z.r*z.r+z.i*z.i);
+ value.r = f*z.r;
+ value.i = -f*z.i;
+ return value;
+}
+
+/* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
+static gem_complex gem_cxdexp(gem_complex z)
+{
+ gem_complex value;
+ double exp_z_r;
+ exp_z_r = exp(z.r);
+ value.r = exp_z_r*cos(z.i);
+ value.i = exp_z_r*sin(z.i);
+ return value;
+}
+
+/* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z), */
+/* where -PI < Arg(z) < PI */
+static gem_complex gem_cxlog(gem_complex z)
+{
+ gem_complex value;
+ double mag2;
+ mag2 = z.r*z.r+z.i*z.i;
+ if (mag2 < 0.)
+ {
+ fprintf(stderr, "ERROR in gem_cxlog func\n");
+ }
+ value.r = log(sqrt(mag2));
+ if (z.r == 0.)
+ {
+ value.i = PI/2.;
+ if (z.i < 0.)
+ {
+ value.i = -value.i;
+ }
+ }
+ else
+ {
+ value.i = atan2(z.i, z.r);
+ }
+ return value;
+}
+
+/* Square root of a complex number z = |z| exp(I*the) -- z^(1/2) */
+/* z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */
+/* where 0 < the < 2*PI */
+static gem_complex gem_cxdsqrt(gem_complex z)
+{
+ gem_complex value;
+ double sq;
+ sq = gem_cx_abs(z);
+ value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */
+ value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */
+ if (z.i < 0.)
+ {
+ value.r = -value.r;
+ }
+ return value;
+}
+
+/* Complex power of a complex number z1^z2 */
+static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2)
+{
+ gem_complex value;
+ value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2));
+ return value;
+}
+
+/* ------------ end of complex.c ------------ */
+
+/* This next part was derived from cubic.c, also received from Omer Markovitch.
+ * ------------- from cubic.c ------------- */
+
+/* Solver for a cubic equation: x^3-a*x^2+b*x-c=0 */
+static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam,
+ double a, double b, double c)
+{
+ double t1, t2, two_3, temp;
+ gem_complex ctemp, ct3;
+
+ two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c;
+ temp = 4.*t1*t1*t1+t2*t2;
+
+ ctemp = gem_cmplx(temp, 0.); ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp));
+ ct3 = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.));
+
+ ctemp = gem_rxcdiv(-two_3*t1/3., ct3);
+ ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3));
+
+ *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp);
+
+ ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.))));
+ ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam));
+ ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam);
+ *al = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
+ *be = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
+}
+
+/* ------------ end of cubic.c ------------ */
+
+/* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
+ * ------------- from [cr]error.c ------------- */
+
+/************************************************************/
+/* Real valued error function and related functions */
+/* x, y : real variables */
+/* erf(x) : error function */
+/* erfc(x) : complementary error function */
+/* omega(x) : exp(x*x)*erfc(x) */
+/* W(x,y) : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
+/************************************************************/
+
+/*---------------------------------------------------------------------------*/
+/* Utilzed the series approximation of erf(x) */
+/* Relative error=|err(x)|/erf(x)<EPS */
+/* Handbook of Mathematical functions, Abramowitz, p 297 */
+/* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
+/*---------------------------------------------------------------------------*/
+/* This gives erfc(z) correctly only upto >10-15 */
+
+static double gem_erf(double x)
+{
+ double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12;
+ temp = x;
+ sum = temp;
+ xm = 26.;
+ x2 = x*x;
+ x4 = x2*x2;
+ x6 = x4*x2;
+ x8 = x6*x2;
+ x10 = x8*x2;
+ x12 = x10*x2;
+ exp2 = exp(-x2);
+ if (x <= xm)
+ {
+ for (n = 1.; n <= 2000.; n += 1.)
+ {
+ temp *= 2.*x2/(2.*n+1.);
+ sum += temp;
+ if (fabs(temp/sum) < 1.E-16)
+ {
+ break;
+ }
+ }
+
+ if (n >= 2000.)
+ {
+ fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n);
+ }
+ sum *= 2./sPI*exp2;
+ }
+ else
+ {
+ /* from the asymptotic expansion of experfc(x) */
+ sum = (1. - 0.5/x2 + 0.75/x4
+ - 1.875/x6 + 6.5625/x8
+ - 29.53125/x10 + 162.421875/x12)
+ / sPI/x;
+ sum *= exp2; /* now sum is erfc(x) */
+ sum = -sum+1.;
+ }
+ return x >= 0.0 ? sum : -sum;
+}
+
+/* Result --> Alex's code for erfc and experfc behaves better than this */
+/* complementray error function. Returns 1.-erf(x) */
+static double gem_erfc(double x)
+{
+ double t, z, ans;
+ z = fabs(x);
+ t = 1.0/(1.0+0.5*z);
+
+ ans = t * exp(-z*z - 1.26551223 +
+ t*(1.00002368 +
+ t*(0.37409196 +
+ t*(0.09678418 +
+ t*(-0.18628806 +
+ t*(0.27886807 +
+ t*(-1.13520398 +
+ t*(1.48851587 +
+ t*(-0.82215223 +
+ t*0.17087277)))))))));
+
+ return x >= 0.0 ? ans : 2.0-ans;
+}
+
+/* omega(x)=exp(x^2)erfc(x) */
+static double gem_omega(double x)
+{
+ double xm, ans, xx, x4, x6, x8, x10, x12;
+ xm = 26;
+ xx = x*x;
+ x4 = xx*xx;
+ x6 = x4*xx;
+ x8 = x6*xx;
+ x10 = x8*xx;
+ x12 = x10*xx;
+
+ if (x <= xm)
+ {
+ ans = exp(xx)*gem_erfc(x);
+ }
+ else
+ {
+ /* Asymptotic expansion */
+ ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x;
+ }
+ return ans;
+}
+
+/*---------------------------------------------------------------------------*/
+/* Utilzed the series approximation of erf(z=x+iy) */
+/* Relative error=|err(z)|/|erf(z)|<EPS */
+/* Handbook of Mathematical functions, Abramowitz, p 299 */
+/* comega(z=x+iy)=cexp(z^2)*cerfc(z) */
+/*---------------------------------------------------------------------------*/
+static gem_complex gem_comega(gem_complex z)
+{
+ gem_complex value;
+ double x, y;
+ double sumr, sumi, n, n2, f, temp, temp1;
+ double x2, y2, cos_2xy, sin_2xy, cosh_2xy, sinh_2xy, cosh_ny, sinh_ny, exp_y2;
+
+ x = z.r;
+ y = z.i;
+ x2 = x*x;
+ y2 = y*y;
+ sumr = 0.;
+ sumi = 0.;
+ cos_2xy = cos(2.*x*y);
+ sin_2xy = sin(2.*x*y);
+ cosh_2xy = cosh(2.*x*y);
+ sinh_2xy = sinh(2.*x*y);
+ exp_y2 = exp(-y2);
+
+ for (n = 1.0, temp = 0.; n <= 2000.; n += 1.0)
+ {
+ n2 = n*n;
+ cosh_ny = cosh(n*y);
+ sinh_ny = sinh(n*y);
+ f = exp(-n2/4.)/(n2+4.*x2);
+ /* if(f<1.E-200) break; */
+ sumr += (2.*x - 2.*x*cosh_ny*cos_2xy + n*sinh_ny*sin_2xy)*f;
+ sumi += (2.*x*cosh_ny*sin_2xy + n*sinh_ny*cos_2xy)*f;
+ temp1 = sqrt(sumr*sumr+sumi*sumi);
+ if (fabs((temp1-temp)/temp1) < 1.E-16)
+ {
+ break;
+ }
+ temp = temp1;
+ }
+ if (n == 2000.)
+ {
+ fprintf(stderr, "iteration exceeds %lg\n", n);
+ }
+ sumr *= 2./PI;
+ sumi *= 2./PI;
+
+ if (x != 0.)
+ {
+ f = 1./2./PI/x;
+ }
+ else
+ {
+ f = 0.;
+ }
+ value.r = gem_omega(x)-(f*(1.-cos_2xy)+sumr);
+ value.i = -(f*sin_2xy+sumi);
+ value = gem_cxmul(value, gem_cmplx(exp_y2*cos_2xy, exp_y2*sin_2xy));
+ return (value);
+}
+
+/* ------------ end of [cr]error.c ------------ */
+
+/*_ REVERSIBLE GEMINATE RECOMBINATION
+ *
+ * Here are the functions for reversible geminate recombination. */
+
+/* Changes the unit from square cm per s to square Ångström per ps,
+ * since Omers code uses the latter units while g_mds outputs the former.
+ * g_hbond expects a diffusion coefficent given in square cm per s. */
+static double sqcm_per_s_to_sqA_per_ps (real D)
+{
+ fprintf(stdout, "Diffusion coefficient is %f A^2/ps\n", D*1e4);
+ return (double)(D*1e4);
+}
+
+
+static double eq10v2(double theoryCt[], double time[], int manytimes,
+ double ka, double kd, t_gemParams *params)
+{
+ /* Finding the 3 roots */
+ int i;
+ double kD, D, r, a, b, c, tsqrt, sumimaginary;
+ gem_complex
+ alpha, beta, gamma,
+ c1, c2, c3, c4,
+ oma, omb, omc,
+ part1, part2, part3, part4;
+
+ kD = params->kD;
+ D = params->D;
+ r = params->sigma;
+ a = (1.0 + ka/kD) * sqrt(D)/r;
+ b = kd;
+ c = kd * sqrt(D)/r;
+
+ gem_solve(&alpha, &beta, &gamma, a, b, c);
+ /* Finding the 3 roots */
+
+ sumimaginary = 0;
+ part1 = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta, gamma), gem_cxsub(beta, gamma))); /* 1(2+3)(2-3) */
+ part2 = gem_cxmul(beta, gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha))); /* 2(3+1)(3-1) */
+ part3 = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta))); /* 3(1+2)(1-2) */
+ part4 = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */
+
+#pragma omp parallel for \
+ private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
+ reduction(+:sumimaginary) \
+ default(shared) \
+ schedule(guided)
+ for (i = 0; i < manytimes; i++)
+ {
+ tsqrt = sqrt(time[i]);
+ oma = gem_comega(gem_cxrmul(alpha, tsqrt));
+ c1 = gem_cxmul(oma, gem_cxdiv(part1, part4));
+ omb = gem_comega(gem_cxrmul(beta, tsqrt));
+ c2 = gem_cxmul(omb, gem_cxdiv(part2, part4));
+ omc = gem_comega(gem_cxrmul(gamma, tsqrt));
+ c3 = gem_cxmul(omc, gem_cxdiv(part3, part4));
+ c4.r = c1.r+c2.r+c3.r;
+ c4.i = c1.i+c2.i+c3.i;
+ theoryCt[i] = c4.r;
+ sumimaginary += c4.i * c4.i;
+ }
+
+ return sumimaginary;
+
+} /* eq10v2 */
+
+/* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
+static double getLogIndex(const int i, const t_gemParams *params)
+{
+ return (exp(((double)(i)) * params->logQuota) -1);
+}
+
+extern t_gemParams *init_gemParams(const double sigma, const double D,
+ const real *t, const int len, const int nFitPoints,
+ const real begFit, const real endFit,
+ const real ballistic, const int nBalExp)
+{
+ double tDelta;
+ t_gemParams *p;
+
+ snew(p, 1);
+
+ /* A few hardcoded things here. For now anyway. */
+/* p->ka_min = 0; */
+/* p->ka_max = 100; */
+/* p->dka = 10; */
+/* p->kd_min = 0; */
+/* p->kd_max = 2; */
+/* p->dkd = 0.2; */
+ p->ka = 0;
+ p->kd = 0;
+/* p->lsq = -1; */
+/* p->lifetime = 0; */
+ p->sigma = sigma*10; /* Omer uses Å, not nm */
+/* p->lsq_old = 99999; */
+ p->D = sqcm_per_s_to_sqA_per_ps(D);
+ p->kD = 4*acos(-1.0)*sigma*p->D;
+
+
+ /* Parameters used by calcsquare(). Better to calculate them
+ * here than in calcsquare every time it's called. */
+ p->len = len;
+/* p->logAfterTime = logAfterTime; */
+ tDelta = (t[len-1]-t[0]) / len;
+ if (tDelta <= 0)
+ {
+ gmx_fatal(FARGS, "Time between frames is non-positive!");
+ }
+ else
+ {
+ p->tDelta = tDelta;
+ }
+
+ p->nExpFit = nBalExp;
+/* p->nLin = logAfterTime / tDelta; */
+ p->nFitPoints = nFitPoints;
+ p->begFit = begFit;
+ p->endFit = endFit;
+ p->logQuota = (double)(log(p->len))/(p->nFitPoints-1);
+/* if (p->nLin <= 0) { */
+/* fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
+/* sfree(p); */
+/* return NULL; */
+/* } */
+/* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
+/* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
+/* p->logPF = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
+/* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */
+
+/* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
+ p->ballistic = ballistic;
+ return p;
+}
+
+/* There was a misunderstanding regarding the fitting. From our
+ * recent correspondence it appears that Omer's code require
+ * the ACF data on a log-scale and does not operate on the raw data.
+ * This needs to be redone in gemFunc_residual() as well as in the
+ * t_gemParams structure. */
+
+static real* d2r(const double *d, const int nn);
+
+extern real fitGemRecomb(double gmx_unused *ct,
+ double gmx_unused *time,
+ double gmx_unused **ctFit,
+ const int gmx_unused nData,
+ t_gemParams gmx_unused *params)
+{
+
+ int nThreads, i, iter, status, maxiter;
+ real size, d2, tol, *dumpdata;
+ size_t p, n;
+ gemFitData *GD;
+ char *dumpstr, dumpname[128];
+
+ missing_code_message();
+ return -1;
+
+}
+
+
+/* Removes the ballistic term from the beginning of the ACF,
+ * just like in Omer's paper.
+ */
+extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative)
+{
+
+ /* Fit with 4 exponentials and one constant term,
+ * subtract the fatest exponential. */
+
+ int nData, i, status, iter;
+ balData *BD;
+ double *guess, /* Initial guess. */
+ *A, /* The fitted parameters. (A1, B1, A2, B2,... C) */
+ a[2],
+ ddt[2];
+ gmx_bool sorted;
+ size_t n;
+ size_t p;
+
+ nData = 0;
+ do
+ {
+ nData++;
+ }
+ while (t[nData] < tMax+t[0] && nData < len);
+
+ p = nexp*2+1; /* Number of parameters. */
+
+ missing_code_message();
+ return;
+}
+
+extern void dumpN(const real *e, const int nn, char *fn)
+{
+ /* For debugging only */
+ int i;
+ FILE *f;
+ char standardName[] = "Nt.xvg";
+ if (fn == NULL)
+ {
+ fn = standardName;
+ }
+
+ f = fopen(fn, "w");
+ fprintf(f,
+ "@ type XY\n"
+ "@ xaxis label \"Frame\"\n"
+ "@ yaxis label \"N\"\n"
+ "@ s0 line type 3\n");
+
+ for (i = 0; i < nn; i++)
+ {
+ fprintf(f, "%-10i %-g\n", i, e[i]);
+ }
+
+ fclose(f);
+}
+
+static real* d2r(const double *d, const int nn)
+{
+ real *r;
+ int i;
+
+ snew(r, nn);
+ for (i = 0; i < nn; i++)
+ {
+ r[i] = (real)d[i];
+ }
+
+ return r;
+}
+
+static void _patchBad(double *ct, int n, double dy)
+{
+ /* Just do lin. interpolation for now. */
+ int i;
+
+ for (i = 1; i < n; i++)
+ {
+ ct[i] = ct[0]+i*dy;
+ }
+}
+
+static void patchBadPart(double *ct, int n)
+{
+ _patchBad(ct, n, (ct[n] - ct[0])/n);
+}
+
+static void patchBadTail(double *ct, int n)
+{
+ _patchBad(ct+1, n-1, ct[1]-ct[0]);
+
+}
+
+extern void fixGemACF(double *ct, int len)
+{
+ int i, j, b, e;
+ gmx_bool bBad;
+
+ /* Let's separate two things:
+ * - identification of bad parts
+ * - patching of bad parts.
+ */
+
+ b = 0; /* Start of a bad stretch */
+ e = 0; /* End of a bad stretch */
+ bBad = FALSE;
+
+ /* An acf of binary data must be one at t=0. */
- fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0));
++ if (fabs(ct[0]-1.0) > 1e-6)
+ {
+ ct[0] = 1.0;
++ fprintf(stderr, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct[0]-1.0));
+ }
+
+ for (i = 0; i < len; i++)
+ {
+
+#ifdef HAS_ISFINITE
+ if (isfinite(ct[i]))
+#elif defined(HAS__ISFINITE)
+ if (_isfinite(ct[i]))
+#else
+ if (1)
+#endif
+ {
+ if (!bBad)
+ {
+ /* Still on a good stretch. Proceed.*/
+ continue;
+ }
+
+ /* Patch up preceding bad stretch. */
+ if (i == (len-1))
+ {
+ /* It's the tail */
+ if (b <= 1)
+ {
+ gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting.");
+ }
+ patchBadTail(&(ct[b-2]), (len-b)+1);
+ }
+
+ e = i;
+ patchBadPart(&(ct[b-1]), (e-b)+1);
+ bBad = FALSE;
+ }
+ else
+ {
+ if (!bBad)
+ {
+ b = i;
+
+ bBad = TRUE;
+ }
+ }
+ }
+}
--- /dev/null
- if (bSplit && i > 0 && abs(x[i]) < 1e-5)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <math.h>
+#include <string.h>
+#include "gromacs/commandline/pargs.h"
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "gromacs/utility/smalloc.h"
+#include "macros.h"
+#include "gmx_fatal.h"
+#include "vec.h"
+#include "pbc.h"
+#include "gromacs/fileio/futil.h"
+#include "index.h"
+#include "gromacs/fileio/pdbio.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/matio.h"
+#include "mshift.h"
+#include "xvgr.h"
+#include "rmpbc.h"
+#include "txtdump.h"
+#include "eigio.h"
+#include "physics.h"
+#include "gmx_ana.h"
+
+#include "gromacs/math/do_fit.h"
+
+static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
+{
+ int i;
+ double hwkT, w, dS, S = 0;
+ double hbar, lambda;
+
+ hbar = PLANCK1/(2*M_PI);
+ for (i = 0; (i < n-nskip); i++)
+ {
+ if (eigval[i] > 0)
+ {
+ lambda = eigval[i]*AMU;
+ w = sqrt(BOLTZMANN*temp/lambda)/NANO;
+ hwkT = (hbar*w)/(BOLTZMANN*temp);
+ dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
+ S += dS;
+ if (debug)
+ {
+ fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
+ i, w, lambda, hwkT, dS);
+ }
+ }
+ else
+ {
+ fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
+ w = 0;
+ }
+ }
+ fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
+ S*RGAS);
+}
+
+static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
+ real *eigval, real temp)
+{
+ double dd, deter;
+ int *indx;
+ int i, j, k, m;
+ char buf[256];
+ double hbar, kt, kteh, S;
+
+ hbar = PLANCK1/(2*M_PI);
+ kt = BOLTZMANN*temp;
+ kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
+ if (debug)
+ {
+ fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
+ }
+
+ deter = 0;
+ for (i = 0; (i < n-nskip); i++)
+ {
+ dd = 1+kteh*eigval[i];
+ deter += log(dd);
+ }
+ S = 0.5*RGAS*deter;
+
+ fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
+}
+
+const char *proj_unit;
+
+static real tick_spacing(real range, int minticks)
+{
+ real sp;
+
+ if (range <= 0)
+ {
+ return 1.0;
+ }
+
+ sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
+ while (range/sp < minticks-1)
+ {
+ sp = sp/2;
+ }
+
+ return sp;
+}
+
+static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
+ const char *title, const char *subtitle,
+ const char *xlabel, const char **ylabel,
+ int n, real *x, real **y, real ***sy,
+ real scale_x, gmx_bool bZero, gmx_bool bSplit,
+ const output_env_t oenv)
+{
+ FILE *out;
+ int g, s, i;
+ real min, max, xsp, ysp;
+
+ out = gmx_ffopen(file, "w");
+ if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
+ {
+ fprintf(out, "@ autoscale onread none\n");
+ }
+ for (g = 0; g < ngraphs; g++)
+ {
+ if (y)
+ {
+ min = y[g][0];
+ max = y[g][0];
+ for (i = 0; i < n; i++)
+ {
+ if (y[g][i] < min)
+ {
+ min = y[g][i];
+ }
+ if (y[g][i] > max)
+ {
+ max = y[g][i];
+ }
+ }
+ }
+ else
+ {
+ min = sy[g][0][0];
+ max = sy[g][0][0];
+ for (s = 0; s < nsetspergraph; s++)
+ {
+ for (i = 0; i < n; i++)
+ {
+ if (sy[g][s][i] < min)
+ {
+ min = sy[g][s][i];
+ }
+ if (sy[g][s][i] > max)
+ {
+ max = sy[g][s][i];
+ }
+ }
+ }
+ }
+ if (bZero)
+ {
+ min = 0;
+ }
+ else
+ {
+ min = min-0.1*(max-min);
+ }
+ max = max+0.1*(max-min);
+ xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
+ ysp = tick_spacing(max-min, 3);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
+ if (g == 0)
+ {
+ fprintf(out, "@ title \"%s\"\n", title);
+ if (subtitle)
+ {
+ fprintf(out, "@ subtitle \"%s\"\n", subtitle);
+ }
+ }
+ if (g == ngraphs-1)
+ {
+ fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
+ }
+ else
+ {
+ fprintf(out, "@ xaxis ticklabel off\n");
+ }
+ if (n > 1)
+ {
+ fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
+ fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
+ fprintf(out, "@ world ymin %g\n", min);
+ fprintf(out, "@ world ymax %g\n", max);
+ }
+ fprintf(out, "@ view xmin 0.15\n");
+ fprintf(out, "@ view xmax 0.85\n");
+ fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
+ fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
+ fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
+ fprintf(out, "@ xaxis tick major %g\n", xsp);
+ fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
+ fprintf(out, "@ xaxis ticklabel start type spec\n");
+ fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
+ fprintf(out, "@ yaxis tick major %g\n", ysp);
+ fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
+ fprintf(out, "@ yaxis ticklabel start type spec\n");
+ fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
+ if ((min < 0) && (max > 0))
+ {
+ fprintf(out, "@ zeroxaxis bar on\n");
+ fprintf(out, "@ zeroxaxis bar linestyle 3\n");
+ }
+ }
+ for (s = 0; s < nsetspergraph; s++)
+ {
+ for (i = 0; i < n; i++)
+ {
- if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5)
++ if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
+ {
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "&\n");
+ }
+ }
+ fprintf(out, "%10.4f %10.5f\n",
+ x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
+ }
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "&\n");
+ }
+ }
+ }
+ gmx_ffclose(out);
+}
+
+static void
+compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
+ real *eigval1, int neig1, real *eigval2, int neig2)
+{
+ int n, nrow;
+ int i, j, k;
+ double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
+
+ n = min(n1, n2);
+
+ n = min(n, min(neig1, neig2));
+ fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
+
+ sum1 = 0;
+ for (i = 0; i < n; i++)
+ {
+ if (eigval1[i] < 0)
+ {
+ eigval1[i] = 0;
+ }
+ sum1 += eigval1[i];
+ eigval1[i] = sqrt(eigval1[i]);
+ }
+ trace1 = sum1;
+ for (i = n; i < neig1; i++)
+ {
+ trace1 += eigval1[i];
+ }
+ sum2 = 0;
+ for (i = 0; i < n; i++)
+ {
+ if (eigval2[i] < 0)
+ {
+ eigval2[i] = 0;
+ }
+ sum2 += eigval2[i];
+ eigval2[i] = sqrt(eigval2[i]);
+ }
+ trace2 = sum2;
+ for (i = n; i < neig2; i++)
+ {
+ trace2 += eigval2[i];
+ }
+
+ fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
+ if (neig1 != n || neig2 != n)
+ {
+ fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
+ (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
+ }
+ fprintf(stdout, "Square root of the traces: %g and %g\n",
+ sqrt(sum1), sqrt(sum2));
+
+ sab = 0;
+ for (i = 0; i < n; i++)
+ {
+ tmp = 0;
+ for (j = 0; j < n; j++)
+ {
+ ip = 0;
+ for (k = 0; k < natoms; k++)
+ {
+ ip += iprod(eigvec1[i][k], eigvec2[j][k]);
+ }
+ tmp += eigval2[j]*ip*ip;
+ }
+ sab += eigval1[i]*tmp;
+ }
+
+ samsb2 = sum1+sum2-2*sab;
+ if (samsb2 < 0)
+ {
+ samsb2 = 0;
+ }
+
+ fprintf(stdout, "The overlap of the covariance matrices:\n");
+ fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
+ tmp = 1-sab/sqrt(sum1*sum2);
+ if (tmp < 0)
+ {
+ tmp = 0;
+ }
+ fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
+}
+
+
+static void inprod_matrix(const char *matfile, int natoms,
+ int nvec1, int *eignr1, rvec **eigvec1,
+ int nvec2, int *eignr2, rvec **eigvec2,
+ gmx_bool bSelect, int noutvec, int *outvec)
+{
+ FILE *out;
+ real **mat;
+ int i, x1, y1, x, y, nlevels;
+ int nx, ny;
+ real inp, *t_x, *t_y, max;
+ t_rgb rlo, rhi;
+
+ snew(t_y, nvec2);
+ if (bSelect)
+ {
+ nx = noutvec;
+ ny = 0;
+ for (y1 = 0; y1 < nx; y1++)
+ {
+ if (outvec[y1] < nvec2)
+ {
+ t_y[ny] = eignr2[outvec[y1]]+1;
+ ny++;
+ }
+ }
+ }
+ else
+ {
+ nx = nvec1;
+ ny = nvec2;
+ for (y = 0; y < ny; y++)
+ {
+ t_y[y] = eignr2[y]+1;
+ }
+ }
+
+ fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
+ nx, nvec2);
+
+ snew(mat, nx);
+ snew(t_x, nx);
+ max = 0;
+ for (x1 = 0; x1 < nx; x1++)
+ {
+ snew(mat[x1], ny);
+ if (bSelect)
+ {
+ x = outvec[x1];
+ }
+ else
+ {
+ x = x1;
+ }
+ t_x[x1] = eignr1[x]+1;
+ fprintf(stderr, " %d", eignr1[x]+1);
+ for (y1 = 0; y1 < ny; y1++)
+ {
+ inp = 0;
+ if (bSelect)
+ {
+ while (outvec[y1] >= nvec2)
+ {
+ y1++;
+ }
+ y = outvec[y1];
+ }
+ else
+ {
+ y = y1;
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ inp += iprod(eigvec1[x][i], eigvec2[y][i]);
+ }
+ mat[x1][y1] = fabs(inp);
+ if (mat[x1][y1] > max)
+ {
+ max = mat[x1][y1];
+ }
+ }
+ }
+ fprintf(stderr, "\n");
+ rlo.r = 1; rlo.g = 1; rlo.b = 1;
+ rhi.r = 0; rhi.g = 0; rhi.b = 0;
+ nlevels = 41;
+ out = gmx_ffopen(matfile, "w");
+ write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
+ nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
+ gmx_ffclose(out);
+}
+
+static void overlap(const char *outfile, int natoms,
+ rvec **eigvec1,
+ int nvec2, int *eignr2, rvec **eigvec2,
+ int noutvec, int *outvec,
+ const output_env_t oenv)
+{
+ FILE *out;
+ int i, v, vec, x;
+ real overlap, inp;
+
+ fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
+ for (i = 0; i < noutvec; i++)
+ {
+ fprintf(stderr, "%d ", outvec[i]+1);
+ }
+ fprintf(stderr, "\n");
+
+ out = xvgropen(outfile, "Subspace overlap",
+ "Eigenvectors of trajectory 2", "Overlap", oenv);
+ fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
+ noutvec);
+ overlap = 0;
+ for (x = 0; x < nvec2; x++)
+ {
+ for (v = 0; v < noutvec; v++)
+ {
+ vec = outvec[v];
+ inp = 0;
+ for (i = 0; i < natoms; i++)
+ {
+ inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
+ }
+ overlap += sqr(inp);
+ }
+ fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
+ }
+
+ gmx_ffclose(out);
+}
+
+static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
+ const char *projfile, const char *twodplotfile,
+ const char *threedplotfile, const char *filterfile, int skip,
+ const char *extremefile, gmx_bool bExtrAll, real extreme,
+ int nextr, t_atoms *atoms, int natoms, atom_id *index,
+ gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
+ real *sqrtm, rvec *xav,
+ int *eignr, rvec **eigvec,
+ int noutvec, int *outvec, gmx_bool bSplit,
+ const output_env_t oenv)
+{
+ FILE *xvgrout = NULL;
+ int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
+ t_trxstatus *out = NULL;
+ t_trxstatus *status;
+ int noutvec_extr, imin, imax;
+ real *pmin, *pmax;
+ atom_id *all_at;
+ matrix box;
+ rvec *xread, *x;
+ real t, inp, **inprod = NULL, min = 0, max = 0;
+ char str[STRLEN], str2[STRLEN], **ylabel, *c;
+ real fact;
+ gmx_rmpbc_t gpbc = NULL;
+
+ snew(x, natoms);
+
+ if (bExtrAll)
+ {
+ noutvec_extr = noutvec;
+ }
+ else
+ {
+ noutvec_extr = 1;
+ }
+
+
+ if (trajfile)
+ {
+ snew(inprod, noutvec+1);
+
+ if (filterfile)
+ {
+ fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
+ filterfile);
+ for (i = 0; i < noutvec; i++)
+ {
+ fprintf(stderr, "%d ", outvec[i]+1);
+ }
+ fprintf(stderr, "\n");
+ out = open_trx(filterfile, "w");
+ }
+ snew_size = 0;
+ nfr = 0;
+ nframes = 0;
+ nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
+ if (nat > atoms->nr)
+ {
+ gmx_fatal(FARGS, "the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)", nat, atoms->nr);
+ }
+ snew(all_at, nat);
+
+ if (top)
+ {
+ gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
+ }
+
+ for (i = 0; i < nat; i++)
+ {
+ all_at[i] = i;
+ }
+ do
+ {
+ if (nfr % skip == 0)
+ {
+ if (top)
+ {
+ gmx_rmpbc(gpbc, nat, box, xread);
+ }
+ if (nframes >= snew_size)
+ {
+ snew_size += 100;
+ for (i = 0; i < noutvec+1; i++)
+ {
+ srenew(inprod[i], snew_size);
+ }
+ }
+ inprod[noutvec][nframes] = t;
+ /* calculate x: a fitted struture of the selected atoms */
+ if (bFit)
+ {
+ reset_x(nfit, ifit, nat, NULL, xread, w_rls);
+ do_fit(nat, w_rls, xref, xread);
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ copy_rvec(xread[index[i]], x[i]);
+ }
+
+ for (v = 0; v < noutvec; v++)
+ {
+ vec = outvec[v];
+ /* calculate (mass-weighted) projection */
+ inp = 0;
+ for (i = 0; i < natoms; i++)
+ {
+ inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
+ eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
+ eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
+ }
+ inprod[v][nframes] = inp;
+ }
+ if (filterfile)
+ {
+ for (i = 0; i < natoms; i++)
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ /* misuse xread for output */
+ xread[index[i]][d] = xav[i][d];
+ for (v = 0; v < noutvec; v++)
+ {
+ xread[index[i]][d] +=
+ inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
+ }
+ }
+ }
+ write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
+ }
+ nframes++;
+ }
+ nfr++;
+ }
+ while (read_next_x(oenv, status, &t, xread, box));
+ close_trx(status);
+ sfree(x);
+ if (filterfile)
+ {
+ close_trx(out);
+ }
+ }
+ else
+ {
+ snew(xread, atoms->nr);
+ }
+
+ if (top)
+ {
+ gmx_rmpbc_done(gpbc);
+ }
+
+
+ if (projfile)
+ {
+ snew(ylabel, noutvec);
+ for (v = 0; v < noutvec; v++)
+ {
+ sprintf(str, "vec %d", eignr[outvec[v]]+1);
+ ylabel[v] = strdup(str);
+ }
+ sprintf(str, "projection on eigenvectors (%s)", proj_unit);
+ write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
+ (const char **)ylabel,
+ nframes, inprod[noutvec], inprod, NULL,
+ output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
+ }
+
+ if (twodplotfile)
+ {
+ sprintf(str, "projection on eigenvector %d (%s)",
+ eignr[outvec[0]]+1, proj_unit);
+ sprintf(str2, "projection on eigenvector %d (%s)",
+ eignr[outvec[noutvec-1]]+1, proj_unit);
+ xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
+ oenv);
+ for (i = 0; i < nframes; i++)
+ {
- if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5)
++ if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
+ {
+ fprintf(xvgrout, "&\n");
+ }
+ fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
+ }
+ gmx_ffclose(xvgrout);
+ }
+
+ if (threedplotfile)
+ {
+ t_atoms atoms;
+ rvec *x;
+ real *b = NULL;
+ matrix box;
+ char *resnm, *atnm, pdbform[STRLEN];
+ gmx_bool bPDB, b4D;
+ FILE *out;
+
+ if (noutvec < 3)
+ {
+ gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
+ }
+
+ /* initialize */
+ bPDB = fn2ftp(threedplotfile) == efPDB;
+ clear_mat(box);
+ box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
+
+ b4D = bPDB && (noutvec >= 4);
+ if (b4D)
+ {
+ fprintf(stderr, "You have selected four or more eigenvectors:\n"
+ "fourth eigenvector will be plotted "
+ "in bfactor field of pdb file\n");
+ sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
+ eignr[outvec[0]]+1, eignr[outvec[1]]+1,
+ eignr[outvec[2]]+1, eignr[outvec[3]]+1);
+ }
+ else
+ {
+ sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
+ eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
+ }
+ init_t_atoms(&atoms, nframes, FALSE);
+ snew(x, nframes);
+ snew(b, nframes);
+ atnm = strdup("C");
+ resnm = strdup("PRJ");
+
+ if (nframes > 10000)
+ {
+ fact = 10000.0/nframes;
+ }
+ else
+ {
+ fact = 1.0;
+ }
+
+ for (i = 0; i < nframes; i++)
+ {
+ atoms.atomname[i] = &atnm;
+ atoms.atom[i].resind = i;
+ atoms.resinfo[i].name = &resnm;
+ atoms.resinfo[i].nr = ceil(i*fact);
+ atoms.resinfo[i].ic = ' ';
+ x[i][XX] = inprod[0][i];
+ x[i][YY] = inprod[1][i];
+ x[i][ZZ] = inprod[2][i];
+ if (b4D)
+ {
+ b[i] = inprod[3][i];
+ }
+ }
+ if ( ( b4D || bSplit ) && bPDB)
+ {
+ strcpy(pdbform, get_pdbformat());
+ strcat(pdbform, "%8.4f%8.4f\n");
+
+ out = gmx_ffopen(threedplotfile, "w");
+ fprintf(out, "HEADER %s\n", str);
+ if (b4D)
+ {
+ fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
+ }
+ j = 0;
+ for (i = 0; i < atoms.nr; i++)
+ {
++ if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
+ {
+ fprintf(out, "TER\n");
+ j = 0;
+ }
+ fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
+ 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i]);
+ if (j > 0)
+ {
+ fprintf(out, "CONECT%5d%5d\n", i, i+1);
+ }
+ j++;
+ }
+ fprintf(out, "TER\n");
+ gmx_ffclose(out);
+ }
+ else
+ {
+ write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
+ }
+ free_t_atoms(&atoms, FALSE);
+ }
+
+ if (extremefile)
+ {
+ snew(pmin, noutvec_extr);
+ snew(pmax, noutvec_extr);
+ if (extreme == 0)
+ {
+ fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
+ fprintf(stderr,
+ "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
+ imin = 0;
+ imax = 0;
+ for (v = 0; v < noutvec_extr; v++)
+ {
+ for (i = 0; i < nframes; i++)
+ {
+ if (inprod[v][i] < inprod[v][imin])
+ {
+ imin = i;
+ }
+ if (inprod[v][i] > inprod[v][imax])
+ {
+ imax = i;
+ }
+ }
+ pmin[v] = inprod[v][imin];
+ pmax[v] = inprod[v][imax];
+ fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
+ eignr[outvec[v]]+1,
+ pmin[v], imin, pmax[v], imax);
+ }
+ }
+ else
+ {
+ pmin[0] = -extreme;
+ pmax[0] = extreme;
+ }
+ /* build format string for filename: */
+ strcpy(str, extremefile); /* copy filename */
+ c = strrchr(str, '.'); /* find where extention begins */
+ strcpy(str2, c); /* get extention */
+ sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
+ for (v = 0; v < noutvec_extr; v++)
+ {
+ /* make filename using format string */
+ if (noutvec_extr == 1)
+ {
+ strcpy(str2, extremefile);
+ }
+ else
+ {
+ sprintf(str2, str, eignr[outvec[v]]+1);
+ }
+ fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
+ nextr, outvec[v]+1, str2);
+ out = open_trx(str2, "w");
+ for (frame = 0; frame < nextr; frame++)
+ {
+ if ((extreme == 0) && (nextr <= 3))
+ {
+ for (i = 0; i < natoms; i++)
+ {
+ atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
+ }
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ xread[index[i]][d] =
+ (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
+ *eigvec[outvec[v]][i][d]/sqrtm[i]);
+ }
+ }
+ write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
+ }
+ close_trx(out);
+ }
+ sfree(pmin);
+ sfree(pmax);
+ }
+ fprintf(stderr, "\n");
+}
+
+static void components(const char *outfile, int natoms,
+ int *eignr, rvec **eigvec,
+ int noutvec, int *outvec,
+ const output_env_t oenv)
+{
+ int g, s, v, i;
+ real *x, ***y;
+ char str[STRLEN], **ylabel;
+
+ fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
+
+ snew(ylabel, noutvec);
+ snew(y, noutvec);
+ snew(x, natoms);
+ for (i = 0; i < natoms; i++)
+ {
+ x[i] = i+1;
+ }
+ for (g = 0; g < noutvec; g++)
+ {
+ v = outvec[g];
+ sprintf(str, "vec %d", eignr[v]+1);
+ ylabel[g] = strdup(str);
+ snew(y[g], 4);
+ for (s = 0; s < 4; s++)
+ {
+ snew(y[g][s], natoms);
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ y[g][0][i] = norm(eigvec[v][i]);
+ for (s = 0; s < 3; s++)
+ {
+ y[g][s+1][i] = eigvec[v][i][s];
+ }
+ }
+ }
+ write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
+ "black: total, red: x, green: y, blue: z",
+ "Atom number", (const char **)ylabel,
+ natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
+ fprintf(stderr, "\n");
+}
+
+static void rmsf(const char *outfile, int natoms, real *sqrtm,
+ int *eignr, rvec **eigvec,
+ int noutvec, int *outvec,
+ real *eigval, int neig,
+ const output_env_t oenv)
+{
+ int nrow, g, v, i;
+ real *x, **y;
+ char str[STRLEN], **ylabel;
+
+ for (i = 0; i < neig; i++)
+ {
+ if (eigval[i] < 0)
+ {
+ eigval[i] = 0;
+ }
+ }
+
+ fprintf(stderr, "Writing rmsf to %s\n", outfile);
+
+ snew(ylabel, noutvec);
+ snew(y, noutvec);
+ snew(x, natoms);
+ for (i = 0; i < natoms; i++)
+ {
+ x[i] = i+1;
+ }
+ for (g = 0; g < noutvec; g++)
+ {
+ v = outvec[g];
+ if (eignr[v] >= neig)
+ {
+ gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
+ }
+ sprintf(str, "vec %d", eignr[v]+1);
+ ylabel[g] = strdup(str);
+ snew(y[g], natoms);
+ for (i = 0; i < natoms; i++)
+ {
+ y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
+ }
+ }
+ write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
+ "Atom number", (const char **)ylabel,
+ natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
+ fprintf(stderr, "\n");
+}
+
+int gmx_anaeig(int argc, char *argv[])
+{
+ static const char *desc[] = {
+ "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
+ "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
+ "([gmx-nmeig]).[PAR]",
+
+ "When a trajectory is projected on eigenvectors, all structures are",
+ "fitted to the structure in the eigenvector file, if present, otherwise",
+ "to the structure in the structure file. When no run input file is",
+ "supplied, periodicity will not be taken into account. Most analyses",
+ "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
+ "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
+
+ "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
+ "[TT]-first[tt] to [TT]-last[tt].[PAR]",
+
+ "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
+ "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
+
+ "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
+ "[TT]-first[tt] to [TT]-last[tt].",
+ "The projections of a trajectory on the eigenvectors of its",
+ "covariance matrix are called principal components (pc's).",
+ "It is often useful to check the cosine content of the pc's,",
+ "since the pc's of random diffusion are cosines with the number",
+ "of periods equal to half the pc index.",
+ "The cosine content of the pc's can be calculated with the program",
+ "[gmx-analyze].[PAR]",
+
+ "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
+ "[TT]-first[tt] and [TT]-last[tt].[PAR]",
+
+ "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
+ "three selected eigenvectors.[PAR]",
+
+ "[TT]-filt[tt]: filter the trajectory to show only the motion along",
+ "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
+
+ "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
+ "on the average structure and interpolate [TT]-nframes[tt] frames",
+ "between them, or set your own extremes with [TT]-max[tt]. The",
+ "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
+ "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
+ "will be written to separate files. Chain identifiers will be added",
+ "when writing a [TT].pdb[tt] file with two or three structures (you",
+ "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
+
+ " Overlap calculations between covariance analysis:[BR]",
+ " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
+
+ "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
+ "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
+ "in file [TT]-v[tt].[PAR]",
+
+ "[TT]-inpr[tt]: calculate a matrix of inner-products between",
+ "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
+ "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
+ "have been set explicitly.[PAR]",
+
+ "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
+ "a single number for the overlap between the covariance matrices is",
+ "generated. The formulas are:[BR]",
+ " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
+ "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
+ " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
+ "where M1 and M2 are the two covariance matrices and tr is the trace",
+ "of a matrix. The numbers are proportional to the overlap of the square",
+ "root of the fluctuations. The normalized overlap is the most useful",
+ "number, it is 1 for identical matrices and 0 when the sampled",
+ "subspaces are orthogonal.[PAR]",
+ "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
+ "computed based on the Quasiharmonic approach and based on",
+ "Schlitter's formula."
+ };
+ static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
+ static real max = 0.0, temp = 298.15;
+ static gmx_bool bSplit = FALSE, bEntropy = FALSE;
+ t_pargs pa[] = {
+ { "-first", FALSE, etINT, {&first},
+ "First eigenvector for analysis (-1 is select)" },
+ { "-last", FALSE, etINT, {&last},
+ "Last eigenvector for analysis (-1 is till the last)" },
+ { "-skip", FALSE, etINT, {&skip},
+ "Only analyse every nr-th frame" },
+ { "-max", FALSE, etREAL, {&max},
+ "Maximum for projection of the eigenvector on the average structure, "
+ "max=0 gives the extremes" },
+ { "-nframes", FALSE, etINT, {&nextr},
+ "Number of frames for the extremes output" },
+ { "-split", FALSE, etBOOL, {&bSplit},
+ "Split eigenvector projections where time is zero" },
+ { "-entropy", FALSE, etBOOL, {&bEntropy},
+ "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
+ { "-temp", FALSE, etREAL, {&temp},
+ "Temperature for entropy calculations" },
+ { "-nevskip", FALSE, etINT, {&nskip},
+ "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." }
+ };
+#define NPA asize(pa)
+
+ FILE *out;
+ int status, trjout;
+ t_topology top;
+ int ePBC = -1;
+ t_atoms *atoms = NULL;
+ rvec *xtop, *xref1, *xref2, *xrefp = NULL;
+ gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
+ int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
+ rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
+ matrix topbox;
+ real xid, totmass, *sqrtm, *w_rls, t, lambda;
+ int natoms, step;
+ char *grpname;
+ const char *indexfile;
+ char title[STRLEN];
+ int i, j, d;
+ int nout, *iout, noutvec, *outvec, nfit;
+ atom_id *index, *ifit;
+ const char *VecFile, *Vec2File, *topfile;
+ const char *EigFile, *Eig2File;
+ const char *CompFile, *RmsfFile, *ProjOnVecFile;
+ const char *TwoDPlotFile, *ThreeDPlotFile;
+ const char *FilterFile, *ExtremeFile;
+ const char *OverlapFile, *InpMatFile;
+ gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
+ gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
+ real *eigval1 = NULL, *eigval2 = NULL;
+ int neig1, neig2;
+ double **xvgdata;
+ output_env_t oenv;
+ gmx_rmpbc_t gpbc;
+
+ t_filenm fnm[] = {
+ { efTRN, "-v", "eigenvec", ffREAD },
+ { efTRN, "-v2", "eigenvec2", ffOPTRD },
+ { efTRX, "-f", NULL, ffOPTRD },
+ { efTPS, NULL, NULL, ffOPTRD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efXVG, "-eig", "eigenval", ffOPTRD },
+ { efXVG, "-eig2", "eigenval2", ffOPTRD },
+ { efXVG, "-comp", "eigcomp", ffOPTWR },
+ { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
+ { efXVG, "-proj", "proj", ffOPTWR },
+ { efXVG, "-2d", "2dproj", ffOPTWR },
+ { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
+ { efTRX, "-filt", "filtered", ffOPTWR },
+ { efTRX, "-extr", "extreme.pdb", ffOPTWR },
+ { efXVG, "-over", "overlap", ffOPTWR },
+ { efXPM, "-inpr", "inprod", ffOPTWR }
+ };
+#define NFILE asize(fnm)
+
+ if (!parse_common_args(&argc, argv,
+ PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
+ NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+
+ indexfile = ftp2fn_null(efNDX, NFILE, fnm);
+
+ VecFile = opt2fn("-v", NFILE, fnm);
+ Vec2File = opt2fn_null("-v2", NFILE, fnm);
+ topfile = ftp2fn(efTPS, NFILE, fnm);
+ EigFile = opt2fn_null("-eig", NFILE, fnm);
+ Eig2File = opt2fn_null("-eig2", NFILE, fnm);
+ CompFile = opt2fn_null("-comp", NFILE, fnm);
+ RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
+ ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
+ TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
+ ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
+ FilterFile = opt2fn_null("-filt", NFILE, fnm);
+ ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
+ OverlapFile = opt2fn_null("-over", NFILE, fnm);
+ InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
+
+ bTop = fn2bTPX(topfile);
+ bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
+ || FilterFile || ExtremeFile;
+ bFirstLastSet =
+ opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
+ bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
+ OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
+ bVec2 = Vec2File || OverlapFile || InpMatFile;
+ bM = RmsfFile || bProj;
+ bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
+ || TwoDPlotFile || ThreeDPlotFile;
+ bIndex = bM || bProj;
+ bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
+ FilterFile || (bIndex && indexfile);
+ bCompare = Vec2File || Eig2File;
+ bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
+
+ read_eigenvectors(VecFile, &natoms, &bFit1,
+ &xref1, &bDMR1, &xav1, &bDMA1,
+ &nvec1, &eignr1, &eigvec1, &eigval1);
+ neig1 = DIM*natoms;
+
+ /* Overwrite eigenvalues from separate files if the user provides them */
+ if (EigFile != NULL)
+ {
+ int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
+ if (neig_tmp != neig1)
+ {
+ fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
+ }
+ neig1 = neig_tmp;
+ srenew(eigval1, neig1);
+ for (j = 0; j < neig1; j++)
+ {
+ real tmp = eigval1[j];
+ eigval1[j] = xvgdata[1][j];
+ if (debug && (eigval1[j] != tmp))
+ {
+ fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
+ j, tmp, eigval1[j]);
+ }
+ }
+ for (j = 0; j < i; j++)
+ {
+ sfree(xvgdata[j]);
+ }
+ sfree(xvgdata);
+ fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
+ }
+
+ if (bEntropy)
+ {
+ if (bDMA1)
+ {
+ gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
+ }
+ calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
+ calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
+ }
+
+ if (bVec2)
+ {
+ if (!Vec2File)
+ {
+ gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
+ }
+ read_eigenvectors(Vec2File, &neig2, &bFit2,
+ &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
+
+ neig2 = DIM*neig2;
+ if (neig2 != neig1)
+ {
+ gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
+ }
+ }
+
+ if (Eig2File != NULL)
+ {
+ neig2 = read_xvg(Eig2File, &xvgdata, &i);
+ srenew(eigval2, neig2);
+ for (j = 0; j < neig2; j++)
+ {
+ eigval2[j] = xvgdata[1][j];
+ }
+ for (j = 0; j < i; j++)
+ {
+ sfree(xvgdata[j]);
+ }
+ sfree(xvgdata);
+ fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
+ }
+
+
+ if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
+ {
+ bM = FALSE;
+ }
+ if ((xref1 == NULL) && (bM || bTraj))
+ {
+ bTPS = TRUE;
+ }
+
+ xtop = NULL;
+ nfit = 0;
+ ifit = NULL;
+ w_rls = NULL;
+
+ if (!bTPS)
+ {
+ bTop = FALSE;
+ }
+ else
+ {
+ bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
+ title, &top, &ePBC, &xtop, NULL, topbox, bM);
+ atoms = &top.atoms;
+ gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
+ gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
+ /* Fitting is only required for the projection */
+ if (bProj && bFit1)
+ {
+ if (xref1 == NULL)
+ {
+ printf("\nNote: the structure in %s should be the same\n"
+ " as the one used for the fit in g_covar\n", topfile);
+ }
+ printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
+ get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
+
+ snew(w_rls, atoms->nr);
+ for (i = 0; (i < nfit); i++)
+ {
+ if (bDMR1)
+ {
+ w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
+ }
+ else
+ {
+ w_rls[ifit[i]] = 1.0;
+ }
+ }
+
+ snew(xrefp, atoms->nr);
+ if (xref1 != NULL)
+ {
+ /* Safety check between selected fit-group and reference structure read from the eigenvector file */
+ if (natoms != nfit)
+ {
+ gmx_fatal(FARGS, "you selected a group with %d elements instead of %d, your selection does not fit the reference structure in the eigenvector file.", nfit, natoms);
+ }
+ for (i = 0; (i < nfit); i++)
+ {
+ copy_rvec(xref1[i], xrefp[ifit[i]]);
+ }
+ }
+ else
+ {
+ /* The top coordinates are the fitting reference */
+ for (i = 0; (i < nfit); i++)
+ {
+ copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
+ }
+ reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
+ }
+ }
+ gmx_rmpbc_done(gpbc);
+ }
+
+ if (bIndex)
+ {
+ printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
+ get_index(atoms, indexfile, 1, &i, &index, &grpname);
+ if (i != natoms)
+ {
+ gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
+ }
+ printf("\n");
+ }
+
+ snew(sqrtm, natoms);
+ if (bM && bDMA1)
+ {
+ proj_unit = "u\\S1/2\\Nnm";
+ for (i = 0; (i < natoms); i++)
+ {
+ sqrtm[i] = sqrt(atoms->atom[index[i]].m);
+ }
+ }
+ else
+ {
+ proj_unit = "nm";
+ for (i = 0; (i < natoms); i++)
+ {
+ sqrtm[i] = 1.0;
+ }
+ }
+
+ if (bVec2)
+ {
+ t = 0;
+ totmass = 0;
+ for (i = 0; (i < natoms); i++)
+ {
+ for (d = 0; (d < DIM); d++)
+ {
+ t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
+ totmass += sqr(sqrtm[i]);
+ }
+ }
+ fprintf(stdout, "RMSD (without fit) between the two average structures:"
+ " %.3f (nm)\n\n", sqrt(t/totmass));
+ }
+
+ if (last == -1)
+ {
+ last = natoms*DIM;
+ }
+ if (first > -1)
+ {
+ if (bFirstToLast)
+ {
+ /* make an index from first to last */
+ nout = last-first+1;
+ snew(iout, nout);
+ for (i = 0; i < nout; i++)
+ {
+ iout[i] = first-1+i;
+ }
+ }
+ else if (ThreeDPlotFile)
+ {
+ /* make an index of first+(0,1,2) and last */
+ nout = bPDB3D ? 4 : 3;
+ nout = min(last-first+1, nout);
+ snew(iout, nout);
+ iout[0] = first-1;
+ iout[1] = first;
+ if (nout > 3)
+ {
+ iout[2] = first+1;
+ }
+ iout[nout-1] = last-1;
+ }
+ else
+ {
+ /* make an index of first and last */
+ nout = 2;
+ snew(iout, nout);
+ iout[0] = first-1;
+ iout[1] = last-1;
+ }
+ }
+ else
+ {
+ printf("Select eigenvectors for output, end your selection with 0\n");
+ nout = -1;
+ iout = NULL;
+
+ do
+ {
+ nout++;
+ srenew(iout, nout+1);
+ if (1 != scanf("%d", &iout[nout]))
+ {
+ gmx_fatal(FARGS, "Error reading user input");
+ }
+ iout[nout]--;
+ }
+ while (iout[nout] >= 0);
+
+ printf("\n");
+ }
+ /* make an index of the eigenvectors which are present */
+ snew(outvec, nout);
+ noutvec = 0;
+ for (i = 0; i < nout; i++)
+ {
+ j = 0;
+ while ((j < nvec1) && (eignr1[j] != iout[i]))
+ {
+ j++;
+ }
+ if ((j < nvec1) && (eignr1[j] == iout[i]))
+ {
+ outvec[noutvec] = j;
+ noutvec++;
+ }
+ }
+ fprintf(stderr, "%d eigenvectors selected for output", noutvec);
+ if (noutvec <= 100)
+ {
+ fprintf(stderr, ":");
+ for (j = 0; j < noutvec; j++)
+ {
+ fprintf(stderr, " %d", eignr1[outvec[j]]+1);
+ }
+ }
+ fprintf(stderr, "\n");
+
+ if (CompFile)
+ {
+ components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
+ }
+
+ if (RmsfFile)
+ {
+ rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
+ neig1, oenv);
+ }
+
+ if (bProj)
+ {
+ project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
+ bTop ? &top : NULL, ePBC, topbox,
+ ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
+ ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
+ bFit1, xrefp, nfit, ifit, w_rls,
+ sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
+ oenv);
+ }
+
+ if (OverlapFile)
+ {
+ overlap(OverlapFile, natoms,
+ eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
+ }
+
+ if (InpMatFile)
+ {
+ inprod_matrix(InpMatFile, natoms,
+ nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
+ bFirstLastSet, noutvec, outvec);
+ }
+
+ if (bCompare)
+ {
+ compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
+ }
+
+
+ if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
+ !bCompare && !bEntropy)
+ {
+ fprintf(stderr, "\nIf you want some output,"
+ " set one (or two or ...) of the output file options\n");
+ }
+
+
+ view_all(oenv, NFILE, fnm);
+
+ return 0;
+}
--- /dev/null
- if (abs(qall) > 0.01)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "gromacs/commandline/pargs.h"
+#include "typedefs.h"
+#include "gromacs/utility/smalloc.h"
+#include "vec.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trxio.h"
+#include "xvgr.h"
+#include "rmpbc.h"
+#include "pbc.h"
+#include "physics.h"
+#include "index.h"
+#include "gromacs/statistics/statistics.h"
+#include "gmx_ana.h"
+#include "macros.h"
+
+#include "gromacs/legacyheaders/gmx_fatal.h"
+
+#define SQR(x) (pow(x, 2.0))
+#define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
+
+static void index_atom2mol(int *n, int *index, t_block *mols)
+{
+ int nat, i, nmol, mol, j;
+
+ nat = *n;
+ i = 0;
+ nmol = 0;
+ mol = 0;
+ while (i < nat)
+ {
+ while (index[i] > mols->index[mol])
+ {
+ mol++;
+ if (mol >= mols->nr)
+ {
+ gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
+ }
+ }
+ for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
+ {
+ if (i >= nat || index[i] != j)
+ {
+ gmx_fatal(FARGS, "The index group does not consist of whole molecules");
+ }
+ i++;
+ }
+ index[nmol++] = mol;
+ }
+
+ fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
+
+ *n = nmol;
+}
+
+static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
+{
+
+ real mtot;
+ real qtot;
+ real qall;
+ int i, j, k, l;
+ int ai, ci;
+ gmx_bool bNEU;
+ ai = 0;
+ ci = 0;
+ qall = 0.0;
+
+
+
+ for (i = 0; i < top.mols.nr; i++)
+ {
+ k = top.mols.index[i];
+ l = top.mols.index[i+1];
+ mtot = 0.0;
+ qtot = 0.0;
+
+ for (j = k; j < l; j++)
+ {
+ mtot += top.atoms.atom[j].m;
+ qtot += top.atoms.atom[j].q;
+ }
+
+ for (j = k; j < l; j++)
+ {
+ top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
+ mass2[j] = top.atoms.atom[j].m/mtot;
+ qmol[j] = qtot;
+ }
+
+
+ qall += qtot;
+
+ if (qtot < 0.0)
+ {
+ ai++;
+ }
+ if (qtot > 0.0)
+ {
+ ci++;
+ }
+ }
+
++ if (fabs(qall) > 0.01)
+ {
+ printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
+ bNEU = FALSE;
+ }
+ else
+ {
+ bNEU = TRUE;
+ }
+
+ return bNEU;
+
+}
+
+static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
+{
+
+ rvec hbox;
+ int d, i, m;
+
+ for (d = 0; d < DIM; d++)
+ {
+ hbox[d] = 0.5*box[d][d];
+ }
+ for (i = 0; i < natoms; i++)
+ {
+ for (m = DIM-1; m >= 0; m--)
+ {
+ while (x[i][m]-xp[i][m] <= -hbox[m])
+ {
+ for (d = 0; d <= m; d++)
+ {
+ x[i][d] += box[m][d];
+ }
+ }
+ while (x[i][m]-xp[i][m] > hbox[m])
+ {
+ for (d = 0; d <= m; d++)
+ {
+ x[i][d] -= box[m][d];
+ }
+ }
+ }
+ }
+}
+
+static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
+ rvec fr[], rvec mj, real mass2[], real qmol[])
+{
+
+ int i, j, k, l;
+ rvec tmp;
+ rvec center;
+ rvec mt1, mt2;
+ t_pbc pbc;
+
+
+ calc_box_center(ecenterRECT, box, center);
+
+ if (!bNoJump)
+ {
+ set_pbc(&pbc, ePBC, box);
+ }
+
+ clear_rvec(tmp);
+
+
+ for (i = 0; i < isize; i++)
+ {
+ clear_rvec(mt1);
+ clear_rvec(mt2);
+ k = top.mols.index[index0[i]];
+ l = top.mols.index[index0[i+1]];
+ for (j = k; j < l; j++)
+ {
+ svmul(mass2[j], fr[j], tmp);
+ rvec_inc(mt1, tmp);
+ }
+
+ if (bNoJump)
+ {
+ svmul(qmol[k], mt1, mt1);
+ }
+ else
+ {
+ pbc_dx(&pbc, mt1, center, mt2);
+ svmul(qmol[k], mt2, mt1);
+ }
+
+ rvec_inc(mj, mt1);
+
+ }
+
+}
+
+
+static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
+{
+
+ /* bCOR determines if the correlation is computed via
+ * static properties (FALSE) or the correlation integral (TRUE)
+ */
+
+ real epsilon = 0.0;
+
+
+ if (bCOR)
+ {
+ epsilon = md2-2.0*cor+mj2;
+ }
+ else
+ {
+ epsilon = md2+mj2+2.0*cor;
+ }
+
+ if (eps_rf == 0.0)
+ {
+ epsilon = 1.0+prefactor*epsilon;
+
+ }
+ else
+ {
+ epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
+ epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
+ }
+
+
+ return epsilon;
+
+}
+
+
+static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
+{
+
+ int i;
+ real corint;
+ real deltat = 0.0;
+ real rfr;
+ real sigma = 0.0;
+ real sigma_ret = 0.0;
+ corint = 0.0;
+
+ if (nfr > 1)
+ {
+ i = 0;
+
+ while (i < nfr)
+ {
+
+ rfr = (real) (nfr/nshift-i/nshift);
+ cacf[i] /= rfr;
+
+ if (time[vfr[i]] <= time[vfr[ei]])
+ {
+ sigma_ret = sigma;
+ }
+
+ fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
+
+ if ((i+1) < nfr)
+ {
+ deltat = (time[vfr[i+1]]-time[vfr[i]]);
+ }
+ corint = 2.0*deltat*cacf[i]*prefactor;
+ if (i == 0 || (i+1) == nfr)
+ {
+ corint *= 0.5;
+ }
+ sigma += corint;
+
+ i++;
+ }
+
+ }
+ else
+ {
+ printf("Too less points.\n");
+ }
+
+ return sigma_ret;
+
+}
+
+static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
+{
+
+ int i;
+ real rtmp;
+ real rfr;
+
+
+ fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
+
+
+
+ for (i = 0; i < nfr; i++)
+ {
+
+ if (refr[i] != 0.0)
+ {
+ dsp2[i] *= prefactor/refr[i];
+ fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
+ }
+
+
+ }
+
+}
+
+
+static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
+ FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
+ int ePBC, t_topology top, t_trxframe fr, real temp,
+ real trust, real bfit, real efit, real bvit, real evit,
+ t_trxstatus *status, int isize, int nmols, int nshift,
+ atom_id *index0, int indexm[], real mass2[],
+ real qmol[], real eps_rf, const output_env_t oenv)
+{
+ int i, j, k, l, f;
+ int valloc, nalloc, nfr, nvfr, m, itrust = 0;
+ int vshfr;
+ real *xshfr = NULL;
+ int *vfr = NULL;
+ real refr = 0.0;
+ real rfr = 0.0;
+ real *cacf = NULL;
+ real *time = NULL;
+ real *djc = NULL;
+ real corint = 0.0;
+ real prefactorav = 0.0;
+ real prefactor = 0.0;
+ real volume;
+ real volume_av = 0.0;
+ real dk_s, dk_d, dk_f;
+ real dm_s, dm_d;
+ real mj = 0.0;
+ real mj2 = 0.0;
+ real mjd = 0.0;
+ real mjdav = 0.0;
+ real md2 = 0.0;
+ real mdav2 = 0.0;
+ real sgk;
+ rvec mja_tmp;
+ rvec mjd_tmp;
+ rvec mdvec;
+ rvec *mu = NULL;
+ rvec *xp = NULL;
+ rvec *v0 = NULL;
+ rvec *mjdsp = NULL;
+ real *dsp2 = NULL;
+ real t0;
+ real rtmp;
+ real qtmp;
+
+
+
+ rvec tmp;
+ rvec *mtrans = NULL;
+
+ /*
+ * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
+ */
+
+ int bi, ei, ie, ii;
+ real rest = 0.0;
+ real sigma = 0.0;
+ real malt = 0.0;
+ real err = 0.0;
+ real *xfit;
+ real *yfit;
+ gmx_rmpbc_t gpbc = NULL;
+
+ /*
+ * indices for EH
+ */
+
+ ei = 0;
+ bi = 0;
+
+ /*
+ * indices for GK
+ */
+
+ ii = 0;
+ ie = 0;
+ t0 = 0;
+ sgk = 0.0;
+
+
+ /* This is the main loop over frames */
+
+
+ nfr = 0;
+
+
+ nvfr = 0;
+ vshfr = 0;
+ nalloc = 0;
+ valloc = 0;
+
+ clear_rvec(mja_tmp);
+ clear_rvec(mjd_tmp);
+ clear_rvec(mdvec);
+ clear_rvec(tmp);
+ gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
+
+ do
+ {
+
+ refr = (real)(nfr+1);
+
+ if (nfr >= nalloc)
+ {
+ nalloc += 100;
+ srenew(time, nalloc);
+ srenew(mu, nalloc);
+ srenew(mjdsp, nalloc);
+ srenew(dsp2, nalloc);
+ srenew(mtrans, nalloc);
+ srenew(xshfr, nalloc);
+
+ for (i = nfr; i < nalloc; i++)
+ {
+ clear_rvec(mjdsp[i]);
+ clear_rvec(mu[i]);
+ clear_rvec(mtrans[i]);
+ dsp2[i] = 0.0;
+ xshfr[i] = 0.0;
+ }
+ }
+
+
+
+ if (nfr == 0)
+ {
+ t0 = fr.time;
+
+ }
+
+ time[nfr] = fr.time-t0;
+
+ if (time[nfr] <= bfit)
+ {
+ bi = nfr;
+ }
+ if (time[nfr] <= efit)
+ {
+ ei = nfr;
+ }
+
+ if (bNoJump)
+ {
+
+ if (xp)
+ {
+ remove_jump(fr.box, fr.natoms, xp, fr.x);
+ }
+ else
+ {
+ snew(xp, fr.natoms);
+ }
+
+ for (i = 0; i < fr.natoms; i++)
+ {
+ copy_rvec(fr.x[i], xp[i]);
+ }
+
+ }
+
+ gmx_rmpbc_trxfr(gpbc, &fr);
+
+ calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
+
+ for (i = 0; i < isize; i++)
+ {
+ j = index0[i];
+ svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
+ rvec_inc(mu[nfr], fr.x[j]);
+ }
+
+ /*if(mod(nfr,nshift)==0){*/
+ if (nfr%nshift == 0)
+ {
+ for (j = nfr; j >= 0; j--)
+ {
+ rvec_sub(mtrans[nfr], mtrans[j], tmp);
+ dsp2[nfr-j] += norm2(tmp);
+ xshfr[nfr-j] += 1.0;
+ }
+ }
+
+ if (fr.bV)
+ {
+ if (nvfr >= valloc)
+ {
+ valloc += 100;
+ srenew(vfr, valloc);
+ if (bINT)
+ {
+ srenew(djc, valloc);
+ }
+ srenew(v0, valloc);
+ if (bACF)
+ {
+ srenew(cacf, valloc);
+ }
+ }
+ if (time[nfr] <= bvit)
+ {
+ ii = nvfr;
+ }
+ if (time[nfr] <= evit)
+ {
+ ie = nvfr;
+ }
+ vfr[nvfr] = nfr;
+ clear_rvec(v0[nvfr]);
+ if (bACF)
+ {
+ cacf[nvfr] = 0.0;
+ }
+ if (bINT)
+ {
+ djc[nvfr] = 0.0;
+ }
+ for (i = 0; i < isize; i++)
+ {
+ j = index0[i];
+ svmul(mass2[j], fr.v[j], fr.v[j]);
+ svmul(qmol[j], fr.v[j], fr.v[j]);
+ rvec_inc(v0[nvfr], fr.v[j]);
+ }
+
+ fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
+ if (bACF || bINT)
+ {
+ /*if(mod(nvfr,nshift)==0){*/
+ if (nvfr%nshift == 0)
+ {
+ for (j = nvfr; j >= 0; j--)
+ {
+ if (bACF)
+ {
+ cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
+ }
+ if (bINT)
+ {
+ djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
+ }
+ }
+ vshfr++;
+ }
+ }
+ nvfr++;
+ }
+
+ volume = det(fr.box);
+ volume_av += volume;
+
+ rvec_inc(mja_tmp, mtrans[nfr]);
+ mjd += iprod(mu[nfr], mtrans[nfr]);
+ rvec_inc(mdvec, mu[nfr]);
+
+ mj2 += iprod(mtrans[nfr], mtrans[nfr]);
+ md2 += iprod(mu[nfr], mu[nfr]);
+
+ fprintf(fmj, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mtrans[nfr][XX], mtrans[nfr][YY], mtrans[nfr][ZZ], mj2/refr, norm(mja_tmp)/refr);
+ fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
+ time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
+
+ nfr++;
+
+ }
+ while (read_next_frame(oenv, status, &fr));
+
+ gmx_rmpbc_done(gpbc);
+
+ volume_av /= refr;
+
+ prefactor = 1.0;
+ prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
+
+
+ prefactorav = E_CHARGE*E_CHARGE;
+ prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
+
+ fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
+
+ calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
+
+ /*
+ * Now we can average and calculate the correlation functions
+ */
+
+
+ mj2 /= refr;
+ mjd /= refr;
+ md2 /= refr;
+
+ svmul(1.0/refr, mdvec, mdvec);
+ svmul(1.0/refr, mja_tmp, mja_tmp);
+
+ mdav2 = norm2(mdvec);
+ mj = norm2(mja_tmp);
+ mjdav = iprod(mdvec, mja_tmp);
+
+
+ printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mja_tmp[XX], mja_tmp[YY], mja_tmp[ZZ], mj2);
+ printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mdvec[XX], mdvec[YY], mdvec[ZZ], md2);
+
+ if (v0 != NULL)
+ {
+ trust *= (real) nvfr;
+ itrust = (int) trust;
+
+ if (bINT)
+ {
+
+ printf("\nCalculating M_D - current correlation integral ... \n");
+ corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
+
+ }
+
+ if (bACF)
+ {
+
+ printf("\nCalculating current autocorrelation ... \n");
+ sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
+
+ if (ie > ii)
+ {
+
+ snew(xfit, ie-ii+1);
+ snew(yfit, ie-ii+1);
+
+ for (i = ii; i <= ie; i++)
+ {
+
+ xfit[i-ii] = log(time[vfr[i]]);
+ rtmp = fabs(cacf[i]);
+ yfit[i-ii] = log(rtmp);
+
+ }
+
+ lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
+
+ malt = exp(malt);
+
+ sigma += 1.0;
+
+ malt *= prefactorav*2.0e12/sigma;
+
+ sfree(xfit);
+ sfree(yfit);
+
+ }
+ }
+ }
+
+
+ /* Calculation of the dielectric constant */
+
+ fprintf(stderr, "\n********************************************\n");
+ dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
+ fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
+ fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
+ fprintf(stderr, "********************************************\n");
+
+
+ dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
+ fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
+ fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
+ fprintf(stderr, "\n********************************************\n");
+ if (bINT)
+ {
+ dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
+ fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
+ fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
+ }
+
+ fprintf(stderr, "\n***************************************************");
+ fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
+ fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
+
+
+
+ if (bACF)
+ {
+ fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
+ fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
+ }
+
+ if (ei > bi)
+ {
+ fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
+ fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
+
+ snew(xfit, ei-bi+1);
+ snew(yfit, ei-bi+1);
+
+ for (i = bi; i <= ei; i++)
+ {
+ xfit[i-bi] = time[i];
+ yfit[i-bi] = dsp2[i];
+ }
+
+ lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
+
+ sigma *= 1e12;
+ dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
+
+ fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
+ fprintf(stderr, "sigma=%.4f\n", sigma);
+ fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
+ fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
+
+ sfree(xfit);
+ sfree(yfit);
+ }
+ else
+ {
+ fprintf(stderr, "Too less points for a fit.\n");
+ }
+
+
+ if (v0 != NULL)
+ {
+ sfree(v0);
+ }
+ if (bACF)
+ {
+ sfree(cacf);
+ }
+ if (bINT)
+ {
+ sfree(djc);
+ }
+
+ sfree(time);
+
+
+ sfree(mjdsp);
+ sfree(mu);
+}
+
+
+
+int gmx_current(int argc, char *argv[])
+{
+
+ static int nshift = 1000;
+ static real temp = 300.0;
+ static real trust = 0.25;
+ static real eps_rf = 0.0;
+ static gmx_bool bNoJump = TRUE;
+ static real bfit = 100.0;
+ static real bvit = 0.5;
+ static real efit = 400.0;
+ static real evit = 5.0;
+ t_pargs pa[] = {
+ { "-sh", FALSE, etINT, {&nshift},
+ "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
+ { "-nojump", FALSE, etBOOL, {&bNoJump},
+ "Removes jumps of atoms across the box."},
+ { "-eps", FALSE, etREAL, {&eps_rf},
+ "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
+ { "-bfit", FALSE, etREAL, {&bfit},
+ "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
+ { "-efit", FALSE, etREAL, {&efit},
+ "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
+ { "-bvit", FALSE, etREAL, {&bvit},
+ "Begin of the fit of the current autocorrelation function to a*t^b."},
+ { "-evit", FALSE, etREAL, {&evit},
+ "End of the fit of the current autocorrelation function to a*t^b."},
+ { "-tr", FALSE, etREAL, {&trust},
+ "Fraction of the trajectory taken into account for the integral."},
+ { "-temp", FALSE, etREAL, {&temp},
+ "Temperature for calculating epsilon."}
+ };
+
+ output_env_t oenv;
+ t_topology top;
+ char title[STRLEN];
+ char **grpname = NULL;
+ const char *indexfn;
+ t_trxframe fr;
+ real *mass2 = NULL;
+ rvec *xtop, *vtop;
+ matrix box;
+ atom_id *index0;
+ int *indexm = NULL;
+ int isize;
+ t_trxstatus *status;
+ int flags = 0;
+ gmx_bool bTop;
+ gmx_bool bNEU;
+ gmx_bool bACF;
+ gmx_bool bINT;
+ int ePBC = -1;
+ int natoms;
+ int nmols;
+ int i, j, k = 0, l;
+ int step;
+ real t;
+ real lambda;
+ real *qmol;
+ FILE *outf = NULL;
+ FILE *outi = NULL;
+ FILE *tfile = NULL;
+ FILE *mcor = NULL;
+ FILE *fmj = NULL;
+ FILE *fmd = NULL;
+ FILE *fmjdsp = NULL;
+ FILE *fcur = NULL;
+ t_filenm fnm[] = {
+ { efTPS, NULL, NULL, ffREAD }, /* this is for the topology */
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
+ { efXVG, "-o", "current.xvg", ffWRITE },
+ { efXVG, "-caf", "caf.xvg", ffOPTWR },
+ { efXVG, "-dsp", "dsp.xvg", ffWRITE },
+ { efXVG, "-md", "md.xvg", ffWRITE },
+ { efXVG, "-mj", "mj.xvg", ffWRITE},
+ { efXVG, "-mc", "mc.xvg", ffOPTWR }
+ };
+
+#define NFILE asize(fnm)
+
+
+ const char *desc[] = {
+ "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
+ "of the rotational and translational dipole moment of the system, and the resulting static",
+ "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
+ "Furthermore, the routine is capable of extracting the static conductivity from the current ",
+ "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
+ "can be used to obtain the static conductivity."
+ "[PAR]",
+ "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
+ "correlation of the rotational and translational part of the dipole moment in the corresponding",
+ "file. However, this option is only available for trajectories containing velocities.",
+ "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
+ "autocorrelation functions. Since averaging proceeds by shifting the starting point",
+ "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
+ "starting points. Towards the end, statistical inaccuracy grows and integrating the",
+ "correlation function only yields reliable values until a certain point, depending on",
+ "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
+ "for calculating the static dielectric constant.",
+ "[PAR]",
+ "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
+ "[PAR]",
+ "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
+ "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
+ "tin-foil boundary conditions).",
+ "[PAR]",
+ "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
+ "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
+ "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
+ "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
+ "option has to be used with care, since only very short time spans fulfill the approximation that the density",
+ "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
+ "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
+ "the translational part of the dielectric constant."
+ };
+
+
+ /* At first the arguments will be parsed and the system information processed */
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+
+ bACF = opt2bSet("-caf", NFILE, fnm);
+ bINT = opt2bSet("-mc", NFILE, fnm);
+
+ bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
+
+
+
+ sfree(xtop);
+ sfree(vtop);
+ indexfn = ftp2fn_null(efNDX, NFILE, fnm);
+ snew(grpname, 1);
+
+ get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
+
+ flags = flags | TRX_READ_X | TRX_READ_V;
+
+ read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
+
+ snew(mass2, top.atoms.nr);
+ snew(qmol, top.atoms.nr);
+
+ bNEU = precalc(top, mass2, qmol);
+
+
+ snew(indexm, isize);
+
+ for (i = 0; i < isize; i++)
+ {
+ indexm[i] = index0[i];
+ }
+
+ nmols = isize;
+
+
+ index_atom2mol(&nmols, indexm, &top.mols);
+
+ if (fr.bV)
+ {
+ if (bACF)
+ {
+ outf = xvgropen(opt2fn("-caf", NFILE, fnm),
+ "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
+ "ACF (e nm/ps)\\S2", oenv);
+ fprintf(outf, "# time\t acf\t average \t std.dev\n");
+ }
+ fcur = xvgropen(opt2fn("-o", NFILE, fnm),
+ "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
+ fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
+ if (bINT)
+ {
+ mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
+ "M\\sD\\N - current autocorrelation function",
+ output_env_get_xvgr_tlabel(oenv),
+ "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
+ fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
+ }
+ }
+
+ fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
+ "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
+ "< M\\sJ\\N > (enm)", oenv);
+ fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
+ fmd = xvgropen(opt2fn("-md", NFILE, fnm),
+ "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
+ "< M\\sD\\N > (enm)", oenv);
+ fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
+
+ fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
+ "MSD of the squared translational dipole moment M",
+ output_env_get_xvgr_tlabel(oenv),
+ "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
+ oenv);
+
+
+ /* System information is read and prepared, dielectric() processes the frames
+ * and calculates the requested quantities */
+
+ dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
+ temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
+ index0, indexm, mass2, qmol, eps_rf, oenv);
+
+ gmx_ffclose(fmj);
+ gmx_ffclose(fmd);
+ gmx_ffclose(fmjdsp);
+ if (bACF)
+ {
+ gmx_ffclose(outf);
+ }
+ if (bINT)
+ {
+ gmx_ffclose(mcor);
+ }
+
+ return 0;
+}
--- /dev/null
- if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+
+#include "sysstuff.h"
+#include <string.h>
+#include "typedefs.h"
+#include "gromacs/utility/smalloc.h"
+#include "macros.h"
+#include "vec.h"
+#include "xvgr.h"
+#include "pbc.h"
+#include "gromacs/fileio/futil.h"
+#include "gromacs/commandline/pargs.h"
+#include "index.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trxio.h"
+#include "rmpbc.h"
+#include "gmx_ana.h"
+
+
+static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
+ real *rmin, real *rmax, int *min_ind)
+{
+#define NSHIFT 26
+ int sx, sy, sz, i, j, s;
+ real sqr_box, r2min, r2max, r2;
+ rvec shift[NSHIFT], d0, d;
+
+ sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));
+
+ s = 0;
+ for (sz = -1; sz <= 1; sz++)
+ {
+ for (sy = -1; sy <= 1; sy++)
+ {
+ for (sx = -1; sx <= 1; sx++)
+ {
+ if (sx != 0 || sy != 0 || sz != 0)
+ {
+ for (i = 0; i < DIM; i++)
+ {
+ shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
+ }
+ s++;
+ }
+ }
+ }
+ }
+
+ r2min = sqr_box;
+ r2max = 0;
+
+ for (i = 0; i < n; i++)
+ {
+ for (j = i+1; j < n; j++)
+ {
+ rvec_sub(x[index[i]], x[index[j]], d0);
+ r2 = norm2(d0);
+ if (r2 > r2max)
+ {
+ r2max = r2;
+ }
+ for (s = 0; s < NSHIFT; s++)
+ {
+ rvec_add(d0, shift[s], d);
+ r2 = norm2(d);
+ if (r2 < r2min)
+ {
+ r2min = r2;
+ min_ind[0] = i;
+ min_ind[1] = j;
+ }
+ }
+ }
+ }
+
+ *rmin = sqrt(r2min);
+ *rmax = sqrt(r2max);
+}
+
+static void periodic_mindist_plot(const char *trxfn, const char *outfn,
+ t_topology *top, int ePBC,
+ int n, atom_id index[], gmx_bool bSplit,
+ const output_env_t oenv)
+{
+ FILE *out;
+ const char *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" };
+ t_trxstatus *status;
+ real t;
+ rvec *x;
+ matrix box;
+ int natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
+ real r, rmin, rmax, rmint, tmint;
+ gmx_bool bFirst;
+ gmx_rmpbc_t gpbc = NULL;
+
+ natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
+
+ check_index(NULL, n, index, NULL, natoms);
+
+ out = xvgropen(outfn, "Minimum distance to periodic image",
+ output_env_get_time_label(oenv), "Distance (nm)", oenv);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(out, "@ subtitle \"and maximum internal distance\"\n");
+ }
+ xvgr_legend(out, 5, leg, oenv);
+
+ rmint = box[XX][XX];
+ tmint = 0;
+
+ if (NULL != top)
+ {
+ gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
+ }
+
+ bFirst = TRUE;
+ do
+ {
+ if (NULL != top)
+ {
+ gmx_rmpbc(gpbc, natoms, box, x);
+ }
+
+ periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
+ if (rmin < rmint)
+ {
+ rmint = rmin;
+ tmint = t;
+ ind_mini = ind_min[0];
+ ind_minj = ind_min[1];
+ }
- if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
++ if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
+ {
+ fprintf(out, "&\n");
+ }
+ fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
+ output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
+ bFirst = FALSE;
+ }
+ while (read_next_x(oenv, status, &t, x, box));
+
+ if (NULL != top)
+ {
+ gmx_rmpbc_done(gpbc);
+ }
+
+ gmx_ffclose(out);
+
+ fprintf(stdout,
+ "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
+ "between atoms %d and %d\n",
+ rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv),
+ index[ind_mini]+1, index[ind_minj]+1);
+}
+
+static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
+ int nx1, int nx2, atom_id index1[], atom_id index2[],
+ gmx_bool bGroup,
+ real *rmin, real *rmax, int *nmin, int *nmax,
+ int *ixmin, int *jxmin, int *ixmax, int *jxmax)
+{
+ int i, j, i0 = 0, j1;
+ int ix, jx;
+ atom_id *index3;
+ rvec dx;
+ real r2, rmin2, rmax2, rcut2;
+ t_pbc pbc;
+ int nmin_j, nmax_j;
+
+ *ixmin = -1;
+ *jxmin = -1;
+ *ixmax = -1;
+ *jxmax = -1;
+ *nmin = 0;
+ *nmax = 0;
+
+ rcut2 = sqr(rcut);
+
+ /* Must init pbc every step because of pressure coupling */
+ if (bPBC)
+ {
+ set_pbc(&pbc, ePBC, box);
+ }
+ if (index2)
+ {
+ i0 = 0;
+ j1 = nx2;
+ index3 = index2;
+ }
+ else
+ {
+ j1 = nx1;
+ index3 = index1;
+ }
+
+ rmin2 = 1e12;
+ rmax2 = -1e12;
+
+ for (j = 0; (j < j1); j++)
+ {
+ jx = index3[j];
+ if (index2 == NULL)
+ {
+ i0 = j + 1;
+ }
+ nmin_j = 0;
+ nmax_j = 0;
+ for (i = i0; (i < nx1); i++)
+ {
+ ix = index1[i];
+ if (ix != jx)
+ {
+ if (bPBC)
+ {
+ pbc_dx(&pbc, x[ix], x[jx], dx);
+ }
+ else
+ {
+ rvec_sub(x[ix], x[jx], dx);
+ }
+ r2 = iprod(dx, dx);
+ if (r2 < rmin2)
+ {
+ rmin2 = r2;
+ *ixmin = ix;
+ *jxmin = jx;
+ }
+ if (r2 > rmax2)
+ {
+ rmax2 = r2;
+ *ixmax = ix;
+ *jxmax = jx;
+ }
+ if (r2 <= rcut2)
+ {
+ nmin_j++;
+ }
+ else if (r2 > rcut2)
+ {
+ nmax_j++;
+ }
+ }
+ }
+ if (bGroup)
+ {
+ if (nmin_j > 0)
+ {
+ (*nmin)++;
+ }
+ if (nmax_j > 0)
+ {
+ (*nmax)++;
+ }
+ }
+ else
+ {
+ *nmin += nmin_j;
+ *nmax += nmax_j;
+ }
+ }
+ *rmin = sqrt(rmin2);
+ *rmax = sqrt(rmax2);
+}
+
+void dist_plot(const char *fn, const char *afile, const char *dfile,
+ const char *nfile, const char *rfile, const char *xfile,
+ real rcut, gmx_bool bMat, t_atoms *atoms,
+ int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit,
+ gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC,
+ gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName,
+ const output_env_t oenv)
+{
+ FILE *atm, *dist, *num;
+ t_trxstatus *trxout;
+ char buf[256];
+ char **leg;
+ real t, dmin, dmax, **mindres = NULL, **maxdres = NULL;
+ int nmin, nmax;
+ t_trxstatus *status;
+ int i = -1, j, k, natoms;
+ int min1, min2, max1, max2, min1r, min2r, max1r, max2r;
+ atom_id oindex[2];
+ rvec *x0;
+ matrix box;
+ t_trxframe frout;
+ gmx_bool bFirst;
+ FILE *respertime = NULL;
+
+ if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
+ {
+ gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
+ }
+
+ sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
+ dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
+ sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
+ num = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL;
+ atm = afile ? gmx_ffopen(afile, "w") : NULL;
+ trxout = xfile ? open_trx(xfile, "w") : NULL;
+
+ if (bMat)
+ {
+ if (ng == 1)
+ {
+ snew(leg, 1);
+ sprintf(buf, "Internal in %s", grpn[0]);
+ leg[0] = strdup(buf);
+ xvgr_legend(dist, 0, (const char**)leg, oenv);
+ if (num)
+ {
+ xvgr_legend(num, 0, (const char**)leg, oenv);
+ }
+ }
+ else
+ {
+ snew(leg, (ng*(ng-1))/2);
+ for (i = j = 0; (i < ng-1); i++)
+ {
+ for (k = i+1; (k < ng); k++, j++)
+ {
+ sprintf(buf, "%s-%s", grpn[i], grpn[k]);
+ leg[j] = strdup(buf);
+ }
+ }
+ xvgr_legend(dist, j, (const char**)leg, oenv);
+ if (num)
+ {
+ xvgr_legend(num, j, (const char**)leg, oenv);
+ }
+ }
+ }
+ else
+ {
+ snew(leg, ng-1);
+ for (i = 0; (i < ng-1); i++)
+ {
+ sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
+ leg[i] = strdup(buf);
+ }
+ xvgr_legend(dist, ng-1, (const char**)leg, oenv);
+ if (num)
+ {
+ xvgr_legend(num, ng-1, (const char**)leg, oenv);
+ }
+ }
+
+ if (bEachResEachTime)
+ {
+ sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
+ respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
+ xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
+ if (bPrintResName)
+ {
+ fprintf(respertime, "# ");
+ }
+ for (j = 0; j < nres; j++)
+ {
+ fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
+ }
+ fprintf(respertime, "\n");
+
+ }
+
+ j = 0;
+ if (nres)
+ {
+ snew(mindres, ng-1);
+ snew(maxdres, ng-1);
+ for (i = 1; i < ng; i++)
+ {
+ snew(mindres[i-1], nres);
+ snew(maxdres[i-1], nres);
+ for (j = 0; j < nres; j++)
+ {
+ mindres[i-1][j] = 1e6;
+ }
+ /* maxdres[*][*] is already 0 */
+ }
+ }
+ bFirst = TRUE;
+ do
+ {
++ if (bSplit && !bFirst && fabs(t/output_env_get_time_factor(oenv)) < 1e-5)
+ {
+ fprintf(dist, "&\n");
+ if (num)
+ {
+ fprintf(num, "&\n");
+ }
+ if (atm)
+ {
+ fprintf(atm, "&\n");
+ }
+ }
+ fprintf(dist, "%12e", output_env_conv_time(oenv, t));
+ if (num)
+ {
+ fprintf(num, "%12e", output_env_conv_time(oenv, t));
+ }
+
+ if (bMat)
+ {
+ if (ng == 1)
+ {
+ calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup,
+ &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
+ fprintf(dist, " %12e", bMin ? dmin : dmax);
+ if (num)
+ {
+ fprintf(num, " %8d", bMin ? nmin : nmax);
+ }
+ }
+ else
+ {
+ for (i = 0; (i < ng-1); i++)
+ {
+ for (k = i+1; (k < ng); k++)
+ {
+ calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k],
+ bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
+ fprintf(dist, " %12e", bMin ? dmin : dmax);
+ if (num)
+ {
+ fprintf(num, " %8d", bMin ? nmin : nmax);
+ }
+ }
+ }
+ }
+ }
+ else
+ {
+ for (i = 1; (i < ng); i++)
+ {
+ calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup,
+ &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
+ fprintf(dist, " %12e", bMin ? dmin : dmax);
+ if (num)
+ {
+ fprintf(num, " %8d", bMin ? nmin : nmax);
+ }
+ if (nres)
+ {
+ for (j = 0; j < nres; j++)
+ {
+ calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i],
+ &(index[0][residue[j]]), index[i], bGroup,
+ &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r);
+ mindres[i-1][j] = min(mindres[i-1][j], dmin);
+ maxdres[i-1][j] = max(maxdres[i-1][j], dmax);
+ }
+ }
+ }
+ }
+ fprintf(dist, "\n");
+ if (num)
+ {
+ fprintf(num, "\n");
+ }
+ if ( (bMin ? min1 : max1) != -1)
+ {
+ if (atm)
+ {
+ fprintf(atm, "%12e %12d %12d\n",
+ output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1),
+ 1+(bMin ? min2 : max2));
+ }
+ }
+
+ if (trxout)
+ {
+ oindex[0] = bMin ? min1 : max1;
+ oindex[1] = bMin ? min2 : max2;
+ write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL);
+ }
+ bFirst = FALSE;
+ /*dmin should be minimum distance for residue and group*/
+ if (bEachResEachTime)
+ {
+ fprintf(respertime, "%12e", t);
+ for (i = 1; i < ng; i++)
+ {
+ for (j = 0; j < nres; j++)
+ {
+ fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
+ /*reset distances for next time point*/
+ mindres[i-1][j] = 1e6;
+ maxdres[i-1][j] = 0;
+ }
+ }
+ fprintf(respertime, "\n");
+ }
+ }
+ while (read_next_x(oenv, status, &t, x0, box));
+
+ close_trj(status);
+ gmx_ffclose(dist);
+ if (num)
+ {
+ gmx_ffclose(num);
+ }
+ if (atm)
+ {
+ gmx_ffclose(atm);
+ }
+ if (trxout)
+ {
+ close_trx(trxout);
+ }
+
+ if (nres && !bEachResEachTime)
+ {
+ FILE *res;
+
+ sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
+ res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
+ xvgr_legend(res, ng-1, (const char**)leg, oenv);
+ for (j = 0; j < nres; j++)
+ {
+ fprintf(res, "%4d", j+1);
+ for (i = 1; i < ng; i++)
+ {
+ fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
+ }
+ fprintf(res, "\n");
+ }
+ }
+
+ sfree(x0);
+}
+
+int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
+{
+ int i;
+ int nres = 0, resnr, presnr;
+ int *residx;
+
+ /* build index of first atom numbers for each residue */
+ presnr = NOTSET;
+ snew(residx, atoms->nres+1);
+ for (i = 0; i < n; i++)
+ {
+ resnr = atoms->atom[index[i]].resind;
+ if (resnr != presnr)
+ {
+ residx[nres] = i;
+ nres++;
+ presnr = resnr;
+ }
+ }
+ if (debug)
+ {
+ printf("Found %d residues out of %d (%d/%d atoms)\n",
+ nres, atoms->nres, atoms->nr, n);
+ }
+ srenew(residx, nres+1);
+ /* mark end of last residue */
+ residx[nres] = n;
+ *resindex = residx;
+ return nres;
+}
+
+void dump_res(FILE *out, int nres, atom_id *resindex, atom_id index[])
+{
+ int i, j;
+
+ for (i = 0; i < nres-1; i++)
+ {
+ fprintf(out, "Res %d (%d):", i, resindex[i+1]-resindex[i]);
+ for (j = resindex[i]; j < resindex[i+1]; j++)
+ {
+ fprintf(out, " %d(%d)", j, index[j]);
+ }
+ fprintf(out, "\n");
+ }
+}
+
+int gmx_mindist(int argc, char *argv[])
+{
+ const char *desc[] = {
+ "[THISMODULE] computes the distance between one group and a number of",
+ "other groups. Both the minimum distance",
+ "(between any pair of atoms from the respective groups)",
+ "and the number of contacts within a given",
+ "distance are written to two separate output files.",
+ "With the [TT]-group[tt] option a contact of an atom in another group",
+ "with multiple atoms in the first group is counted as one contact",
+ "instead of as multiple contacts.",
+ "With [TT]-or[tt], minimum distances to each residue in the first",
+ "group are determined and plotted as a function of residue number.[PAR]",
+ "With option [TT]-pi[tt] the minimum distance of a group to its",
+ "periodic image is plotted. This is useful for checking if a protein",
+ "has seen its periodic image during a simulation. Only one shift in",
+ "each direction is considered, giving a total of 26 shifts.",
+ "It also plots the maximum distance within the group and the lengths",
+ "of the three box vectors.[PAR]",
+ "Also [gmx-distance] calculates distances."
+ };
+ const char *bugs[] = {
+ "The [TT]-pi[tt] option is very slow."
+ };
+
+ static gmx_bool bMat = FALSE, bPI = FALSE, bSplit = FALSE, bMax = FALSE, bPBC = TRUE;
+ static gmx_bool bGroup = FALSE;
+ static real rcutoff = 0.6;
+ static int ng = 1;
+ static gmx_bool bEachResEachTime = FALSE, bPrintResName = FALSE;
+ t_pargs pa[] = {
+ { "-matrix", FALSE, etBOOL, {&bMat},
+ "Calculate half a matrix of group-group distances" },
+ { "-max", FALSE, etBOOL, {&bMax},
+ "Calculate *maximum* distance instead of minimum" },
+ { "-d", FALSE, etREAL, {&rcutoff},
+ "Distance for contacts" },
+ { "-group", FALSE, etBOOL, {&bGroup},
+ "Count contacts with multiple atoms in the first group as one" },
+ { "-pi", FALSE, etBOOL, {&bPI},
+ "Calculate minimum distance with periodic images" },
+ { "-split", FALSE, etBOOL, {&bSplit},
+ "Split graph where time is zero" },
+ { "-ng", FALSE, etINT, {&ng},
+ "Number of secondary groups to compute distance to a central group" },
+ { "-pbc", FALSE, etBOOL, {&bPBC},
+ "Take periodic boundary conditions into account" },
+ { "-respertime", FALSE, etBOOL, {&bEachResEachTime},
+ "When writing per-residue distances, write distance for each time point" },
+ { "-printresname", FALSE, etBOOL, {&bPrintResName},
+ "Write residue names" }
+ };
+ output_env_t oenv;
+ t_topology *top = NULL;
+ int ePBC = -1;
+ char title[256];
+ real t;
+ rvec *x;
+ matrix box;
+ gmx_bool bTop = FALSE;
+
+ FILE *atm;
+ int i, j, nres = 0;
+ const char *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
+ char **grpname;
+ int *gnx;
+ atom_id **index, *residues = NULL;
+ t_filenm fnm[] = {
+ { efTRX, "-f", NULL, ffREAD },
+ { efTPS, NULL, NULL, ffOPTRD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efXVG, "-od", "mindist", ffWRITE },
+ { efXVG, "-on", "numcont", ffOPTWR },
+ { efOUT, "-o", "atm-pair", ffOPTWR },
+ { efTRO, "-ox", "mindist", ffOPTWR },
+ { efXVG, "-or", "mindistres", ffOPTWR }
+ };
+#define NFILE asize(fnm)
+
+ if (!parse_common_args(&argc, argv,
+ PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+
+ trxfnm = ftp2fn(efTRX, NFILE, fnm);
+ ndxfnm = ftp2fn_null(efNDX, NFILE, fnm);
+ distfnm = opt2fn("-od", NFILE, fnm);
+ numfnm = opt2fn_null("-on", NFILE, fnm);
+ atmfnm = ftp2fn_null(efOUT, NFILE, fnm);
+ oxfnm = opt2fn_null("-ox", NFILE, fnm);
+ resfnm = opt2fn_null("-or", NFILE, fnm);
+ if (bPI || resfnm != NULL)
+ {
+ /* We need a tps file */
+ tpsfnm = ftp2fn(efTPS, NFILE, fnm);
+ }
+ else
+ {
+ tpsfnm = ftp2fn_null(efTPS, NFILE, fnm);
+ }
+
+ if (!tpsfnm && !ndxfnm)
+ {
+ gmx_fatal(FARGS, "You have to specify either the index file or a tpr file");
+ }
+
+ if (bPI)
+ {
+ ng = 1;
+ fprintf(stderr, "Choose a group for distance calculation\n");
+ }
+ else if (!bMat)
+ {
+ ng++;
+ }
+
+ snew(gnx, ng);
+ snew(index, ng);
+ snew(grpname, ng);
+
+ if (tpsfnm || resfnm || !ndxfnm)
+ {
+ snew(top, 1);
+ bTop = read_tps_conf(tpsfnm, title, top, &ePBC, &x, NULL, box, FALSE);
+ if (bPI && !bTop)
+ {
+ printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
+ }
+ }
+ get_index(top ? &(top->atoms) : NULL, ndxfnm, ng, gnx, index, grpname);
+
+ if (bMat && (ng == 1))
+ {
+ ng = gnx[0];
+ printf("Special case: making distance matrix between all atoms in group %s\n",
+ grpname[0]);
+ srenew(gnx, ng);
+ srenew(index, ng);
+ srenew(grpname, ng);
+ for (i = 1; (i < ng); i++)
+ {
+ gnx[i] = 1;
+ grpname[i] = grpname[0];
+ snew(index[i], 1);
+ index[i][0] = index[0][i];
+ }
+ gnx[0] = 1;
+ }
+
+ if (resfnm)
+ {
+ nres = find_residues(top ? &(top->atoms) : NULL,
+ gnx[0], index[0], &residues);
+ if (debug)
+ {
+ dump_res(debug, nres, residues, index[0]);
+ }
+ }
+
+ if (bPI)
+ {
+ periodic_mindist_plot(trxfnm, distfnm, top, ePBC, gnx[0], index[0], bSplit, oenv);
+ }
+ else
+ {
+ dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
+ rcutoff, bMat, top ? &(top->atoms) : NULL,
+ ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
+ bGroup, bEachResEachTime, bPrintResName, oenv);
+ }
+
+ do_view(oenv, distfnm, "-nxy");
+ if (!bPI)
+ {
+ do_view(oenv, numfnm, "-nxy");
+ }
+
+ return 0;
+}
--- /dev/null
- abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "gromacs/utility/smalloc.h"
+#include <math.h>
+#include "macros.h"
+#include "typedefs.h"
+#include "xvgr.h"
+#include "copyrite.h"
+#include "gromacs/commandline/pargs.h"
+#include "vec.h"
+#include "index.h"
+#include "gmx_fatal.h"
+#include "gromacs/fileio/futil.h"
+#include "princ.h"
+#include "rmpbc.h"
+#include "gromacs/fileio/matio.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trxio.h"
+#include "cmat.h"
+#include "viewit.h"
+#include "gmx_ana.h"
+
+#include "gromacs/math/do_fit.h"
+
+static void norm_princ(t_atoms *atoms, int isize, atom_id *index, int natoms,
+ rvec *x)
+{
+ int i, m;
+ rvec princ, vec;
+
+ /* equalize principal components: */
+ /* orient principal axes, get principal components */
+ orient_princ(atoms, isize, index, natoms, x, NULL, princ);
+
+ /* calc our own principal components */
+ clear_rvec(vec);
+ for (m = 0; m < DIM; m++)
+ {
+ for (i = 0; i < isize; i++)
+ {
+ vec[m] += sqr(x[index[i]][m]);
+ }
+ vec[m] = sqrt(vec[m] / isize);
+ /* calculate scaling constants */
+ vec[m] = 1 / (sqrt(3) * vec[m]);
+ }
+
+ /* scale coordinates */
+ for (i = 0; i < natoms; i++)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ x[i][m] *= vec[m];
+ }
+ }
+}
+
+int gmx_rms(int argc, char *argv[])
+{
+ const char *desc[] =
+ {
+ "[THISMODULE] compares two structures by computing the root mean square",
+ "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
+ "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
+ "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
+ "This is selected by [TT]-what[tt].[PAR]"
+
+ "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
+ "reference structure. The reference structure",
+ "is taken from the structure file ([TT]-s[tt]).[PAR]",
+
+ "With option [TT]-mir[tt] also a comparison with the mirror image of",
+ "the reference structure is calculated.",
+ "This is useful as a reference for 'significant' values, see",
+ "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
+
+ "Option [TT]-prev[tt] produces the comparison with a previous frame",
+ "the specified number of frames ago.[PAR]",
+
+ "Option [TT]-m[tt] produces a matrix in [TT].xpm[tt] format of",
+ "comparison values of each structure in the trajectory with respect to",
+ "each other structure. This file can be visualized with for instance",
+ "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]",
+
+ "Option [TT]-fit[tt] controls the least-squares fitting of",
+ "the structures on top of each other: complete fit (rotation and",
+ "translation), translation only, or no fitting at all.[PAR]",
+
+ "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
+ "If you select the option (default) and ",
+ "supply a valid [TT].tpr[tt] file masses will be taken from there, ",
+ "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
+ "[TT]GMXLIB[tt]. This is fine for proteins, but not",
+ "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
+ "is assigned to unknown atoms. You can check whether this happend by",
+ "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
+
+ "With [TT]-f2[tt], the 'other structures' are taken from a second",
+ "trajectory, this generates a comparison matrix of one trajectory",
+ "versus the other.[PAR]",
+
+ "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
+
+ "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
+ "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
+ "comparison group are considered."
+ };
+ static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
+ static gmx_bool bDeltaLog = FALSE;
+ static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
+ static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1,
+ bond_user_min = -1, delta_maxy = 0.0;
+ /* strings and things for selecting difference method */
+ enum
+ {
+ ewSel, ewRMSD, ewRho, ewRhoSc, ewNR
+ };
+ int ewhat;
+ const char *what[ewNR + 1] =
+ { NULL, "rmsd", "rho", "rhosc", NULL };
+ const char *whatname[ewNR] =
+ { NULL, "RMSD", "Rho", "Rho sc" };
+ const char *whatlabel[ewNR] =
+ { NULL, "RMSD (nm)", "Rho", "Rho sc" };
+ const char *whatxvgname[ewNR] =
+ { NULL, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
+ const char *whatxvglabel[ewNR] =
+ { NULL, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
+ /* strings and things for fitting methods */
+ enum
+ {
+ efSel, efFit, efReset, efNone, efNR
+ };
+ int efit;
+ const char *fit[efNR + 1] =
+ { NULL, "rot+trans", "translation", "none", NULL };
+ const char *fitgraphlabel[efNR + 1] =
+ { NULL, "lsq fit", "translational fit", "no fit" };
+ static int nrms = 1;
+ static gmx_bool bMassWeighted = TRUE;
+ t_pargs pa[] =
+ {
+ { "-what", FALSE, etENUM,
+ { what }, "Structural difference measure" },
+ { "-pbc", FALSE, etBOOL,
+ { &bPBC }, "PBC check" },
+ { "-fit", FALSE, etENUM,
+ { fit }, "Fit to reference structure" },
+ { "-prev", FALSE, etINT,
+ { &prev }, "Compare with previous frame" },
+ { "-split", FALSE, etBOOL,
+ { &bSplit }, "Split graph where time is zero" },
+ { "-fitall", FALSE, etBOOL,
+ { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
+ { "-skip", FALSE, etINT,
+ { &freq }, "Only write every nr-th frame to matrix" },
+ { "-skip2", FALSE, etINT,
+ { &freq2 }, "Only write every nr-th frame to matrix" },
+ { "-max", FALSE, etREAL,
+ { &rmsd_user_max }, "Maximum level in comparison matrix" },
+ { "-min", FALSE, etREAL,
+ { &rmsd_user_min }, "Minimum level in comparison matrix" },
+ { "-bmax", FALSE, etREAL,
+ { &bond_user_max }, "Maximum level in bond angle matrix" },
+ { "-bmin", FALSE, etREAL,
+ { &bond_user_min }, "Minimum level in bond angle matrix" },
+ { "-mw", FALSE, etBOOL,
+ { &bMassWeighted }, "Use mass weighting for superposition" },
+ { "-nlevels", FALSE, etINT,
+ { &nlevels }, "Number of levels in the matrices" },
+ { "-ng", FALSE, etINT,
+ { &nrms }, "Number of groups to compute RMS between" },
+ { "-dlog", FALSE, etBOOL,
+ { &bDeltaLog },
+ "HIDDENUse a log x-axis in the delta t matrix" },
+ { "-dmax", FALSE, etREAL,
+ { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
+ { "-aver", FALSE, etINT,
+ { &avl },
+ "HIDDENAverage over this distance in the RMSD matrix" }
+ };
+ int natoms_trx, natoms_trx2, natoms;
+ int i, j, k, m, teller, teller2, tel_mat, tel_mat2;
+#define NFRAME 5000
+ int maxframe = NFRAME, maxframe2 = NFRAME;
+ real t, *w_rls, *w_rms, *w_rls_m = NULL, *w_rms_m = NULL;
+ gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
+ gmx_bool bFit, bReset;
+ t_topology top;
+ int ePBC;
+ t_iatom *iatom = NULL;
+
+ matrix box;
+ rvec *x, *xp, *xm = NULL, **mat_x = NULL, **mat_x2, *mat_x2_j = NULL, vec1,
+ vec2;
+ t_trxstatus *status;
+ char buf[256], buf2[256];
+ int ncons = 0;
+ FILE *fp;
+ real rlstot = 0, **rls, **rlsm = NULL, *time, *time2, *rlsnorm = NULL,
+ **rmsd_mat = NULL, **bond_mat = NULL, *axis, *axis2, *del_xaxis,
+ *del_yaxis, rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
+ real **rmsdav_mat = NULL, av_tot, weight, weight_tot;
+ real **delta = NULL, delta_max, delta_scalex = 0, delta_scaley = 0,
+ *delta_tot;
+ int delta_xsize = 0, del_lev = 100, mx, my, abs_my;
+ gmx_bool bA1, bA2, bPrev, bTop, *bInMat = NULL;
+ int ifit, *irms, ibond = 0, *ind_bond1 = NULL, *ind_bond2 = NULL, n_ind_m =
+ 0;
+ atom_id *ind_fit, **ind_rms, *ind_m = NULL, *rev_ind_m = NULL, *ind_rms_m =
+ NULL;
+ char *gn_fit, **gn_rms;
+ t_rgb rlo, rhi;
+ output_env_t oenv;
+ gmx_rmpbc_t gpbc = NULL;
+
+ t_filenm fnm[] =
+ {
+ { efTPS, NULL, NULL, ffREAD },
+ { efTRX, "-f", NULL, ffREAD },
+ { efTRX, "-f2", NULL, ffOPTRD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efXVG, NULL, "rmsd", ffWRITE },
+ { efXVG, "-mir", "rmsdmir", ffOPTWR },
+ { efXVG, "-a", "avgrp", ffOPTWR },
+ { efXVG, "-dist", "rmsd-dist", ffOPTWR },
+ { efXPM, "-m", "rmsd", ffOPTWR },
+ { efDAT, "-bin", "rmsd", ffOPTWR },
+ { efXPM, "-bm", "bond", ffOPTWR }
+ };
+#define NFILE asize(fnm)
+
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW
+ | PCA_BE_NICE, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
+ &oenv))
+ {
+ return 0;
+ }
+ /* parse enumerated options: */
+ ewhat = nenum(what);
+ if (ewhat == ewRho || ewhat == ewRhoSc)
+ {
+ please_cite(stdout, "Maiorov95");
+ }
+ efit = nenum(fit);
+ bFit = efit == efFit;
+ bReset = efit == efReset;
+ if (bFit)
+ {
+ bReset = TRUE; /* for fit, reset *must* be set */
+ }
+ else
+ {
+ bFitAll = FALSE;
+ }
+
+ /* mark active cmdline options */
+ bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
+ bFile2 = opt2bSet("-f2", NFILE, fnm);
+ bMat = opt2bSet("-m", NFILE, fnm);
+ bBond = opt2bSet("-bm", NFILE, fnm);
+ bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
+ * your RMSD matrix (hidden option */
+ bNorm = opt2bSet("-a", NFILE, fnm);
+ bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
+ if (freq <= 0)
+ {
+ fprintf(stderr, "The number of frames to skip is <= 0. "
+ "Writing out all frames.\n\n");
+ freq = 1;
+ }
+ if (!bFreq2)
+ {
+ freq2 = freq;
+ }
+ else if (bFile2 && freq2 <= 0)
+ {
+ fprintf(stderr,
+ "The number of frames to skip in second trajectory is <= 0.\n"
+ " Writing out all frames.\n\n");
+ freq2 = 1;
+ }
+
+ bPrev = (prev > 0);
+ if (bPrev)
+ {
+ prev = abs(prev);
+ if (freq != 1)
+ {
+ fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
+ }
+ }
+
+ if (bFile2 && !bMat && !bBond)
+ {
+ fprintf(
+ stderr,
+ "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
+ " will not read from %s\n", opt2fn("-f2", NFILE,
+ fnm));
+ bFile2 = FALSE;
+ }
+
+ if (bDelta)
+ {
+ bMat = TRUE;
+ if (bFile2)
+ {
+ fprintf(stderr,
+ "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
+ " will not read from %s\n", opt2fn("-f2",
+ NFILE, fnm));
+ bFile2 = FALSE;
+ }
+ }
+
+ bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xp,
+ NULL, box, TRUE);
+ snew(w_rls, top.atoms.nr);
+ snew(w_rms, top.atoms.nr);
+
+ if (!bTop && bBond)
+ {
+ fprintf(stderr,
+ "WARNING: Need a run input file for bond angle matrix,\n"
+ " will not calculate bond angle matrix.\n");
+ bBond = FALSE;
+ }
+
+ if (bReset)
+ {
+ fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares"
+ : "translational");
+ get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit,
+ &ind_fit, &gn_fit);
+ }
+ else
+ {
+ ifit = 0;
+ }
+
+ if (bReset)
+ {
+ if (bFit && ifit < 3)
+ {
+ gmx_fatal(FARGS, "Need >= 3 points to fit!\n" );
+ }
+
+ bMass = FALSE;
+ for (i = 0; i < ifit; i++)
+ {
+ if (bMassWeighted)
+ {
+ w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
+ }
+ else
+ {
+ w_rls[ind_fit[i]] = 1;
+ }
+ bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
+ }
+ if (!bMass)
+ {
+ fprintf(stderr, "All masses in the fit group are 0, using masses of 1\n");
+ for (i = 0; i < ifit; i++)
+ {
+ w_rls[ind_fit[i]] = 1;
+ }
+ }
+ }
+
+ if (bMat || bBond)
+ {
+ nrms = 1;
+ }
+
+ snew(gn_rms, nrms);
+ snew(ind_rms, nrms);
+ snew(irms, nrms);
+
+ fprintf(stderr, "Select group%s for %s calculation\n",
+ (nrms > 1) ? "s" : "", whatname[ewhat]);
+ get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
+ nrms, irms, ind_rms, gn_rms);
+
+ if (bNorm)
+ {
+ snew(rlsnorm, irms[0]);
+ }
+ snew(rls, nrms);
+ for (j = 0; j < nrms; j++)
+ {
+ snew(rls[j], maxframe);
+ }
+ if (bMirror)
+ {
+ snew(rlsm, nrms);
+ for (j = 0; j < nrms; j++)
+ {
+ snew(rlsm[j], maxframe);
+ }
+ }
+ snew(time, maxframe);
+ for (j = 0; j < nrms; j++)
+ {
+ bMass = FALSE;
+ for (i = 0; i < irms[j]; i++)
+ {
+ if (bMassWeighted)
+ {
+ w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
+ }
+ else
+ {
+ w_rms[ind_rms[j][i]] = 1.0;
+ }
+ bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
+ }
+ if (!bMass)
+ {
+ fprintf(stderr, "All masses in group %d are 0, using masses of 1\n", j);
+ for (i = 0; i < irms[j]; i++)
+ {
+ w_rms[ind_rms[j][i]] = 1;
+ }
+ }
+ }
+ /* Prepare reference frame */
+ if (bPBC)
+ {
+ gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
+ gmx_rmpbc(gpbc, top.atoms.nr, box, xp);
+ }
+ if (bReset)
+ {
+ reset_x(ifit, ind_fit, top.atoms.nr, NULL, xp, w_rls);
+ }
+ if (bMirror)
+ {
+ /* generate reference structure mirror image: */
+ snew(xm, top.atoms.nr);
+ for (i = 0; i < top.atoms.nr; i++)
+ {
+ copy_rvec(xp[i], xm[i]);
+ xm[i][XX] = -xm[i][XX];
+ }
+ }
+ if (ewhat == ewRhoSc)
+ {
+ norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
+ }
+
+ /* read first frame */
+ natoms_trx = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
+ if (natoms_trx != top.atoms.nr)
+ {
+ fprintf(stderr,
+ "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
+ top.atoms.nr, natoms_trx);
+ }
+ natoms = min(top.atoms.nr, natoms_trx);
+ if (bMat || bBond || bPrev)
+ {
+ snew(mat_x, NFRAME);
+
+ if (bPrev)
+ {
+ /* With -prev we use all atoms for simplicity */
+ n_ind_m = natoms;
+ }
+ else
+ {
+ /* Check which atoms we need (fit/rms) */
+ snew(bInMat, natoms);
+ for (i = 0; i < ifit; i++)
+ {
+ bInMat[ind_fit[i]] = TRUE;
+ }
+ n_ind_m = ifit;
+ for (i = 0; i < irms[0]; i++)
+ {
+ if (!bInMat[ind_rms[0][i]])
+ {
+ bInMat[ind_rms[0][i]] = TRUE;
+ n_ind_m++;
+ }
+ }
+ }
+ /* Make an index of needed atoms */
+ snew(ind_m, n_ind_m);
+ snew(rev_ind_m, natoms);
+ j = 0;
+ for (i = 0; i < natoms; i++)
+ {
+ if (bPrev || bInMat[i])
+ {
+ ind_m[j] = i;
+ rev_ind_m[i] = j;
+ j++;
+ }
+ }
+ snew(w_rls_m, n_ind_m);
+ snew(ind_rms_m, irms[0]);
+ snew(w_rms_m, n_ind_m);
+ for (i = 0; i < ifit; i++)
+ {
+ w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
+ }
+ for (i = 0; i < irms[0]; i++)
+ {
+ ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
+ w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
+ }
+ sfree(bInMat);
+ }
+
+ if (bBond)
+ {
+ ncons = 0;
+ for (k = 0; k < F_NRE; k++)
+ {
+ if (IS_CHEMBOND(k))
+ {
+ iatom = top.idef.il[k].iatoms;
+ ncons += top.idef.il[k].nr/3;
+ }
+ }
+ fprintf(stderr, "Found %d bonds in topology\n", ncons);
+ snew(ind_bond1, ncons);
+ snew(ind_bond2, ncons);
+ ibond = 0;
+ for (k = 0; k < F_NRE; k++)
+ {
+ if (IS_CHEMBOND(k))
+ {
+ iatom = top.idef.il[k].iatoms;
+ ncons = top.idef.il[k].nr/3;
+ for (i = 0; i < ncons; i++)
+ {
+ bA1 = FALSE;
+ bA2 = FALSE;
+ for (j = 0; j < irms[0]; j++)
+ {
+ if (iatom[3*i+1] == ind_rms[0][j])
+ {
+ bA1 = TRUE;
+ }
+ if (iatom[3*i+2] == ind_rms[0][j])
+ {
+ bA2 = TRUE;
+ }
+ }
+ if (bA1 && bA2)
+ {
+ ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
+ ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
+ ibond++;
+ }
+ }
+ }
+ }
+ fprintf(stderr, "Using %d bonds for bond angle matrix\n", ibond);
+ if (ibond == 0)
+ {
+ gmx_fatal(FARGS, "0 bonds found");
+ }
+ }
+
+ /* start looping over frames: */
+ tel_mat = 0;
+ teller = 0;
+ do
+ {
+ if (bPBC)
+ {
+ gmx_rmpbc(gpbc, natoms, box, x);
+ }
+
+ if (bReset)
+ {
+ reset_x(ifit, ind_fit, natoms, NULL, x, w_rls);
+ }
+ if (ewhat == ewRhoSc)
+ {
+ norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
+ }
+
+ if (bFit)
+ {
+ /*do the least squares fit to original structure*/
+ do_fit(natoms, w_rls, xp, x);
+ }
+
+ if (teller % freq == 0)
+ {
+ /* keep frame for matrix calculation */
+ if (bMat || bBond || bPrev)
+ {
+ if (tel_mat >= NFRAME)
+ {
+ srenew(mat_x, tel_mat+1);
+ }
+ snew(mat_x[tel_mat], n_ind_m);
+ for (i = 0; i < n_ind_m; i++)
+ {
+ copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]);
+ }
+ }
+ tel_mat++;
+ }
+
+ /*calculate energy of root_least_squares*/
+ if (bPrev)
+ {
+ j = tel_mat-prev-1;
+ if (j < 0)
+ {
+ j = 0;
+ }
+ for (i = 0; i < n_ind_m; i++)
+ {
+ copy_rvec(mat_x[j][i], xp[ind_m[i]]);
+ }
+ if (bReset)
+ {
+ reset_x(ifit, ind_fit, natoms, NULL, xp, w_rls);
+ }
+ if (bFit)
+ {
+ do_fit(natoms, w_rls, x, xp);
+ }
+ }
+ for (j = 0; (j < nrms); j++)
+ {
+ rls[j][teller] =
+ calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
+ }
+ if (bNorm)
+ {
+ for (j = 0; (j < irms[0]); j++)
+ {
+ rlsnorm[j] +=
+ calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
+ }
+ }
+
+ if (bMirror)
+ {
+ if (bFit)
+ {
+ /*do the least squares fit to mirror of original structure*/
+ do_fit(natoms, w_rls, xm, x);
+ }
+
+ for (j = 0; j < nrms; j++)
+ {
+ rlsm[j][teller] =
+ calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
+ }
+ }
+ time[teller] = output_env_conv_time(oenv, t);
+
+ teller++;
+ if (teller >= maxframe)
+ {
+ maxframe += NFRAME;
+ srenew(time, maxframe);
+ for (j = 0; (j < nrms); j++)
+ {
+ srenew(rls[j], maxframe);
+ }
+ if (bMirror)
+ {
+ for (j = 0; (j < nrms); j++)
+ {
+ srenew(rlsm[j], maxframe);
+ }
+ }
+ }
+ }
+ while (read_next_x(oenv, status, &t, x, box));
+ close_trj(status);
+
+ if (bFile2)
+ {
+ snew(time2, maxframe2);
+
+ fprintf(stderr, "\nWill read second trajectory file\n");
+ snew(mat_x2, NFRAME);
+ natoms_trx2 =
+ read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
+ if (natoms_trx2 != natoms_trx)
+ {
+ gmx_fatal(FARGS,
+ "Second trajectory (%d atoms) does not match the first one"
+ " (%d atoms)", natoms_trx2, natoms_trx);
+ }
+ tel_mat2 = 0;
+ teller2 = 0;
+ do
+ {
+ if (bPBC)
+ {
+ gmx_rmpbc(gpbc, natoms, box, x);
+ }
+
+ if (bReset)
+ {
+ reset_x(ifit, ind_fit, natoms, NULL, x, w_rls);
+ }
+ if (ewhat == ewRhoSc)
+ {
+ norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
+ }
+
+ if (bFit)
+ {
+ /*do the least squares fit to original structure*/
+ do_fit(natoms, w_rls, xp, x);
+ }
+
+ if (teller2 % freq2 == 0)
+ {
+ /* keep frame for matrix calculation */
+ if (bMat)
+ {
+ if (tel_mat2 >= NFRAME)
+ {
+ srenew(mat_x2, tel_mat2+1);
+ }
+ snew(mat_x2[tel_mat2], n_ind_m);
+ for (i = 0; i < n_ind_m; i++)
+ {
+ copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]);
+ }
+ }
+ tel_mat2++;
+ }
+
+ time2[teller2] = output_env_conv_time(oenv, t);
+
+ teller2++;
+ if (teller2 >= maxframe2)
+ {
+ maxframe2 += NFRAME;
+ srenew(time2, maxframe2);
+ }
+ }
+ while (read_next_x(oenv, status, &t, x, box));
+ close_trj(status);
+ }
+ else
+ {
+ mat_x2 = mat_x;
+ time2 = time;
+ tel_mat2 = tel_mat;
+ freq2 = freq;
+ }
+ gmx_rmpbc_done(gpbc);
+
+ if (bMat || bBond)
+ {
+ /* calculate RMS matrix */
+ fprintf(stderr, "\n");
+ if (bMat)
+ {
+ fprintf(stderr, "Building %s matrix, %dx%d elements\n",
+ whatname[ewhat], tel_mat, tel_mat2);
+ snew(rmsd_mat, tel_mat);
+ }
+ if (bBond)
+ {
+ fprintf(stderr, "Building bond angle matrix, %dx%d elements\n",
+ tel_mat, tel_mat2);
+ snew(bond_mat, tel_mat);
+ }
+ snew(axis, tel_mat);
+ snew(axis2, tel_mat2);
+ rmsd_max = 0;
+ if (bFile2)
+ {
+ rmsd_min = 1e10;
+ }
+ else
+ {
+ rmsd_min = 0;
+ }
+ rmsd_avg = 0;
+ bond_max = 0;
+ bond_min = 1e10;
+ for (j = 0; j < tel_mat2; j++)
+ {
+ axis2[j] = time2[freq2*j];
+ }
+ if (bDelta)
+ {
+ if (bDeltaLog)
+ {
+ delta_scalex = 8.0/log(2.0);
+ delta_xsize = (int)(log(tel_mat/2)*delta_scalex+0.5)+1;
+ }
+ else
+ {
+ delta_xsize = tel_mat/2;
+ }
+ delta_scaley = 1.0/delta_maxy;
+ snew(delta, delta_xsize);
+ for (j = 0; j < delta_xsize; j++)
+ {
+ snew(delta[j], del_lev+1);
+ }
+ if (avl > 0)
+ {
+ snew(rmsdav_mat, tel_mat);
+ for (j = 0; j < tel_mat; j++)
+ {
+ snew(rmsdav_mat[j], tel_mat);
+ }
+ }
+ }
+
+ if (bFitAll)
+ {
+ snew(mat_x2_j, natoms);
+ }
+ for (i = 0; i < tel_mat; i++)
+ {
+ axis[i] = time[freq*i];
+ fprintf(stderr, "\r element %5d; time %5.2f ", i, axis[i]);
+ if (bMat)
+ {
+ snew(rmsd_mat[i], tel_mat2);
+ }
+ if (bBond)
+ {
+ snew(bond_mat[i], tel_mat2);
+ }
+ for (j = 0; j < tel_mat2; j++)
+ {
+ if (bFitAll)
+ {
+ for (k = 0; k < n_ind_m; k++)
+ {
+ copy_rvec(mat_x2[j][k], mat_x2_j[k]);
+ }
+ do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j);
+ }
+ else
+ {
+ mat_x2_j = mat_x2[j];
+ }
+ if (bMat)
+ {
+ if (bFile2 || (i < j))
+ {
+ rmsd_mat[i][j] =
+ calc_similar_ind(ewhat != ewRMSD, irms[0], ind_rms_m,
+ w_rms_m, mat_x[i], mat_x2_j);
+ if (rmsd_mat[i][j] > rmsd_max)
+ {
+ rmsd_max = rmsd_mat[i][j];
+ }
+ if (rmsd_mat[i][j] < rmsd_min)
+ {
+ rmsd_min = rmsd_mat[i][j];
+ }
+ rmsd_avg += rmsd_mat[i][j];
+ }
+ else
+ {
+ rmsd_mat[i][j] = rmsd_mat[j][i];
+ }
+ }
+ if (bBond)
+ {
+ if (bFile2 || (i <= j))
+ {
+ ang = 0.0;
+ for (m = 0; m < ibond; m++)
+ {
+ rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1);
+ rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2);
+ ang += acos(cos_angle(vec1, vec2));
+ }
+ bond_mat[i][j] = ang*180.0/(M_PI*ibond);
+ if (bond_mat[i][j] > bond_max)
+ {
+ bond_max = bond_mat[i][j];
+ }
+ if (bond_mat[i][j] < bond_min)
+ {
+ bond_min = bond_mat[i][j];
+ }
+ }
+ else
+ {
+ bond_mat[i][j] = bond_mat[j][i];
+ }
+ }
+ }
+ }
+ if (bFile2)
+ {
+ rmsd_avg /= tel_mat*tel_mat2;
+ }
+ else
+ {
+ rmsd_avg /= tel_mat*(tel_mat - 1)/2;
+ }
+ if (bMat && (avl > 0))
+ {
+ rmsd_max = 0.0;
+ rmsd_min = 0.0;
+ rmsd_avg = 0.0;
+ for (j = 0; j < tel_mat-1; j++)
+ {
+ for (i = j+1; i < tel_mat; i++)
+ {
+ av_tot = 0;
+ weight_tot = 0;
+ for (my = -avl; my <= avl; my++)
+ {
+ if ((j+my >= 0) && (j+my < tel_mat))
+ {
+ abs_my = abs(my);
+ for (mx = -avl; mx <= avl; mx++)
+ {
+ if ((i+mx >= 0) && (i+mx < tel_mat))
+ {
+ weight = (real)(avl+1-max(abs(mx), abs_my));
+ av_tot += weight*rmsd_mat[i+mx][j+my];
+ weight_tot += weight;
+ }
+ }
+ }
+ }
+ rmsdav_mat[i][j] = av_tot/weight_tot;
+ rmsdav_mat[j][i] = rmsdav_mat[i][j];
+ if (rmsdav_mat[i][j] > rmsd_max)
+ {
+ rmsd_max = rmsdav_mat[i][j];
+ }
+ }
+ }
+ rmsd_mat = rmsdav_mat;
+ }
+
+ if (bMat)
+ {
+ fprintf(stderr, "\n%s: Min %f, Max %f, Avg %f\n",
+ whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg);
+ rlo.r = 1; rlo.g = 1; rlo.b = 1;
+ rhi.r = 0; rhi.g = 0; rhi.b = 0;
+ if (rmsd_user_max != -1)
+ {
+ rmsd_max = rmsd_user_max;
+ }
+ if (rmsd_user_min != -1)
+ {
+ rmsd_min = rmsd_user_min;
+ }
+ if ((rmsd_user_max != -1) || (rmsd_user_min != -1))
+ {
+ fprintf(stderr, "Min and Max value set to resp. %f and %f\n",
+ rmsd_min, rmsd_max);
+ }
+ sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]);
+ write_xpm(opt2FILE("-m", NFILE, fnm, "w"), 0, buf, whatlabel[ewhat],
+ output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
+ axis, axis2, rmsd_mat, rmsd_min, rmsd_max, rlo, rhi, &nlevels);
+ /* Print the distribution of RMSD values */
+ if (opt2bSet("-dist", NFILE, fnm))
+ {
+ low_rmsd_dist(opt2fn("-dist", NFILE, fnm), rmsd_max, tel_mat, rmsd_mat, oenv);
+ }
+
+ if (bDelta)
+ {
+ snew(delta_tot, delta_xsize);
+ for (j = 0; j < tel_mat-1; j++)
+ {
+ for (i = j+1; i < tel_mat; i++)
+ {
+ mx = i-j;
+ if (mx < tel_mat/2)
+ {
+ if (bDeltaLog)
+ {
+ mx = (int)(log(mx)*delta_scalex+0.5);
+ }
+ my = (int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
+ delta_tot[mx] += 1.0;
+ if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy))
+ {
+ delta[mx][my] += 1.0;
+ }
+ }
+ }
+ }
+ delta_max = 0;
+ for (i = 0; i < delta_xsize; i++)
+ {
+ if (delta_tot[i] > 0.0)
+ {
+ delta_tot[i] = 1.0/delta_tot[i];
+ for (j = 0; j <= del_lev; j++)
+ {
+ delta[i][j] *= delta_tot[i];
+ if (delta[i][j] > delta_max)
+ {
+ delta_max = delta[i][j];
+ }
+ }
+ }
+ }
+ fprintf(stderr, "Maximum in delta matrix: %f\n", delta_max);
+ snew(del_xaxis, delta_xsize);
+ snew(del_yaxis, del_lev+1);
+ for (i = 0; i < delta_xsize; i++)
+ {
+ del_xaxis[i] = axis[i]-axis[0];
+ }
+ for (i = 0; i < del_lev+1; i++)
+ {
+ del_yaxis[i] = delta_maxy*i/del_lev;
+ }
+ sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]);
+ fp = gmx_ffopen("delta.xpm", "w");
+ write_xpm(fp, 0, buf, "density", output_env_get_time_label(oenv), whatlabel[ewhat],
+ delta_xsize, del_lev+1, del_xaxis, del_yaxis,
+ delta, 0.0, delta_max, rlo, rhi, &nlevels);
+ gmx_ffclose(fp);
+ }
+ if (opt2bSet("-bin", NFILE, fnm))
+ {
+ /* NB: File must be binary if we use fwrite */
+ fp = ftp2FILE(efDAT, NFILE, fnm, "wb");
+ for (i = 0; i < tel_mat; i++)
+ {
+ if (fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp) != tel_mat2)
+ {
+ gmx_fatal(FARGS, "Error writing to output file");
+ }
+ }
+ gmx_ffclose(fp);
+ }
+ }
+ if (bBond)
+ {
+ fprintf(stderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max);
+ if (bond_user_max != -1)
+ {
+ bond_max = bond_user_max;
+ }
+ if (bond_user_min != -1)
+ {
+ bond_min = bond_user_min;
+ }
+ if ((bond_user_max != -1) || (bond_user_min != -1))
+ {
+ fprintf(stderr, "Bond angle Min and Max set to:\n"
+ "Min. angle: %f, Max. angle: %f\n", bond_min, bond_max);
+ }
+ rlo.r = 1; rlo.g = 1; rlo.b = 1;
+ rhi.r = 0; rhi.g = 0; rhi.b = 0;
+ sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);
+ write_xpm(opt2FILE("-bm", NFILE, fnm, "w"), 0, buf, "degrees",
+ output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
+ axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels);
+ }
+ }
+
+ bAv = opt2bSet("-a", NFILE, fnm);
+
+ /* Write the RMSD's to file */
+ if (!bPrev)
+ {
+ sprintf(buf, "%s", whatxvgname[ewhat]);
+ }
+ else
+ {
+ sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat],
+ time[prev*freq]-time[0], output_env_get_time_label(oenv));
+ }
+ fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
+ whatxvglabel[ewhat], oenv);
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n",
+ (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit],
+ bFit ? " to " : "", bFit ? gn_fit : "");
+ }
+ if (nrms != 1)
+ {
+ xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
+ }
+ for (i = 0; (i < teller); i++)
+ {
+ if (bSplit && i > 0 &&
- if (bSplit && i > 0 && abs(time[i]) < 1e-5)
++ fabs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
+ {
+ fprintf(fp, "&\n");
+ }
+ fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);
+ for (j = 0; (j < nrms); j++)
+ {
+ fprintf(fp, " %12.7f", rls[j][i]);
+ if (bAv)
+ {
+ rlstot += rls[j][i];
+ }
+ }
+ fprintf(fp, "\n");
+ }
+ gmx_ffclose(fp);
+
+ if (bMirror)
+ {
+ /* Write the mirror RMSD's to file */
+ sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
+ sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
+ fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
+ buf2, oenv);
+ if (nrms == 1)
+ {
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
+ gn_rms[0], gn_fit);
+ }
+ }
+ else
+ {
+ if (output_env_get_print_xvgr_codes(oenv))
+ {
+ fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", gn_fit);
+ }
+ xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
+ }
+ for (i = 0; (i < teller); i++)
+ {
++ if (bSplit && i > 0 && fabs(time[i]) < 1e-5)
+ {
+ fprintf(fp, "&\n");
+ }
+ fprintf(fp, "%12.7f", time[i]);
+ for (j = 0; (j < nrms); j++)
+ {
+ fprintf(fp, " %12.7f", rlsm[j][i]);
+ }
+ fprintf(fp, "\n");
+ }
+ gmx_ffclose(fp);
+ }
+
+ if (bAv)
+ {
+ sprintf(buf, "Average %s", whatxvgname[ewhat]);
+ sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
+ fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
+ for (j = 0; (j < nrms); j++)
+ {
+ fprintf(fp, "%10d %10g\n", j, rlstot/teller);
+ }
+ gmx_ffclose(fp);
+ }
+
+ if (bNorm)
+ {
+ fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
+ for (j = 0; (j < irms[0]); j++)
+ {
+ fprintf(fp, "%10d %10g\n", j, rlsnorm[j]/teller);
+ }
+ gmx_ffclose(fp);
+ }
+ do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
+ do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn_null("-mir", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn_null("-m", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn_null("-bm", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn_null("-dist", NFILE, fnm), NULL);
+
+ return 0;
+}
--- /dev/null
- abs(mat[i].axis_x[x]) <
- 0.1*abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "typedefs.h"
+#include "macros.h"
+#include "gmx_fatal.h"
+#include "gromacs/utility/smalloc.h"
+#include "viewit.h"
+#include "gmx_ana.h"
+
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/futil.h"
+#include "gromacs/fileio/matio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/writeps.h"
+
+#define FUDGE 1.2
+#define DDD 2
+
+typedef struct {
+ real major;
+ real minor;
+ real offset;
+ gmx_bool first;
+ int lineatzero;
+ real majorticklen;
+ real minorticklen;
+ char label[STRLEN];
+ real fontsize;
+ char font[STRLEN];
+ real tickfontsize;
+ char tickfont[STRLEN];
+} t_axisdef;
+
+typedef struct {
+ int bw;
+ real linewidth;
+ real xoffs, yoffs;
+ gmx_bool bTitle;
+ gmx_bool bTitleOnce;
+ gmx_bool bYonce;
+ real titfontsize;
+ char titfont[STRLEN];
+ gmx_bool legend;
+ real legfontsize;
+ char legfont[STRLEN];
+ char leglabel[STRLEN];
+ char leg2label[STRLEN];
+ real xboxsize;
+ real yboxsize;
+ real boxspacing;
+ real boxlinewidth;
+ real ticklinewidth;
+ real zerolinewidth;
+ t_axisdef X, Y;
+} t_psrec;
+
+/* MUST correspond to char *legend[] in main() */
+enum {
+ elSel, elBoth, elFirst, elSecond, elNone, elNR
+};
+
+/* MUST correspond to char *combine[] in main() */
+enum {
+ ecSel, ecHalves, ecAdd, ecSub, ecMult, ecDiv, ecNR
+};
+
+void get_params(const char *mpin, const char *mpout, t_psrec *psr)
+{
+ static const char *gmx_bools[BOOL_NR+1] = { "no", "yes", NULL };
+ /* this must correspond to t_rgb *linecolors[] below */
+ static const char *colors[] = { "none", "black", "white", NULL };
+ warninp_t wi;
+ t_inpfile *inp;
+ const char *tmp;
+ int ninp = 0;
+
+ wi = init_warning(FALSE, 0);
+
+ if (mpin != NULL)
+ {
+ inp = read_inpfile(mpin, &ninp, wi);
+ }
+ else
+ {
+ inp = NULL;
+ }
+ ETYPE("black&white", psr->bw, gmx_bools);
+ RTYPE("linewidth", psr->linewidth, 1.0);
+ STYPE("titlefont", psr->titfont, "Helvetica");
+ RTYPE("titlefontsize", psr->titfontsize, 20.0);
+ ETYPE("legend", psr->legend, gmx_bools);
+ STYPE("legendfont", psr->legfont, psr->titfont);
+ STYPE("legendlabel", psr->leglabel, "");
+ STYPE("legend2label", psr->leg2label, psr->leglabel);
+ RTYPE("legendfontsize", psr->legfontsize, 14.0);
+ RTYPE("xbox", psr->xboxsize, 0.0);
+ RTYPE("ybox", psr->yboxsize, 0.0);
+ RTYPE("matrixspacing", psr->boxspacing, 20.0);
+ RTYPE("xoffset", psr->xoffs, 0.0);
+ RTYPE("yoffset", psr->yoffs, psr->xoffs);
+ RTYPE("boxlinewidth", psr->boxlinewidth, psr->linewidth);
+ RTYPE("ticklinewidth", psr->ticklinewidth, psr->linewidth);
+ RTYPE("zerolinewidth", psr->zerolinewidth, psr->ticklinewidth);
+ ETYPE("x-lineat0value", psr->X.lineatzero, colors);
+ RTYPE("x-major", psr->X.major, NOTSET);
+ RTYPE("x-minor", psr->X.minor, NOTSET);
+ RTYPE("x-firstmajor", psr->X.offset, 0.0);
+ ETYPE("x-majorat0", psr->X.first, gmx_bools);
+ RTYPE("x-majorticklen", psr->X.majorticklen, 8.0);
+ RTYPE("x-minorticklen", psr->X.minorticklen, 4.0);
+ STYPE("x-label", psr->X.label, "");
+ RTYPE("x-fontsize", psr->X.fontsize, 16.0);
+ STYPE("x-font", psr->X.font, psr->titfont);
+ RTYPE("x-tickfontsize", psr->X.tickfontsize, 10.0);
+ STYPE("x-tickfont", psr->X.tickfont, psr->X.font);
+ ETYPE("y-lineat0value", psr->Y.lineatzero, colors);
+ RTYPE("y-major", psr->Y.major, psr->X.major);
+ RTYPE("y-minor", psr->Y.minor, psr->X.minor);
+ RTYPE("y-firstmajor", psr->Y.offset, psr->X.offset);
+ ETYPE("y-majorat0", psr->Y.first, gmx_bools);
+ RTYPE("y-majorticklen", psr->Y.majorticklen, psr->X.majorticklen);
+ RTYPE("y-minorticklen", psr->Y.minorticklen, psr->X.minorticklen);
+ STYPE("y-label", psr->Y.label, psr->X.label);
+ RTYPE("y-fontsize", psr->Y.fontsize, psr->X.fontsize);
+ STYPE("y-font", psr->Y.font, psr->X.font);
+ RTYPE("y-tickfontsize", psr->Y.tickfontsize, psr->X.tickfontsize);
+ STYPE("y-tickfont", psr->Y.tickfont, psr->Y.font);
+
+ if (mpout != NULL)
+ {
+ write_inpfile(mpout, ninp, inp, TRUE, wi);
+ }
+
+ done_warning(wi, FARGS);
+}
+
+t_rgb black = { 0, 0, 0 };
+t_rgb white = { 1, 1, 1 };
+t_rgb red = { 1, 0, 0 };
+t_rgb blue = { 0, 0, 1 };
+#define BLACK (&black)
+/* this must correspond to *colors[] in get_params */
+t_rgb *linecolors[] = { NULL, &black, &white, NULL };
+
+gmx_bool diff_maps(int nmap1, t_mapping *map1, int nmap2, t_mapping *map2)
+{
+ int i;
+ gmx_bool bDiff, bColDiff = FALSE;
+
+ if (nmap1 != nmap2)
+ {
+ bDiff = TRUE;
+ }
+ else
+ {
+ bDiff = FALSE;
+ for (i = 0; i < nmap1; i++)
+ {
+ if (!matelmt_cmp(map1[i].code, map2[i].code))
+ {
+ bDiff = TRUE;
+ }
+ if (strcmp(map1[i].desc, map2[i].desc) != 0)
+ {
+ bDiff = TRUE;
+ }
+ if ((map1[i].rgb.r != map2[i].rgb.r) ||
+ (map1[i].rgb.g != map2[i].rgb.g) ||
+ (map1[i].rgb.b != map2[i].rgb.b))
+ {
+ bColDiff = TRUE;
+ }
+ }
+ if (!bDiff && bColDiff)
+ {
+ fprintf(stderr, "Warning: two colormaps differ only in RGB value, using one colormap.\n");
+ }
+ }
+
+ return bDiff;
+}
+
+void leg_discrete(t_psdata ps, real x0, real y0, char *label,
+ real fontsize, char *font, int nmap, t_mapping map[])
+{
+ int i;
+ real yhh;
+ real boxhh;
+
+ boxhh = fontsize+DDD;
+ /* LANDSCAPE */
+ ps_rgb(ps, BLACK);
+ ps_strfont(ps, font, fontsize);
+ yhh = y0+fontsize+3*DDD;
+ if ((int)strlen(label) > 0)
+ {
+ ps_ctext(ps, x0, yhh, label, eXLeft);
+ }
+ ps_moveto(ps, x0, y0);
+ for (i = 0; (i < nmap); i++)
+ {
+ ps_setorigin(ps);
+ ps_rgb(ps, &(map[i].rgb));
+ ps_fillbox(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
+ ps_rgb(ps, BLACK);
+ ps_box(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
+ ps_ctext(ps, boxhh+2*DDD, fontsize/3, map[i].desc, eXLeft);
+ ps_unsetorigin(ps);
+ ps_moverel(ps, DDD, -fontsize/3);
+ }
+}
+
+void leg_continuous(t_psdata ps, real x0, real x, real y0, char *label,
+ real fontsize, char *font,
+ int nmap, t_mapping map[],
+ int mapoffset)
+{
+ int i;
+ real xx0;
+ real yhh, boxxh, boxyh;
+
+ boxyh = fontsize;
+ if (x < 8*fontsize)
+ {
+ x = 8*fontsize;
+ }
+ boxxh = (real)x/(real)(nmap-mapoffset);
+ if (boxxh > fontsize)
+ {
+ boxxh = fontsize;
+ }
+
+ /* LANDSCAPE */
+ xx0 = x0-((nmap-mapoffset)*boxxh)/2.0;
+
+ for (i = 0; (i < nmap-mapoffset); i++)
+ {
+ ps_rgb(ps, &(map[i+mapoffset].rgb));
+ ps_fillbox(ps, xx0+i*boxxh, y0, xx0+(i+1)*boxxh, y0+boxyh);
+ }
+ ps_strfont(ps, font, fontsize);
+ ps_rgb(ps, BLACK);
+ ps_box(ps, xx0, y0, xx0+(nmap-mapoffset)*boxxh, y0+boxyh);
+
+ yhh = y0+boxyh+3*DDD;
+ ps_ctext(ps, xx0+boxxh/2, yhh, map[0].desc, eXCenter);
+ if ((int)strlen(label) > 0)
+ {
+ ps_ctext(ps, x0, yhh, label, eXCenter);
+ }
+ ps_ctext(ps, xx0+((nmap-mapoffset)*boxxh)
+ - boxxh/2, yhh, map[nmap-1].desc, eXCenter);
+}
+
+void leg_bicontinuous(t_psdata ps, real x0, real x, real y0, char *label1,
+ char *label2, real fontsize, char *font,
+ int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
+{
+ real xx1, xx2, x1, x2;
+
+ x1 = x/(nmap1+nmap2)*nmap1; /* width of legend 1 */
+ x2 = x/(nmap1+nmap2)*nmap2; /* width of legend 2 */
+ xx1 = x0-(x2/2.0)-fontsize; /* center of legend 1 */
+ xx2 = x0+(x1/2.0)+fontsize; /* center of legend 2 */
+ x1 -= fontsize/2; /* adjust width */
+ x2 -= fontsize/2; /* adjust width */
+ leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, nmap1, map1, 0);
+ leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, nmap2, map2, 0);
+}
+
+static real box_height(t_matrix *mat, t_psrec *psr)
+{
+ return mat->ny*psr->yboxsize;
+}
+
+static real box_dh(t_psrec *psr)
+{
+ return psr->boxspacing;
+}
+
+#define IS_ONCE (i == nmat-1)
+static real box_dh_top(gmx_bool bOnce, t_psrec *psr)
+{
+ real dh;
+
+ if (psr->bTitle || (psr->bTitleOnce && bOnce) )
+ {
+ dh = 2*psr->titfontsize;
+ }
+ else
+ {
+ dh = 0;
+ }
+
+ return dh;
+}
+
+static gmx_bool box_do_all_x_maj_ticks(t_psrec *psr)
+{
+ return (psr->boxspacing > (1.5*psr->X.majorticklen));
+}
+
+static gmx_bool box_do_all_x_min_ticks(t_psrec *psr)
+{
+ return (psr->boxspacing > (1.5*psr->X.minorticklen));
+}
+
+static void draw_boxes(t_psdata ps, real x0, real y0, real w,
+ int nmat, t_matrix mat[], t_psrec *psr)
+{
+ char buf[128];
+ char *mylab;
+ real xxx;
+ char **xtick, **ytick;
+ real xx, yy, dy, xx00, yy00, offset_x, offset_y;
+ int i, j, x, y, ntx, nty, strlength;
+
+ /* Only necessary when there will be no y-labels */
+ strlength = 0;
+
+ /* Draw the box */
+ ps_rgb(ps, BLACK);
+ ps_linewidth(ps, psr->boxlinewidth);
+ yy00 = y0;
+ for (i = 0; (i < nmat); i++)
+ {
+ dy = box_height(&(mat[i]), psr);
+ ps_box(ps, x0-1, yy00-1, x0+w+1, yy00+dy+1);
+ yy00 += dy+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+ }
+
+ /* Draw the ticks on the axes */
+ ps_linewidth(ps, psr->ticklinewidth);
+ xx00 = x0-1;
+ yy00 = y0-1;
+ for (i = 0; (i < nmat); i++)
+ {
+ if (mat[i].flags & MAT_SPATIAL_X)
+ {
+ ntx = mat[i].nx + 1;
+ offset_x = 0.1;
+ }
+ else
+ {
+ ntx = mat[i].nx;
+ offset_x = 0.6;
+ }
+ if (mat[i].flags & MAT_SPATIAL_Y)
+ {
+ nty = mat[i].ny + 1;
+ offset_y = 0.1;
+ }
+ else
+ {
+ nty = mat[i].ny;
+ offset_y = 0.6;
+ }
+ snew(xtick, ntx);
+ for (j = 0; (j < ntx); j++)
+ {
+ sprintf(buf, "%g", mat[i].axis_x[j]);
+ xtick[j] = strdup(buf);
+ }
+ ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
+ for (x = 0; (x < ntx); x++)
+ {
+ xx = xx00 + (x + offset_x)*psr->xboxsize;
+ if ( ( bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) ||
+ (psr->X.first && (x == 0))) &&
+ ( (i == 0) || box_do_all_x_maj_ticks(psr) ) )
+ {
+ /* Longer tick marks */
+ ps_line (ps, xx, yy00, xx, yy00-psr->X.majorticklen);
+ /* Plot label on lowest graph only */
+ if (i == 0)
+ {
+ ps_ctext(ps, xx,
+ yy00-DDD-psr->X.majorticklen-psr->X.tickfontsize*0.8,
+ xtick[x], eXCenter);
+ }
+ }
+ else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.minor) &&
+ ( (i == 0) || box_do_all_x_min_ticks(psr) ) )
+ {
+ /* Shorter tick marks */
+ ps_line(ps, xx, yy00, xx, yy00-psr->X.minorticklen);
+ }
+ else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) )
+ {
+ /* Even shorter marks, only each X.major */
+ ps_line(ps, xx, yy00, xx, yy00-(psr->boxspacing/2));
+ }
+ }
+ ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
+ snew(ytick, nty);
+ for (j = 0; (j < nty); j++)
+ {
+ sprintf(buf, "%g", mat[i].axis_y[j]);
+ ytick[j] = strdup(buf);
+ }
+
+ for (y = 0; (y < nty); y++)
+ {
+ yy = yy00 + (y + offset_y)*psr->yboxsize;
+ if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.major) ||
+ (psr->Y.first && (y == 0)))
+ {
+ /* Major ticks */
+ strlength = max(strlength, (int)strlen(ytick[y]));
+ ps_line (ps, xx00, yy, xx00-psr->Y.majorticklen, yy);
+ ps_ctext(ps, xx00-psr->Y.majorticklen-DDD,
+ yy-psr->Y.tickfontsize/3.0, ytick[y], eXRight);
+ }
+ else if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.minor) )
+ {
+ /* Minor ticks */
+ ps_line(ps, xx00, yy, xx00-psr->Y.minorticklen, yy);
+ }
+ }
+ sfree(xtick);
+ sfree(ytick);
+
+ /* Label on Y-axis */
+ if (!psr->bYonce || i == nmat/2)
+ {
+ if (strlen(psr->Y.label) > 0)
+ {
+ mylab = psr->Y.label;
+ }
+ else
+ {
+ mylab = mat[i].label_y;
+ }
+ if (strlen(mylab) > 0)
+ {
+ ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
+ ps_flip(ps, TRUE);
+ xxx = x0-psr->X.majorticklen-psr->X.tickfontsize*strlength-DDD;
+ ps_ctext(ps, yy00+box_height(&mat[i], psr)/2.0, 612.5-xxx,
+ mylab, eXCenter);
+ ps_flip(ps, FALSE);
+ }
+ }
+
+ yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+ }
+ /* Label on X-axis */
+ if (strlen(psr->X.label) > 0)
+ {
+ mylab = psr->X.label;
+ }
+ else
+ {
+ mylab = mat[0].label_x;
+ }
+ if (strlen(mylab) > 0)
+ {
+ ps_strfont(ps, psr->X.font, psr->X.fontsize);
+ ps_ctext(ps, x0+w/2, y0-DDD-psr->X.majorticklen-psr->X.tickfontsize*FUDGE-
+ psr->X.fontsize, mylab, eXCenter);
+ }
+}
+
+static void draw_zerolines(t_psdata out, real x0, real y0, real w,
+ int nmat, t_matrix mat[], t_psrec *psr)
+{
+ real xx, yy, dy, xx00, yy00;
+ int i, x, y;
+
+ xx00 = x0-1.5;
+ yy00 = y0-1.5;
+ ps_linewidth(out, psr->zerolinewidth);
+ for (i = 0; (i < nmat); i++)
+ {
+ dy = box_height(&(mat[i]), psr);
+ /* mat[i].axis_x and _y were already set by draw_boxes */
+ if (psr->X.lineatzero)
+ {
+ ps_rgb(out, linecolors[psr->X.lineatzero]);
+ for (x = 0; (x < mat[i].nx); x++)
+ {
+ xx = xx00+(x+0.7)*psr->xboxsize;
+ /* draw lines whenever tick label almost zero (e.g. next trajectory) */
+ if (x != 0 && x < mat[i].nx-1 &&
- abs(mat[i].axis_y[y]) <
- 0.1*abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
++ fabs(mat[i].axis_x[x]) <
++ 0.1*fabs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
+ {
+ ps_line (out, xx, yy00, xx, yy00+dy+2);
+ }
+ }
+ }
+ if (psr->Y.lineatzero)
+ {
+ ps_rgb(out, linecolors[psr->Y.lineatzero]);
+ for (y = 0; (y < mat[i].ny); y++)
+ {
+ yy = yy00+(y+0.7)*psr->yboxsize;
+ /* draw lines whenever tick label almost zero (e.g. next trajectory) */
+ if (y != 0 && y < mat[i].ny-1 &&
- if (abs(mats[i].axis_x[x+1]) < 1e-5)
++ fabs(mat[i].axis_y[y]) <
++ 0.1*fabs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
+ {
+ ps_line (out, xx00, yy, xx00+w+2, yy);
+ }
+ }
+ }
+ yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+ }
+}
+
+static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
+ int elegend, gmx_bool bFrame,
+ real *w, real *h, real *dw, real *dh)
+{
+ int i, maxytick;
+ real ww, hh, dww, dhh;
+
+ hh = dww = dhh = 0;
+ maxytick = 0;
+
+ ww = 0;
+ for (i = 0; (i < nmat); i++)
+ {
+ ww = max(ww, mat[i].nx*psr->xboxsize);
+ hh += box_height(&(mat[i]), psr);
+ maxytick = max(maxytick, mat[i].nx);
+ }
+ if (bFrame)
+ {
+ if (mat[0].label_y[0])
+ {
+ dww += 2.0*(psr->Y.fontsize+DDD);
+ }
+ if (psr->Y.major > 0)
+ {
+ dww += psr->Y.majorticklen + DDD +
+ psr->Y.tickfontsize*(log(maxytick)/log(10.0));
+ }
+ else if (psr->Y.minor > 0)
+ {
+ dww += psr->Y.minorticklen;
+ }
+
+ if (mat[0].label_x[0])
+ {
+ dhh += psr->X.fontsize+2*DDD;
+ }
+ if ( /* fool emacs auto-indent */
+ (elegend == elBoth && (mat[0].legend[0] || mat2[0].legend[0])) ||
+ (elegend == elFirst && mat[0].legend[0]) ||
+ (elegend == elSecond && mat2[0].legend[0]) )
+ {
+ dhh += 2*(psr->legfontsize*FUDGE+2*DDD);
+ }
+ else
+ {
+ dhh += psr->legfontsize*FUDGE+2*DDD;
+ }
+ if (psr->X.major > 0)
+ {
+ dhh += psr->X.tickfontsize*FUDGE+2*DDD+psr->X.majorticklen;
+ }
+ else if (psr->X.minor > 0)
+ {
+ dhh += psr->X.minorticklen;
+ }
+
+ hh += (nmat-1)*box_dh(psr);
+ hh += box_dh_top(TRUE, psr);
+ if (nmat > 1)
+ {
+ hh += (nmat-1)*box_dh_top(FALSE, psr);
+ }
+ }
+ *w = ww;
+ *h = hh;
+ *dw = dww;
+ *dh = dhh;
+}
+
+int add_maps(t_mapping **newmap,
+ int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
+{
+ static char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
+ int nsymbols;
+ int nmap, j, k;
+ t_mapping *map;
+
+ nsymbols = strlen(mapper);
+ nmap = nmap1+nmap2;
+ if (nmap > nsymbols*nsymbols)
+ {
+ gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
+ }
+ printf("Combining colormaps of %d and %d elements into one of %d elements\n",
+ nmap1, nmap2, nmap);
+ snew(map, nmap);
+ for (j = 0; j < nmap1; j++)
+ {
+ map[j].code.c1 = mapper[j % nsymbols];
+ if (nmap > nsymbols)
+ {
+ map[j].code.c2 = mapper[j/nsymbols];
+ }
+ map[j].rgb.r = map1[j].rgb.r;
+ map[j].rgb.g = map1[j].rgb.g;
+ map[j].rgb.b = map1[j].rgb.b;
+ map[j].desc = map1[j].desc;
+ }
+ for (j = 0; j < nmap2; j++)
+ {
+ k = j+nmap1;
+ map[k].code.c1 = mapper[k % nsymbols];
+ if (nmap > nsymbols)
+ {
+ map[k].code.c2 = mapper[k/nsymbols];
+ }
+ map[k].rgb.r = map2[j].rgb.r;
+ map[k].rgb.g = map2[j].rgb.g;
+ map[k].rgb.b = map2[j].rgb.b;
+ map[k].desc = map2[j].desc;
+ }
+
+ *newmap = map;
+ return nmap;
+}
+
+void xpm_mat(const char *outf, int nmat, t_matrix *mat, t_matrix *mat2,
+ gmx_bool bDiag, gmx_bool bFirstDiag)
+{
+ FILE *out;
+ char buf[100];
+ int i, j, k, x, y, col;
+ int nmap;
+ t_mapping *map = NULL;
+
+ out = gmx_ffopen(outf, "w");
+
+ for (i = 0; i < nmat; i++)
+ {
+ if (!mat2 || !diff_maps(mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map))
+ {
+ write_xpm_m(out, mat[0]);
+ }
+ else
+ {
+ nmap = add_maps(&map, mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map);
+ for (x = 0; (x < mat[i].nx); x++)
+ {
+ for (y = 0; (y < mat[i].nx); y++)
+ {
+ if ((x < y) || ((x == y) && bFirstDiag)) /* upper left -> map1 */
+ {
+ col = mat[i].matrix[x][y];
+ }
+ else /* lower right -> map2 */
+ {
+ col = mat[i].nmap+mat[i].matrix[x][y];
+ }
+ if ((bDiag) || (x != y))
+ {
+ mat[i].matrix[x][y] = col;
+ }
+ else
+ {
+ mat[i].matrix[x][y] = 0;
+ }
+ }
+ }
+ sfree(mat[i].map);
+ mat[i].nmap = nmap;
+ mat[i].map = map;
+ if (mat2 && (strcmp(mat[i].title, mat2[i].title) != 0))
+ {
+ sprintf(mat[i].title+strlen(mat[i].title), " / %s", mat2[i].title);
+ }
+ if (mat2 && (strcmp(mat[i].legend, mat2[i].legend) != 0))
+ {
+ sprintf(mat[i].legend+strlen(mat[i].legend), " / %s", mat2[i].legend);
+ }
+ write_xpm_m(out, mat[i]);
+ }
+ }
+ gmx_ffclose(out);
+}
+
+static void tick_spacing(int n, real axis[], real offset, char axisnm,
+ real *major, real *minor)
+{
+ real space;
+ gmx_bool bTryAgain, bFive;
+ int i, j, t, f = 0, ten;
+#define NFACT 4
+ real major_fact[NFACT] = {5, 4, 2, 1};
+ real minor_fact[NFACT] = {5, 4, 4, 5};
+
+ /* start with interval between 10 matrix points: */
+ space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
+ /* get power of 10 */
+ ten = (int)ceil(log(space)/log(10))-1;
+ bTryAgain = TRUE;
+ for (t = ten+2; t > ten-3 && bTryAgain; t--)
+ {
+ for (f = 0; f < NFACT && bTryAgain; f++)
+ {
+ space = pow(10, t) * major_fact[f];
+ /* count how many ticks we would get: */
+ i = 0;
+ for (j = 0; j < n; j++)
+ {
+ if (bRmod(axis[j], offset, space) )
+ {
+ i++;
+ }
+ }
+ /* do we have a reasonable number of ticks ? */
+ bTryAgain = (i > min(10, n-1)) || (i < 5);
+ }
+ }
+ if (bTryAgain)
+ {
+ space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
+ fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n",
+ axisnm, space);
+ }
+ *major = space;
+ *minor = space / minor_fact[f-1];
+ fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n",
+ axisnm, *major, *minor);
+}
+
+void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
+ gmx_bool bFrame, gmx_bool bDiag, gmx_bool bFirstDiag,
+ gmx_bool bTitle, gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
+ real size, real boxx, real boxy, const char *m2p, const char *m2pout,
+ int mapoffset)
+{
+ const char *libm2p;
+ char buf[256], *legend;
+ t_psdata out;
+ t_psrec psrec, *psr;
+ int W, H;
+ int i, j, x, y, col, leg = 0;
+ real x0, y0, xx;
+ real w, h, dw, dh;
+ int nmap1 = 0, nmap2 = 0, leg_nmap;
+ t_mapping *map1 = NULL, *map2 = NULL, *leg_map;
+ gmx_bool bMap1, bNextMap1, bDiscrete;
+
+ /* memory leak: */
+ libm2p = m2p ? gmxlibfn(m2p) : m2p;
+ get_params(libm2p, m2pout, &psrec);
+
+ psr = &psrec;
+
+ if (psr->X.major <= 0)
+ {
+ tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
+ mat[0].axis_x, psr->X.offset, 'X',
+ &(psr->X.major), &(psr->X.minor) );
+ }
+ if (psr->X.minor <= 0)
+ {
+ psr->X.minor = psr->X.major / 2;
+ }
+ if (psr->Y.major <= 0)
+ {
+ tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
+ mat[0].axis_y, psr->Y.offset, 'Y',
+ &(psr->Y.major), &(psr->Y.minor) );
+ }
+ if (psr->Y.minor <= 0)
+ {
+ psr->Y.minor = psr->Y.major / 2;
+ }
+
+ if (boxx > 0)
+ {
+ psr->xboxsize = boxx;
+ psr->yboxsize = boxx;
+ }
+ if (boxy > 0)
+ {
+ psr->yboxsize = boxy;
+ }
+
+ if (psr->xboxsize == 0)
+ {
+ psr->xboxsize = size/mat[0].nx;
+ printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
+ }
+ if (psr->yboxsize == 0)
+ {
+ psr->yboxsize = size/mat[0].nx;
+ printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
+ }
+
+ nmap1 = 0;
+ for (i = 0; (i < nmat); i++)
+ {
+ if (mat[i].nmap > nmap1)
+ {
+ nmap1 = mat[i].nmap;
+ map1 = mat[i].map;
+ leg = i+1;
+ }
+ }
+ if (leg != 1)
+ {
+ printf("Selected legend of matrix # %d for display\n", leg);
+ }
+ if (mat2)
+ {
+ nmap2 = 0;
+ for (i = 0; (i < nmat); i++)
+ {
+ if (mat2[i].nmap > nmap2)
+ {
+ nmap2 = mat2[i].nmap;
+ map2 = mat2[i].map;
+ leg = i+1;
+ }
+ }
+ if (leg != 1)
+ {
+ printf("Selected legend of matrix # %d for second display\n", leg);
+ }
+ }
+ if ( (mat[0].legend[0] == 0) && psr->legend)
+ {
+ strcpy(mat[0].legend, psr->leglabel);
+ }
+
+ bTitle = bTitle && mat[nmat-1].title[0];
+ bTitleOnce = bTitleOnce && mat[nmat-1].title[0];
+ psr->bTitle = bTitle;
+ psr->bTitleOnce = bTitleOnce;
+ psr->bYonce = bYonce;
+
+ /* Set up size of box for nice colors */
+ box_dim(nmat, mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
+
+ /* Set up bounding box */
+ W = w+dw;
+ H = h+dh;
+
+ /* Start box at */
+ x0 = dw;
+ y0 = dh;
+ x = W+psr->xoffs;
+ y = H+psr->yoffs;
+ if (bFrame)
+ {
+ x += 5*DDD;
+ y += 4*DDD;
+ }
+ out = ps_open(outf, 0, 0, x, y);
+ ps_linewidth(out, psr->linewidth);
+ ps_init_rgb_box(out, psr->xboxsize, psr->yboxsize);
+ ps_init_rgb_nbox(out, psr->xboxsize, psr->yboxsize);
+ ps_translate(out, psr->xoffs, psr->yoffs);
+
+ if (bFrame)
+ {
+ ps_comment(out, "Here starts the BOX drawing");
+ draw_boxes(out, x0, y0, w, nmat, mat, psr);
+ }
+
+ for (i = 0; (i < nmat); i++)
+ {
+ if (bTitle || (bTitleOnce && i == nmat-1) )
+ {
+ /* Print title, if any */
+ ps_rgb(out, BLACK);
+ ps_strfont(out, psr->titfont, psr->titfontsize);
+ if (!mat2 || (strcmp(mat[i].title, mat2[i].title) == 0))
+ {
+ strcpy(buf, mat[i].title);
+ }
+ else
+ {
+ sprintf(buf, "%s / %s", mat[i].title, mat2[i].title);
+ }
+ ps_ctext(out, x0+w/2, y0+box_height(&(mat[i]), psr)+psr->titfontsize,
+ buf, eXCenter);
+ }
+ sprintf(buf, "Here starts the filling of box #%d", i);
+ ps_comment(out, buf);
+ for (x = 0; (x < mat[i].nx); x++)
+ {
+ int nexty;
+ int nextcol;
+
+ xx = x0+x*psr->xboxsize;
+ ps_moveto(out, xx, y0);
+ y = 0;
+ bMap1 = (!mat2 || (x < y || (x == y && bFirstDiag)));
+ if ((bDiag) || (x != y))
+ {
+ col = mat[i].matrix[x][y];
+ }
+ else
+ {
+ col = -1;
+ }
+ for (nexty = 1; (nexty <= mat[i].ny); nexty++)
+ {
+ bNextMap1 = (!mat2 || (x < nexty || (x == nexty && bFirstDiag)));
+ /* TRUE: upper left -> map1 */
+ /* FALSE: lower right -> map2 */
+ if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
+ {
+ nextcol = -1;
+ }
+ else
+ {
+ nextcol = mat[i].matrix[x][nexty];
+ }
+ if ( (nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1) )
+ {
+ if (col >= 0)
+ {
+ if (bMap1)
+ {
+ ps_rgb_nbox(out, &(mat[i].map[col].rgb), nexty-y);
+ }
+ else
+ {
+ ps_rgb_nbox(out, &(mat2[i].map[col].rgb), nexty-y);
+ }
+ }
+ else
+ {
+ ps_moverel(out, 0, psr->yboxsize);
+ }
+ y = nexty;
+ bMap1 = bNextMap1;
+ col = nextcol;
+ }
+ }
+ }
+ y0 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
+ }
+
+ if (psr->X.lineatzero || psr->Y.lineatzero)
+ {
+ /* reset y0 for first box */
+ y0 = dh;
+ ps_comment(out, "Here starts the zero lines drawing");
+ draw_zerolines(out, x0, y0, w, nmat, mat, psr);
+ }
+
+ if (elegend != elNone)
+ {
+ ps_comment(out, "Now it's legend time!");
+ ps_linewidth(out, psr->linewidth);
+ if (mat2 == NULL || elegend != elSecond)
+ {
+ bDiscrete = mat[0].bDiscrete;
+ legend = mat[0].legend;
+ leg_nmap = nmap1;
+ leg_map = map1;
+ }
+ else
+ {
+ bDiscrete = mat2[0].bDiscrete;
+ legend = mat2[0].legend;
+ leg_nmap = nmap2;
+ leg_map = map2;
+ }
+ if (bDiscrete)
+ {
+ leg_discrete(out, psr->legfontsize, DDD, legend,
+ psr->legfontsize, psr->legfont, leg_nmap, leg_map);
+ }
+ else
+ {
+ if (elegend != elBoth)
+ {
+ leg_continuous(out, x0+w/2, w/2, DDD, legend,
+ psr->legfontsize, psr->legfont, leg_nmap, leg_map,
+ mapoffset);
+ }
+ else
+ {
+ leg_bicontinuous(out, x0+w/2, w, DDD, mat[0].legend, mat2[0].legend,
+ psr->legfontsize, psr->legfont, nmap1, map1, nmap2, map2);
+ }
+ }
+ ps_comment(out, "Were there, dude");
+ }
+
+ ps_close(out);
+}
+
+void make_axis_labels(int nmat, t_matrix *mat)
+{
+ int i, j;
+
+ for (i = 0; (i < nmat); i++)
+ {
+ /* Make labels for x axis */
+ if (mat[i].axis_x == NULL)
+ {
+ snew(mat[i].axis_x, mat[i].nx);
+ for (j = 0; (j < mat[i].nx); j++)
+ {
+ mat[i].axis_x[j] = j;
+ }
+ }
+ /* Make labels for y axis */
+ if (mat[i].axis_y == NULL)
+ {
+ snew(mat[i].axis_y, mat[i].ny);
+ for (j = 0; (j < mat[i].ny); j++)
+ {
+ mat[i].axis_y[j] = j;
+ }
+ }
+ }
+}
+
+void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
+{
+ int i, x, y, xs, ys;
+
+ for (i = 0; i < nmat; i++)
+ {
+ fprintf(stderr, "converting %dx%d matrix to %dx%d\n",
+ mat[i].nx, mat[i].ny,
+ (mat[i].nx+skip-1)/skip, (mat[i].ny+skip-1)/skip);
+ /* walk through matrix */
+ xs = 0;
+ for (x = 0; (x < mat[i].nx); x++)
+ {
+ if (x % skip == 0)
+ {
+ mat[i].axis_x[xs] = mat[i].axis_x[x];
+ if (mat2)
+ {
+ mat2[i].axis_x[xs] = mat2[i].axis_x[x];
+ }
+ ys = 0;
+ for (y = 0; (y < mat[i].ny); y++)
+ {
+ if (x == 0)
+ {
+ mat[i].axis_y[ys] = mat[i].axis_y[y];
+ if (mat2)
+ {
+ mat2[i].axis_y[ys] = mat2[i].axis_y[y];
+ }
+ }
+ if (y % skip == 0)
+ {
+ mat[i].matrix[xs][ys] = mat[i].matrix[x][y];
+ if (mat2)
+ {
+ mat2[i].matrix[xs][ys] = mat2[i].matrix[x][y];
+ }
+ ys++;
+ }
+ }
+ xs++;
+ }
+ }
+ /* adjust parameters */
+ mat[i].nx = (mat[i].nx+skip-1)/skip;
+ mat[i].ny = (mat[i].ny+skip-1)/skip;
+ if (mat2)
+ {
+ mat2[i].nx = (mat2[i].nx+skip-1)/skip;
+ mat2[i].ny = (mat2[i].ny+skip-1)/skip;
+ }
+ }
+}
+
+void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
+{
+ int i, x, y, m;
+ t_matrix *mats;
+
+ for (i = 0; i < nmat; i++)
+ {
+ for (m = 0; m < (mat2 ? 2 : 1); m++)
+ {
+ if (m == 0)
+ {
+ mats = mat;
+ }
+ else
+ {
+ mats = mat2;
+ }
+ for (x = 0; x < mats[i].nx-1; x++)
+ {
- if (abs(mats[i].axis_y[y+1]) < 1e-5)
++ if (fabs(mats[i].axis_x[x+1]) < 1e-5)
+ {
+ for (y = 0; y < mats[i].ny; y++)
+ {
+ mats[i].matrix[x][y] = 0;
+ }
+ }
+ }
+ for (y = 0; y < mats[i].ny-1; y++)
+ {
++ if (fabs(mats[i].axis_y[y+1]) < 1e-5)
+ {
+ for (x = 0; x < mats[i].nx; x++)
+ {
+ mats[i].matrix[x][y] = 0;
+ }
+ }
+ }
+ }
+ }
+}
+
+void write_combined_matrix(int ecombine, const char *fn,
+ int nmat, t_matrix *mat1, t_matrix *mat2,
+ real *cmin, real *cmax)
+{
+ int i, j, k, nlevels;
+ t_mapping *map = NULL;
+ FILE *out;
+ real **rmat1, **rmat2;
+ real rhi, rlo;
+
+ out = gmx_ffopen(fn, "w");
+ for (k = 0; k < nmat; k++)
+ {
+ if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
+ {
+ gmx_fatal(FARGS, "Size of frame %d in 1st (%dx%d) and 2nd matrix (%dx%d) do"
+ " not match.\n", k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
+ }
+ printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
+ rmat1 = matrix2real(&mat1[k], NULL);
+ rmat2 = matrix2real(&mat2[k], NULL);
+ if (NULL == rmat1 || NULL == rmat2)
+ {
+ gmx_fatal(FARGS, "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
+ "g_rms and g_mdmat provide such data, but not do_dssp.\n",
+ (NULL == rmat1 && NULL == rmat2) ? "both" : "one of the" );
+ }
+ rlo = 1e38;
+ rhi = -1e38;
+ for (j = 0; j < mat1[k].ny; j++)
+ {
+ for (i = 0; i < mat1[k].nx; i++)
+ {
+ switch (ecombine)
+ {
+ case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
+ case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
+ case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
+ case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
+ default:
+ gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
+ }
+ rlo = min(rlo, rmat1[i][j]);
+ rhi = max(rhi, rmat1[i][j]);
+ }
+ }
+ if (cmin)
+ {
+ rlo = *cmin;
+ }
+ if (cmax)
+ {
+ rhi = *cmax;
+ }
+ nlevels = max(mat1[k].nmap, mat2[k].nmap);
+ if (rhi == rlo)
+ {
+ fprintf(stderr,
+ "combination results in uniform matrix (%g), no output\n", rhi);
+ }
+ /*
+ else if (rlo>=0 || rhi<=0)
+ write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
+ mat1[k].label_x, mat1[k].label_y,
+ mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
+ rmat1, rlo, rhi, rhi<=0?red:white, rhi<=0?white:blue,
+ &nlevels);
+ else
+ write_xpm3(out, mat2[k].flags, mat1[k].title, mat1[k].legend,
+ mat1[k].label_x, mat1[k].label_y,
+ mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
+ rmat1, rlo, 0, rhi, red, white, blue, &nlevels);
+ */
+ else
+ {
+ write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
+ mat1[k].label_x, mat1[k].label_y,
+ mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
+ rmat1, rlo, rhi, white, black, &nlevels);
+ }
+ }
+ gmx_ffclose(out);
+}
+
+void do_mat(int nmat, t_matrix *mat, t_matrix *mat2,
+ gmx_bool bFrame, gmx_bool bZeroLine, gmx_bool bDiag, gmx_bool bFirstDiag, gmx_bool bTitle,
+ gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
+ real size, real boxx, real boxy,
+ const char *epsfile, const char *xpmfile, const char *m2p,
+ const char *m2pout, int skip, int mapoffset)
+{
+ int i, j, k;
+
+ if (mat2)
+ {
+ for (k = 0; (k < nmat); k++)
+ {
+ if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
+ {
+ gmx_fatal(FARGS, "WAKE UP!! Size of frame %d in 2nd matrix file (%dx%d) does not match size of 1st matrix (%dx%d) or the other way around.\n",
+ k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
+ }
+ for (j = 0; (j < mat[k].ny); j++)
+ {
+ for (i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
+ {
+ mat[k].matrix[i][j] = mat2[k].matrix[i][j];
+ }
+ }
+ }
+ }
+ for (i = 0; (i < nmat); i++)
+ {
+ fprintf(stderr, "Matrix %d is %d x %d\n", i, mat[i].nx, mat[i].ny);
+ }
+
+ make_axis_labels(nmat, mat);
+
+ if (skip > 1)
+ {
+ prune_mat(nmat, mat, mat2, skip);
+ }
+
+ if (bZeroLine)
+ {
+ zero_lines(nmat, mat, mat);
+ }
+
+ if (epsfile != NULL)
+ {
+ ps_mat(epsfile, nmat, mat, mat2, bFrame, bDiag, bFirstDiag,
+ bTitle, bTitleOnce, bYonce, elegend,
+ size, boxx, boxy, m2p, m2pout, mapoffset);
+ }
+ if (xpmfile != NULL)
+ {
+ xpm_mat(xpmfile, nmat, mat, mat2, bDiag, bFirstDiag);
+ }
+}
+
+void gradient_map(rvec grad, int nmap, t_mapping map[])
+{
+ int i;
+ real c;
+
+ for (i = 0; i < nmap; i++)
+ {
+ c = i/(nmap-1.0);
+ map[i].rgb.r = 1-c*(1-grad[XX]);
+ map[i].rgb.g = 1-c*(1-grad[YY]);
+ map[i].rgb.b = 1-c*(1-grad[ZZ]);
+ }
+}
+
+void gradient_mat(rvec grad, int nmat, t_matrix mat[])
+{
+ int m;
+
+ for (m = 0; m < nmat; m++)
+ {
+ gradient_map(grad, mat[m].nmap, mat[m].map);
+ }
+}
+
+void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
+{
+ int i;
+ real c, r, g, b;
+
+ for (i = 0; i < nmap; i++)
+ {
+ c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
+ if (c > 1)
+ {
+ c = 1;
+ }
+ if (bBlue)
+ {
+ c = 1 - c;
+ }
+ if (c <= 0.25) /* 0-0.25 */
+ {
+ r = 0;
+ g = pow(4*c, 0.666);
+ b = 1;
+ }
+ else if (c <= 0.5) /* 0.25-0.5 */
+ {
+ r = 0;
+ g = 1;
+ b = pow(2-4*c, 0.666);
+ }
+ else if (c <= 0.75) /* 0.5-0.75 */
+ {
+ r = pow(4*c-2, 0.666);
+ g = 1;
+ b = 0;
+ }
+ else /* 0.75-1 */
+ {
+ r = 1;
+ g = pow(4-4*c, 0.666);
+ b = 0;
+ }
+ map[i].rgb.r = r;
+ map[i].rgb.g = g;
+ map[i].rgb.b = b;
+ }
+}
+
+void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
+{
+ int m;
+
+ for (m = 0; m < nmat; m++)
+ {
+ rainbow_map(bBlue, mat[m].nmap, mat[m].map);
+ }
+}
+
+int gmx_xpm2ps(int argc, char *argv[])
+{
+ const char *desc[] = {
+ "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
+ "Labels and axis can be displayed, when they are supplied",
+ "in the correct matrix format.",
+ "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
+ "[gmx-mdmat].[PAR]",
+ "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
+ "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
+ "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
+ "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
+ "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
+ "sets yfont, ytickfont and xtickfont.[PAR]",
+ "When no [TT].m2p[tt] file is supplied, many settings are taken from",
+ "command line options. The most important option is [TT]-size[tt],",
+ "which sets the size of the whole matrix in postscript units.",
+ "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
+ "options (and the corresponding parameters in the [TT].m2p[tt] file),",
+ "which set the size of a single matrix element.[PAR]",
+ "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
+ "files will be read simultaneously and the upper left half of the",
+ "first one ([TT]-f[tt]) is plotted together with the lower right",
+ "half of the second one ([TT]-f2[tt]). The diagonal will contain",
+ "values from the matrix file selected with [TT]-diag[tt].",
+ "Plotting of the diagonal values can be suppressed altogether by",
+ "setting [TT]-diag[tt] to [TT]none[tt].",
+ "In this case, a new color map will be generated with",
+ "a red gradient for negative numbers and a blue for positive.",
+ "If the color coding and legend labels of both matrices are identical,",
+ "only one legend will be displayed, else two separate legends are",
+ "displayed.",
+ "With [TT]-combine[tt], an alternative operation can be selected",
+ "to combine the matrices. The output range is automatically set",
+ "to the actual range of the combined matrix. This can be overridden",
+ "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
+ "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
+ "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
+ "the [IT]y[it]-axis).[PAR]",
+ "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
+ "into attractive color pictures.[PAR]",
+ "Merged or rainbowed matrices can be written to an XPixelMap file with",
+ "the [TT]-xpm[tt] option."
+ };
+
+ output_env_t oenv;
+ const char *fn, *epsfile = NULL, *xpmfile = NULL;
+ int i, nmat, nmat2, etitle, elegend, ediag, erainbow, ecombine;
+ t_matrix *mat = NULL, *mat2 = NULL;
+ gmx_bool bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
+ static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE, bAdd = FALSE;
+ static real size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
+ static rvec grad = {0, 0, 0};
+ enum {
+ etSel, etTop, etOnce, etYlabel, etNone, etNR
+ };
+ const char *title[] = { NULL, "top", "once", "ylabel", "none", NULL };
+ /* MUST correspond to enum elXxx as defined at top of file */
+ const char *legend[] = { NULL, "both", "first", "second", "none", NULL };
+ enum {
+ edSel, edFirst, edSecond, edNone, edNR
+ };
+ const char *diag[] = { NULL, "first", "second", "none", NULL };
+ enum {
+ erSel, erNo, erBlue, erRed, erNR
+ };
+ const char *rainbow[] = { NULL, "no", "blue", "red", NULL };
+ /* MUST correspond to enum ecXxx as defined at top of file */
+ const char *combine[] = {
+ NULL, "halves", "add", "sub", "mult", "div", NULL
+ };
+ static int skip = 1, mapoffset = 0;
+ t_pargs pa[] = {
+ { "-frame", FALSE, etBOOL, {&bFrame},
+ "Display frame, ticks, labels, title and legend" },
+ { "-title", FALSE, etENUM, {title}, "Show title at" },
+ { "-yonce", FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
+ { "-legend", FALSE, etENUM, {legend}, "Show legend" },
+ { "-diag", FALSE, etENUM, {diag}, "Diagonal" },
+ { "-size", FALSE, etREAL, {&size},
+ "Horizontal size of the matrix in ps units" },
+ { "-bx", FALSE, etREAL, {&boxx},
+ "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
+ { "-by", FALSE, etREAL, {&boxy}, "Element y-size" },
+ { "-rainbow", FALSE, etENUM, {rainbow},
+ "Rainbow colors, convert white to" },
+ { "-gradient", FALSE, etRVEC, {grad},
+ "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
+ { "-skip", FALSE, etINT, {&skip},
+ "only write out every nr-th row and column" },
+ { "-zeroline", FALSE, etBOOL, {&bZeroLine},
+ "insert line in [TT].xpm[tt] matrix where axis label is zero"},
+ { "-legoffset", FALSE, etINT, {&mapoffset},
+ "Skip first N colors from [TT].xpm[tt] file for the legend" },
+ { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
+ { "-cmin", FALSE, etREAL, {&cmin}, "Minimum for combination output" },
+ { "-cmax", FALSE, etREAL, {&cmax}, "Maximum for combination output" }
+ };
+#define NPA asize(pa)
+ t_filenm fnm[] = {
+ { efXPM, "-f", NULL, ffREAD },
+ { efXPM, "-f2", "root2", ffOPTRD },
+ { efM2P, "-di", NULL, ffLIBOPTRD },
+ { efM2P, "-do", "out", ffOPTWR },
+ { efEPS, "-o", NULL, ffOPTWR },
+ { efXPM, "-xpm", NULL, ffOPTWR }
+ };
+#define NFILE asize(fnm)
+
+ if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
+ NFILE, fnm, NPA, pa,
+ asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+
+ etitle = nenum(title);
+ elegend = nenum(legend);
+ ediag = nenum(diag);
+ erainbow = nenum(rainbow);
+ ecombine = nenum(combine);
+ bGrad = opt2parg_bSet("-gradient", NPA, pa);
+ for (i = 0; i < DIM; i++)
+ {
+ if (grad[i] < 0 || grad[i] > 1)
+ {
+ gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
+ }
+ }
+ if (!bFrame)
+ {
+ etitle = etNone;
+ elegend = elNone;
+ }
+
+ epsfile = ftp2fn_null(efEPS, NFILE, fnm);
+ xpmfile = opt2fn_null("-xpm", NFILE, fnm);
+ if (epsfile == NULL && xpmfile == NULL)
+ {
+ if (ecombine != ecHalves)
+ {
+ xpmfile = opt2fn("-xpm", NFILE, fnm);
+ }
+ else
+ {
+ epsfile = ftp2fn(efEPS, NFILE, fnm);
+ }
+ }
+ if (ecombine != ecHalves && epsfile)
+ {
+ fprintf(stderr,
+ "WARNING: can only write result of arithmetic combination "
+ "of two matrices to .xpm file\n"
+ " file %s will not be written\n", epsfile);
+ epsfile = NULL;
+ }
+
+ bDiag = ediag != edNone;
+ bFirstDiag = ediag != edSecond;
+
+ fn = opt2fn("-f", NFILE, fnm);
+ nmat = read_xpm_matrix(fn, &mat);
+ fprintf(stderr, "There %s %d matri%s in %s\n", (nmat > 1) ? "are" : "is", nmat, (nmat > 1) ? "ces" : "x", fn);
+ fn = opt2fn_null("-f2", NFILE, fnm);
+ if (fn)
+ {
+ nmat2 = read_xpm_matrix(fn, &mat2);
+ fprintf(stderr, "There %s %d matri%s in %s\n", (nmat2 > 1) ? "are" : "is", nmat2, (nmat2 > 1) ? "ces" : "x", fn);
+ if (nmat != nmat2)
+ {
+ fprintf(stderr, "Different number of matrices, using the smallest number.\n");
+ nmat = nmat2 = min(nmat, nmat2);
+ }
+ }
+ else
+ {
+ if (ecombine != ecHalves)
+ {
+ fprintf(stderr,
+ "WARNING: arithmetic matrix combination selected (-combine), "
+ "but no second matrix (-f2) supplied\n"
+ " no matrix combination will be performed\n");
+ }
+ ecombine = 0;
+ nmat2 = 0;
+ }
+ bTitle = etitle == etTop;
+ bTitleOnce = etitle == etOnce;
+ if (etitle == etYlabel)
+ {
+ for (i = 0; (i < nmat); i++)
+ {
+ strcpy(mat[i].label_y, mat[i].title);
+ if (mat2)
+ {
+ strcpy(mat2[i].label_y, mat2[i].title);
+ }
+ }
+ }
+ if (bGrad)
+ {
+ gradient_mat(grad, nmat, mat);
+ if (mat2)
+ {
+ gradient_mat(grad, nmat2, mat2);
+ }
+ }
+ if (erainbow != erNo)
+ {
+ rainbow_mat(erainbow == erBlue, nmat, mat);
+ if (mat2)
+ {
+ rainbow_mat(erainbow == erBlue, nmat2, mat2);
+ }
+ }
+
+ if ((mat2 == NULL) && (elegend != elNone))
+ {
+ elegend = elFirst;
+ }
+
+ if (ecombine && ecombine != ecHalves)
+ {
+ write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
+ opt2parg_bSet("-cmin", NPA, pa) ? &cmin : NULL,
+ opt2parg_bSet("-cmax", NPA, pa) ? &cmax : NULL);
+ }
+ else
+ {
+ do_mat(nmat, mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
+ bTitle, bTitleOnce, bYonce,
+ elegend, size, boxx, boxy, epsfile, xpmfile,
+ opt2fn_null("-di", NFILE, fnm), opt2fn_null("-do", NFILE, fnm), skip,
+ mapoffset);
+ }
+
+ view_all(oenv, NFILE, fnm);
+
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff