src/tools/eigio.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_eigio_c = "$Id$";
  30
  31 #include "smalloc.h"
  32 #include "vec.h"
  33 #include "eigio.h"
  34 #include "trnio.h"
  35 #include "tpxio.h"
  36
  37 void read_eigenvectors(char *file,int *natoms,bool *bFit,
  38                        rvec **xref,bool *bDMR,
  39                        rvec **xav,bool *bDMA,
  40                        int *nvec, int **eignr, rvec ***eigvec)
  41 {
  42   t_trnheader head;
  43   int    status,i,snew_size;
  44   rvec   *x;
  45   matrix box;
  46   bool   bOK;
  47
  48   *bDMR=FALSE;
  49
  50   /* read (reference (t=-1) and) average (t=0) structure */
  51   status=open_trn(file,"r");
  52   fread_trnheader(status,&head,&bOK);
  53   *natoms=head.natoms;
  54   snew(*xav,*natoms);
  55   fread_htrn(status,&head,box,*xav,NULL,NULL);
  56   if ((head.t>=-1.1) && (head.t<=-0.9)) {
  57     snew(*xref,*natoms);
  58     for(i=0; i<*natoms; i++)
  59       copy_rvec((*xav)[i],(*xref)[i]);
  60     *bDMR = (head.lambda > 0.5);
  61     *bFit = (head.lambda > -0.5);
  62     if (*bFit)
  63       fprintf(stderr,"Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ",*natoms,file);
  64     else {
  65        fprintf(stderr,"Eigenvectors in %s were determined without fitting\n",
  66                file);
  67        sfree(*xref);
  68        *xref=NULL;
  69     }
  70     fread_trnheader(status,&head,&bOK);
  71     fread_htrn(status,&head,box,*xav,NULL,NULL);
  72   }
  73   else {
  74     *bFit=TRUE;
  75     *xref=NULL;
  76   }
  77   *bDMA = (head.lambda > 0.5);
  78   if ((head.t<=-0.01) || (head.t>=0.01))
  79     fprintf(stderr,"WARNING: %s does not start with t=0, which should be the "
  80             "average structure. This might not be a eigenvector file. "
  81             "Some things might go wrong.\n",
  82             file);
  83   else
  84     fprintf(stderr,
  85             "Read %smass weighted average/minimum structure with %d atoms from %s\n",
  86             *bDMA ? "" : "non ",*natoms,file);
  87
  88   snew(x,*natoms);
  89   snew_size=0;
  90   *nvec=0;
  91   while (fread_trnheader(status,&head,&bOK)) {
  92     fread_htrn(status,&head,box,x,NULL,NULL);
  93     if (*nvec >= snew_size) {
  94       snew_size+=10;
  95       srenew(*eignr,snew_size);
  96       srenew(*eigvec,snew_size);
  97     }
  98     i=(int)(head.t+0.01);
  99     if ((head.t-i<=-0.01) || (head.t-i>=0.01))
 100       fatal_error(0,"%s contains a frame with non-integer time (%f), this "
 101                   "time should be an eigenvector index. "
 102                   "This might not be a eigenvector file.",file,head.t);
 103     (*eignr)[*nvec]=i-1;
 104     snew((*eigvec)[*nvec],*natoms);
 105     for(i=0; i<*natoms; i++)
 106       copy_rvec(x[i],(*eigvec)[*nvec][i]);
 107     (*nvec)++;
 108   }
 109   sfree(x);
 110   fprintf(stderr,"Read %d eigenvectors (dim=%d)\n\n",*nvec,*natoms*DIM);
 111 }
 112
 113 void write_eigenvectors(char *trnname,int natoms,real mat[],
 114                         bool bReverse,int begin,int end,
 115                         int WriteXref,rvec *xref,bool bDMR,
 116                         rvec xav[],bool bDMA)
 117 {
 118   int    trnout;
 119   int    ndim,i,j,d,vec;
 120   matrix zerobox;
 121   rvec   *x;
 122
 123   ndim = natoms*DIM;
 124   clear_mat(zerobox);
 125   snew(x,natoms);
 126
 127   fprintf (stderr,
 128            "\nWriting %saverage structure\nand eigenvectors 1 to %d to %s\n",
 129            (WriteXref==eWXR_YES) ? "reference and " : "",
 130            end,trnname);
 131
 132   trnout = open_tpx(trnname,"w");
 133   if (WriteXref == eWXR_YES)
 134     /* misuse lambda: 0/1 mass weighted fit no/yes */
 135     fwrite_trn(trnout,-1,-1,bDMR ? 1.0 : 0.0,zerobox,natoms,xref,NULL,NULL);
 136   else if (WriteXref == eWXR_NOFIT)
 137     /* misuse lambda: -1 no fit */
 138     fwrite_trn(trnout,-1,-1,-1.0,zerobox,natoms,x,NULL,NULL);
 139
 140   /* misuse lambda: 0/1 mass weighted analysis no/yes */
 141   fwrite_trn(trnout,0,0,bDMA ? 1.0 : 0.0,zerobox,natoms,xav,NULL,NULL);
 142
 143   for(i=begin; i<=end; i++) {
 144     if (!bReverse)
 145       vec = i-1;
 146     else
 147       vec = ndim-i;
 148     for (j=0; j<natoms; j++)
 149       for(d=0; d<DIM; d++)
 150         x[j][d]=mat[vec*ndim+DIM*j+d];
 151     fwrite_trn(trnout,i,(real)i,0,zerobox,natoms,x,NULL,NULL);
 152   }
 153   close_trn(trnout);
 154
 155   sfree(x);
 156 }