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