src/gromacs/gmxana/eigio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "eigio.h"
  40
  41 #include "gromacs/fileio/tpxio.h"
  42 #include "gromacs/fileio/trrio.h"
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/utility/futil.h"
  45 #include "gromacs/utility/smalloc.h"
  46
  47 void read_eigenvectors(const char *file, int *natoms, gmx_bool *bFit,
  48                        rvec **xref, gmx_bool *bDMR,
  49                        rvec **xav, gmx_bool *bDMA,
  50                        int *nvec, int **eignr,
  51                        rvec ***eigvec, real **eigval)
  52 {
  53     gmx_trr_header_t   head;
  54     int                i, snew_size;
  55     struct t_fileio   *status;
  56     rvec              *x;
  57     matrix             box;
  58     gmx_bool           bOK;
  59
  60     *bDMR = FALSE;
  61
  62     /* read (reference (t=-1) and) average (t=0) structure */
  63     status = gmx_trr_open(file, "r");
  64     gmx_trr_read_frame_header(status, &head, &bOK);
  65     *natoms = head.natoms;
  66     snew(*xav, *natoms);
  67     gmx_trr_read_frame_data(status, &head, box, *xav, nullptr, nullptr);
  68
  69     if ((head.t >= -1.1) && (head.t <= -0.9))
  70     {
  71         snew(*xref, *natoms);
  72         for (i = 0; i < *natoms; i++)
  73         {
  74             copy_rvec((*xav)[i], (*xref)[i]);
  75         }
  76         *bDMR = (head.lambda > 0.5);
  77         *bFit = (head.lambda > -0.5);
  78         if (*bFit)
  79         {
  80             fprintf(stderr, "Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ", *natoms, file);
  81         }
  82         else
  83         {
  84             fprintf(stderr, "Eigenvectors in %s were determined without fitting\n", file);
  85             sfree(*xref);
  86             *xref = nullptr;
  87         }
  88         gmx_trr_read_frame_header(status, &head, &bOK);
  89         gmx_trr_read_frame_data(status, &head, box, *xav, nullptr, nullptr);
  90     }
  91     else
  92     {
  93         *bFit = TRUE;
  94         *xref = nullptr;
  95     }
  96     *bDMA = (head.lambda > 0.5);
  97     if ((head.t <= -0.01) || (head.t >= 0.01))
  98     {
  99         fprintf(stderr, "WARNING: %s does not start with t=0, which should be the "
 100                 "average structure. This might not be a eigenvector file. "
 101                 "Some things might go wrong.\n",
 102                 file);
 103     }
 104     else
 105     {
 106         fprintf(stderr,
 107                 "Read %smass weighted average/minimum structure with %d atoms from %s\n",
 108                 *bDMA ? "" : "non ", *natoms, file);
 109     }
 110
 111     snew(x, *natoms);
 112     snew_size = 10;
 113     snew(*eignr, snew_size);
 114     snew(*eigval, snew_size);
 115     snew(*eigvec, snew_size);
 116
 117     *nvec = 0;
 118     while (gmx_trr_read_frame_header(status, &head, &bOK))
 119     {
 120         gmx_trr_read_frame_data(status, &head, box, x, nullptr, nullptr);
 121         if (*nvec >= snew_size)
 122         {
 123             snew_size += 10;
 124             srenew(*eignr, snew_size);
 125             srenew(*eigval, snew_size);
 126             srenew(*eigvec, snew_size);
 127         }
 128         i                = head.step;
 129         (*eigval)[*nvec] = head.t;
 130         (*eignr)[*nvec]  = i-1;
 131         snew((*eigvec)[*nvec], *natoms);
 132         for (i = 0; i < *natoms; i++)
 133         {
 134             copy_rvec(x[i], (*eigvec)[*nvec][i]);
 135         }
 136         (*nvec)++;
 137     }
 138     sfree(x);
 139     gmx_trr_close(status);
 140     fprintf(stderr, "Read %d eigenvectors (for %d atoms)\n\n", *nvec, *natoms);
 141 }
 142
 143
 144 void write_eigenvectors(const char *trrname, int natoms, const real mat[],
 145                         gmx_bool bReverse, int begin, int end,
 146                         int WriteXref, const rvec *xref, gmx_bool bDMR,
 147                         const rvec xav[], gmx_bool bDMA, const real eigval[])
 148 {
 149     struct t_fileio *trrout;
 150     int              ndim, i, j, d, vec;
 151     matrix           zerobox;
 152     rvec            *x;
 153
 154     ndim = natoms*DIM;
 155     clear_mat(zerobox);
 156     snew(x, natoms);
 157
 158     fprintf (stderr,
 159              "\nWriting %saverage structure & eigenvectors %d--%d to %s\n",
 160              (WriteXref == eWXR_YES) ? "reference, " : "",
 161              begin, end, trrname);
 162
 163     trrout = gmx_trr_open(trrname, "w");
 164     if (WriteXref == eWXR_YES)
 165     {
 166         /* misuse lambda: 0/1 mass weighted fit no/yes */
 167         gmx_trr_write_frame(trrout, -1, -1, bDMR ? 1.0 : 0.0, zerobox, natoms, xref, nullptr, nullptr);
 168     }
 169     else if (WriteXref == eWXR_NOFIT)
 170     {
 171         /* misuse lambda: -1 no fit */
 172         gmx_trr_write_frame(trrout, -1, -1, -1.0, zerobox, natoms, x, nullptr, nullptr);
 173     }
 174
 175     /* misuse lambda: 0/1 mass weighted analysis no/yes */
 176     gmx_trr_write_frame(trrout, 0, 0, bDMA ? 1.0 : 0.0, zerobox, natoms, xav, nullptr, nullptr);
 177
 178     for (i = 0; i <= (end-begin); i++)
 179     {
 180
 181         if (!bReverse)
 182         {
 183             vec = i;
 184         }
 185         else
 186         {
 187             vec = ndim-i-1;
 188         }
 189
 190         for (j = 0; j < natoms; j++)
 191         {
 192             for (d = 0; d < DIM; d++)
 193             {
 194                 x[j][d] = mat[vec*ndim+DIM*j+d];
 195             }
 196         }
 197
 198         /* Store the eigenvalue in the time field */
 199         gmx_trr_write_frame(trrout, begin+i, eigval[vec], 0, zerobox, natoms, x, nullptr, nullptr);
 200     }
 201     gmx_trr_close(trrout);
 202
 203     sfree(x);
 204 }