src/gromacs/gmxana/eigio.c

   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, 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 "config.h"
  40
  41 #include "gromacs/utility/smalloc.h"
  42 #include "gromacs/math/vec.h"
  43 #include "eigio.h"
  44 #include "gromacs/fileio/trnio.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/utility/futil.h"
  47
  48 void read_eigenvectors(const char *file, int *natoms, gmx_bool *bFit,
  49                        rvec **xref, gmx_bool *bDMR,
  50                        rvec **xav, gmx_bool *bDMA,
  51                        int *nvec, int **eignr,
  52                        rvec ***eigvec, real **eigval)
  53 {
  54     t_trnheader head;
  55     int         i, snew_size;
  56     t_fileio   *status;
  57     rvec       *x;
  58     matrix      box;
  59     gmx_bool    bOK;
  60
  61     *bDMR = FALSE;
  62
  63     /* read (reference (t=-1) and) average (t=0) structure */
  64     status = open_trn(file, "r");
  65     fread_trnheader(status, &head, &bOK);
  66     *natoms = head.natoms;
  67     snew(*xav, *natoms);
  68     fread_htrn(status, &head, box, *xav, NULL, NULL);
  69
  70     if ((head.t >= -1.1) && (head.t <= -0.9))
  71     {
  72         snew(*xref, *natoms);
  73         for (i = 0; i < *natoms; i++)
  74         {
  75             copy_rvec((*xav)[i], (*xref)[i]);
  76         }
  77         *bDMR = (head.lambda > 0.5);
  78         *bFit = (head.lambda > -0.5);
  79         if (*bFit)
  80         {
  81             fprintf(stderr, "Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ", *natoms, file);
  82         }
  83         else
  84         {
  85             fprintf(stderr, "Eigenvectors in %s were determined without fitting\n", file);
  86             sfree(*xref);
  87             *xref = NULL;
  88         }
  89         fread_trnheader(status, &head, &bOK);
  90         fread_htrn(status, &head, box, *xav, NULL, NULL);
  91     }
  92     else
  93     {
  94         *bFit = TRUE;
  95         *xref = NULL;
  96     }
  97     *bDMA = (head.lambda > 0.5);
  98     if ((head.t <= -0.01) || (head.t >= 0.01))
  99     {
 100         fprintf(stderr, "WARNING: %s does not start with t=0, which should be the "
 101                 "average structure. This might not be a eigenvector file. "
 102                 "Some things might go wrong.\n",
 103                 file);
 104     }
 105     else
 106     {
 107         fprintf(stderr,
 108                 "Read %smass weighted average/minimum structure with %d atoms from %s\n",
 109                 *bDMA ? "" : "non ", *natoms, file);
 110     }
 111
 112     snew(x, *natoms);
 113     snew_size = 10;
 114     snew(*eignr, snew_size);
 115     snew(*eigval, snew_size);
 116     snew(*eigvec, snew_size);
 117
 118     *nvec = 0;
 119     while (fread_trnheader(status, &head, &bOK))
 120     {
 121         fread_htrn(status, &head, box, x, NULL, NULL);
 122         if (*nvec >= snew_size)
 123         {
 124             snew_size += 10;
 125             srenew(*eignr, snew_size);
 126             srenew(*eigval, snew_size);
 127             srenew(*eigvec, snew_size);
 128         }
 129         i                = head.step;
 130         (*eigval)[*nvec] = head.t;
 131         (*eignr)[*nvec]  = i-1;
 132         snew((*eigvec)[*nvec], *natoms);
 133         for (i = 0; i < *natoms; i++)
 134         {
 135             copy_rvec(x[i], (*eigvec)[*nvec][i]);
 136         }
 137         (*nvec)++;
 138     }
 139     sfree(x);
 140     fprintf(stderr, "Read %d eigenvectors (for %d atoms)\n\n", *nvec, *natoms);
 141 }
 142
 143
 144 void write_eigenvectors(const char *trnname, int natoms, real mat[],
 145                         gmx_bool bReverse, int begin, int end,
 146                         int WriteXref, rvec *xref, gmx_bool bDMR,
 147                         rvec xav[], gmx_bool bDMA, real eigval[])
 148 {
 149     t_fileio *trnout;
 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, trnname);
 162
 163     trnout = open_tpx(trnname, "w");
 164     if (WriteXref == eWXR_YES)
 165     {
 166         /* misuse lambda: 0/1 mass weighted fit no/yes */
 167         fwrite_trn(trnout, -1, -1, bDMR ? 1.0 : 0.0, zerobox, natoms, xref, NULL, NULL);
 168     }
 169     else if (WriteXref == eWXR_NOFIT)
 170     {
 171         /* misuse lambda: -1 no fit */
 172         fwrite_trn(trnout, -1, -1, -1.0, zerobox, natoms, x, NULL, NULL);
 173     }
 174
 175     /* misuse lambda: 0/1 mass weighted analysis no/yes */
 176     fwrite_trn(trnout, 0, 0, bDMA ? 1.0 : 0.0, zerobox, natoms, xav, NULL, NULL);
 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         fwrite_trn(trnout, begin+i, eigval[vec], 0, zerobox, natoms, x, NULL, NULL);
 200     }
 201     close_trn(trnout);
 202
 203     sfree(x);
 204 }