From: Roland Schulz Date: Fri, 9 May 2014 19:08:58 +0000 (-0400) Subject: Merge release-4-6 into release-5-0 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=73eb6fd38236ec350495096244e82ca44f873561;p=alexxy%2Fgromacs.git Merge release-4-6 into release-5-0 Conflicts: CMakeLists.txt (applied to cmake/gmxManageSharedLibraries.cmake) Replaced gmx_llabs by llabs as mentioned in TODO. Change-Id: I37e8cf5d6725cc9836e6dc29ca44e7ef915417c4 --- 73eb6fd38236ec350495096244e82ca44f873561 diff --cc cmake/gmxManageSharedLibraries.cmake index cc73de31b4,0000000000..fced521e95 mode 100644,000000..100644 --- a/cmake/gmxManageSharedLibraries.cmake +++ b/cmake/gmxManageSharedLibraries.cmake @@@ -1,140 -1,0 +1,140 @@@ +# +# 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") AND NOT GMX_BUILD_MDRUN_ONLY) ++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() diff --cc src/gromacs/fileio/libxdrf.c index 219e64ead9,0000000000..14c7501900 mode 100644,000000..100644 --- a/src/gromacs/fileio/libxdrf.c +++ b/src/gromacs/fileio/libxdrf.c @@@ -1,1669 -1,0 +1,1667 @@@ +/* + * 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 +#endif + +#include +#include +#include +#include +#include + +#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; +} + - +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 (fr != frame && 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 ((((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 ((((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 (abs(low - high) <= header_size) ++ 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; +} diff --cc src/gromacs/gmxana/geminate.c index 79dc949fd3,0000000000..a9ee8113de mode 100644,000000..100644 --- a/src/gromacs/gmxana/geminate.c +++ b/src/gromacs/gmxana/geminate.c @@@ -1,767 -1,0 +1,767 @@@ +/* + * 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 +#endif + +#include +#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)=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)|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. */ - if (abs(ct[0]-1.0) > 1e-6) ++ if (fabs(ct[0]-1.0) > 1e-6) + { + ct[0] = 1.0; - fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(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; + } + } + } +} diff --cc src/gromacs/gmxana/gmx_anaeig.c index f59b1130d8,0000000000..053e9c585f mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_anaeig.c +++ b/src/gromacs/gmxana/gmx_anaeig.c @@@ -1,1487 -1,0 +1,1487 @@@ +/* + * 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 +#endif +#include +#include +#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(x[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 (bSplit && i > 0 && 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 && abs(inprod[noutvec][i]) < 1e-5) ++ 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; +} diff --cc src/gromacs/gmxana/gmx_current.c index c706b40be5,0000000000..ec1a988f5d mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_current.c +++ b/src/gromacs/gmxana/gmx_current.c @@@ -1,1008 -1,0 +1,1008 @@@ +/* + * 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 +#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 (abs(qall) > 0.01) ++ 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_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; +} diff --cc src/gromacs/gmxana/gmx_mindist.c index 2ec58643ef,0000000000..311e3d6873 mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_mindist.c +++ b/src/gromacs/gmxana/gmx_mindist.c @@@ -1,798 -1,0 +1,798 @@@ +/* + * 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 +#endif + +#include +#include + +#include "sysstuff.h" +#include +#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 && abs(t/output_env_get_time_factor(oenv)) < 1e-5) ++ 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; +} diff --cc src/gromacs/gmxana/gmx_rms.c index 13aeb6c807,0000000000..3962646be6 mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_rms.c +++ b/src/gromacs/gmxana/gmx_rms.c @@@ -1,1223 -1,0 +1,1223 @@@ +/* + * 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 +#endif + +#include "gromacs/utility/smalloc.h" +#include +#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 && - abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 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 && abs(time[i]) < 1e-5) ++ 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; +} diff --cc src/gromacs/gmxana/gmx_xpm2ps.c index 064e667c3c,0000000000..d8c5ebe56d mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_xpm2ps.c +++ b/src/gromacs/gmxana/gmx_xpm2ps.c @@@ -1,1619 -1,0 +1,1619 @@@ +/* + * 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 +#endif + +#include +#include +#include +#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_x[x]) < - 0.1*abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) ) ++ 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 && - abs(mat[i].axis_y[y]) < - 0.1*abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) ) ++ 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_x[x+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 (abs(mats[i].axis_y[y+1]) < 1e-5) ++ 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; +}