src/tools/g_nmeig.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_g_nmeig_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "statutil.h"
  34 #include "sysstuff.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "vec.h"
  39 #include "pbc.h"
  40 #include "copyrite.h"
  41 #include "futil.h"
  42 #include "statutil.h"
  43 #include "rdgroup.h"
  44 #include "mshift.h"
  45 #include "xvgr.h"
  46 #include "gstat.h"
  47 #include "txtdump.h"
  48 #include "ql77.h"
  49 #include "eigio.h"
  50
  51 int main(int argc,char *argv[])
  52 {
  53   static char *desc[] = {
  54     "g_nmeig calculates the eigenvectors/values of a (Hessian) matrix,",
  55     "which can be calculated with [TT]nmrun[tt].",
  56     "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
  57     "The structure is written first with t=0. The eigenvectors",
  58     "are written as frames with the eigenvector number as timestamp.",
  59     "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
  60     "An ensemble of structures can be generated from the eigenvectors with",
  61     "[TT]g_nmens[tt]."
  62   };
  63   static bool bM=TRUE;
  64   static int  begin=1,end=100;
  65   t_pargs pa[] = {
  66     { "-m",  FALSE, etBOOL, {&bM},
  67       "Divide elements of Hessian by product of sqrt(mass) of involved atoms prior to diagonalization. This should be used for 'Normal Modes' analyses" },
  68     { "-first", FALSE, etINT, {&begin},
  69       "First eigenvector to write away" },
  70     { "-last",  FALSE, etINT, {&end},
  71       "Last eigenvector to write away" }
  72   };
  73   FILE       *out;
  74   int        status,trjout;
  75   t_topology top;
  76   t_inputrec ir;
  77   rvec       *top_x,*x;
  78   matrix     box;
  79   real       t,*hess,*eigv,rdum,mass_fac;
  80   int        natoms,ndim,count;
  81   char       *grpname,title[256];
  82   int        i,j,k,l,d,gnx;
  83   bool       bSuck;
  84   atom_id    *index;
  85   t_filenm fnm[] = {
  86     { efMTX, "-f", "hessian",  ffREAD  },
  87     { efTPS, NULL, NULL,       ffREAD  },
  88     { efXVG, NULL, "eigenval", ffWRITE },
  89     { efTRN, "-v", "eigenvec", ffWRITE }
  90   };
  91 #define NFILE asize(fnm)
  92
  93   CopyRight(stderr,argv[0]);
  94   parse_common_args(&argc,argv,PCA_SET_NPRI,TRUE,
  95                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
  96
  97   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&top_x,NULL,box,bM);
  98
  99   /*open Hessian matrix and read first 'line' */
 100
 101   natoms=read_first_x(&status,ftp2fn(efMTX,NFILE,fnm),&t,&x,box);
 102   ndim=natoms*DIM;
 103
 104   fprintf(stderr,"Dimensionality of matrix: %d\n",ndim);
 105
 106   /* store two dimensional Hessian matrix in one dimensional array
 107    * to facilitate communication with fortran */
 108
 109   snew(hess,ndim*ndim);
 110
 111   snew(eigv,ndim);
 112
 113   for (i=0; (i<natoms); i++) {
 114     for (j=0; (j<DIM); j++) {
 115       hess[i*DIM+j]=0.0-x[i][j];
 116     }
 117   }
 118
 119   /* read in rest of Hessian */
 120
 121   for (i=1; (i<ndim); i++) {
 122     if (read_next_x(status,&t,natoms,x,box)) {
 123       for (j=0; (j<natoms); j++) {
 124         for (k=0; (k<DIM); k++) {
 125           hess[i*ndim+j*DIM+k]=0.0-x[j][k];
 126         }
 127       }
 128     }
 129     else {
 130       fatal_error(0,"Premature end of file, hessian matrix");
 131     }
 132   }
 133   close_trj(status);
 134
 135
 136   /* check if matrix is approximately symmetric */
 137   for (i=0; (i<ndim); i++) {
 138     for (j=0; (j<i); j++) {
 139       if (hess[i*ndim+j] != 0.0 && hess[j*ndim+i] != 0.0)
 140         rdum=(hess[i*ndim+j]-hess[j*ndim+i])/(hess[i*ndim+j]+hess[j*ndim+i]);
 141       else {
 142         if (fabs(hess[i*ndim+j] - hess[j*ndim+i]) > 1.0e-5)
 143           rdum=1.0;
 144         else
 145           rdum=0.0;
 146       }
 147       if (abs(rdum)>1.0e-5) {
 148         fprintf(stderr,"possible non-symmetric matrix:\n");
 149         fprintf(stderr,"x[%d][%d]=%e,x[%d][%d]=%e\n",i,j,\
 150                 hess[i*ndim+j],j,i,hess[j*ndim+i]);
 151       }
 152       /* make sure that it is symmetric
 153        * this is not very nice, but we did warn, did't we? */
 154       hess[j*ndim+i]=hess[i*ndim+j];
 155     }
 156   }
 157
 158   /* divide elements hess[i][j]
 159      by sqrt(mas[i])*sqrt(mas[j]) when required */
 160
 161   if (bM) {
 162     for (i=0; (i<natoms); i++) {
 163       for (j=0; (j<DIM); j++) {
 164         for (k=0; (k<natoms); k++) {
 165           mass_fac=invsqrt(top.atoms.atom[i].m*top.atoms.atom[k].m);
 166           for (l=0; (l<DIM); l++)
 167             hess[(i*DIM+j)*ndim+k*DIM+l]*=mass_fac;
 168         }
 169       }
 170     }
 171   }
 172
 173   /* call diagonalization routine. Tested only fortran double precision */
 174
 175   fprintf(stderr,"\nDiagonalizing...\n");
 176   fflush(stderr);
 177
 178   ql77 (ndim,hess,eigv);
 179
 180   /* check the output, first 6 eigenvalues should be quite small */
 181
 182   bSuck=FALSE;
 183   for (i=0; (i<6); i++) {
 184     if (abs(eigv[i]) > 1.0e-3)
 185       bSuck=TRUE;
 186   }
 187   if (bSuck) {
 188     fprintf(stderr,"\nOne of the first eigenvalues has a non-zero value.\n");
 189     fprintf(stderr,"This could mean that the reference structure was not\n");
 190     fprintf(stderr,"properly energy minimised\n");
 191   }
 192
 193   /* now write the output */
 194   fprintf (stderr,"Writing eigenvalues...\n");
 195   out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 196                "Eigenvalues","Eigenvector index","Eigenvalue");
 197   if (bM)
 198     fprintf(out,"@ subtitle \"of mass weighted Hessian matrix\"\n");
 199   else
 200     fprintf(out,"@ subtitle \"of Hessian matrix\"\n");
 201   for (i=0; i<ndim; i++)
 202     fprintf (out,"%6d %15g\n",i+1,eigv[i]);
 203   fclose(out);
 204
 205   if (begin<1)
 206     begin=1;
 207   if (end>ndim)
 208     end=ndim;
 209
 210   write_eigenvectors(opt2fn("-v",NFILE,fnm),natoms,hess,FALSE,begin,end,
 211                      eWXR_NO,NULL,FALSE,top_x,bM);
 212
 213   thanx(stdout);
 214
 215   return 0;
 216 }