src/tools/eigio.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "smalloc.h"
  40 #include "vec.h"
  41 #include "eigio.h"
  42 #include "trnio.h"
  43 #include "tpxio.h"
  44 #include "futil.h"
  45
  46 void read_eigenvectors(const char *file,int *natoms,bool *bFit,
  47                        rvec **xref,bool *bDMR,
  48                        rvec **xav,bool *bDMA,
  49                        int *nvec, int **eignr,
  50                        rvec ***eigvec,real **eigval)
  51 {
  52   t_trnheader head;
  53   int    status,i,snew_size;
  54   rvec   *x;
  55   matrix box;
  56   bool   bOK;
  57
  58   *bDMR=FALSE;
  59
  60   /* read (reference (t=-1) and) average (t=0) structure */
  61   status=open_trn(file,"r");
  62   fread_trnheader(status,&head,&bOK);
  63   *natoms=head.natoms;
  64   snew(*xav,*natoms);
  65   fread_htrn(status,&head,box,*xav,NULL,NULL);
  66
  67   if ((head.t>=-1.1) && (head.t<=-0.9))
  68   {
  69       snew(*xref,*natoms);
  70       for(i=0; i<*natoms; i++)
  71           copy_rvec((*xav)[i],(*xref)[i]);
  72       *bDMR = (head.lambda > 0.5);
  73       *bFit = (head.lambda > -0.5);
  74       if (*bFit)
  75       {
  76           fprintf(stderr,"Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ",*natoms,file);
  77       }
  78       else
  79       {
  80           fprintf(stderr,"Eigenvectors in %s were determined without fitting\n",file);
  81           sfree(*xref);
  82           *xref=NULL;
  83       }
  84       fread_trnheader(status,&head,&bOK);
  85       fread_htrn(status,&head,box,*xav,NULL,NULL);
  86   }
  87   else
  88   {
  89       *bFit=TRUE;
  90       *xref=NULL;
  91   }
  92   *bDMA = (head.lambda > 0.5);
  93   if ((head.t<=-0.01) || (head.t>=0.01))
  94   {
  95       fprintf(stderr,"WARNING: %s does not start with t=0, which should be the "
  96               "average structure. This might not be a eigenvector file. "
  97               "Some things might go wrong.\n",
  98               file);
  99   }
 100   else
 101   {
 102       fprintf(stderr,
 103               "Read %smass weighted average/minimum structure with %d atoms from %s\n",
 104               *bDMA ? "" : "non ",*natoms,file);
 105   }
 106
 107   snew(x,*natoms);
 108   snew_size=10;
 109   snew(*eignr,snew_size);
 110   snew(*eigval,snew_size);
 111   snew(*eigvec,snew_size);
 112
 113   *nvec=0;
 114   while (fread_trnheader(status,&head,&bOK))
 115   {
 116       fread_htrn(status,&head,box,x,NULL,NULL);
 117       if (*nvec >= snew_size)
 118       {
 119           snew_size+=10;
 120           srenew(*eignr,snew_size);
 121           srenew(*eigval,snew_size);
 122           srenew(*eigvec,snew_size);
 123       }
 124       i=head.step;
 125       (*eigval)[*nvec]=head.t;
 126       (*eignr)[*nvec]=i-1;
 127       snew((*eigvec)[*nvec],*natoms);
 128       for(i=0; i<*natoms; i++)
 129       {
 130           copy_rvec(x[i],(*eigvec)[*nvec][i]);
 131       }
 132       (*nvec)++;
 133   }
 134   sfree(x);
 135   fprintf(stderr,"Read %d eigenvectors (for %d atoms)\n\n",*nvec,*natoms);
 136 }
 137
 138
 139 void write_eigenvectors(const char *trnname,int natoms,real mat[],
 140                         bool bReverse,int begin,int end,
 141                         int WriteXref,rvec *xref,bool bDMR,
 142                         rvec xav[], bool bDMA,real eigval[])
 143 {
 144     int    trnout;
 145     int    ndim,i,j,d,vec;
 146     matrix zerobox;
 147     rvec   *x;
 148
 149     ndim = natoms*DIM;
 150     clear_mat(zerobox);
 151     snew(x,natoms);
 152
 153     fprintf (stderr,
 154              "\nWriting %saverage structure & eigenvectors %d--%d to %s\n",
 155              (WriteXref==eWXR_YES) ? "reference, " : "",
 156              begin,end,trnname);
 157
 158     trnout = open_tpx(trnname,"w");
 159     if (WriteXref == eWXR_YES)
 160     {
 161         /* misuse lambda: 0/1 mass weighted fit no/yes */
 162         fwrite_trn(trnout,-1,-1,bDMR ? 1.0 : 0.0,zerobox,natoms,xref,NULL,NULL);
 163     }
 164     else if (WriteXref == eWXR_NOFIT)
 165     {
 166         /* misuse lambda: -1 no fit */
 167         fwrite_trn(trnout,-1,-1,-1.0,zerobox,natoms,x,NULL,NULL);
 168     }
 169
 170     /* misuse lambda: 0/1 mass weighted analysis no/yes */
 171     fwrite_trn(trnout,0,0,bDMA ? 1.0 : 0.0,zerobox,natoms,xav,NULL,NULL);
 172
 173     for(i=0; i<=(end-begin); i++)
 174     {
 175
 176         if (!bReverse)
 177             vec = i;
 178         else
 179             vec = ndim-i-1;
 180
 181         for (j=0; j<natoms; j++)
 182             for(d=0; d<DIM; d++)
 183                 x[j][d]=mat[vec*ndim+DIM*j+d];
 184
 185         /* Store the eigenvalue in the time field */
 186         fwrite_trn(trnout,begin+i,eigval[vec],0,zerobox,natoms,x,NULL,NULL);
 187     }
 188     close_trn(trnout);
 189
 190     sfree(x);
 191 }
 192
 193