src/tools/gmx_covar.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 <math.h>
  40 #include <string.h>
  41 #include <time.h>
  42
  43 #ifdef HAVE_SYS_TIME_H
  44 #include <sys/time.h>
  45 #endif
  46
  47 #include "statutil.h"
  48 #include "sysstuff.h"
  49 #include "typedefs.h"
  50 #include "smalloc.h"
  51 #include "macros.h"
  52 #include "vec.h"
  53 #include "pbc.h"
  54 #include "copyrite.h"
  55 #include "futil.h"
  56 #include "statutil.h"
  57 #include "index.h"
  58 #include "confio.h"
  59 #include "trnio.h"
  60 #include "mshift.h"
  61 #include "xvgr.h"
  62 #include "do_fit.h"
  63 #include "rmpbc.h"
  64 #include "txtdump.h"
  65 #include "matio.h"
  66 #include "eigio.h"
  67 #include "physics.h"
  68 #include "gmx_ana.h"
  69 #include "string2.h"
  70
  71 #include "gromacs/linearalgebra/eigensolver.h"
  72
  73 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
  74 char *
  75 gmx_ctime_r(const time_t *clock,char *buf, int n);
  76
  77
  78 int gmx_covar(int argc,char *argv[])
  79 {
  80   const char *desc[] = {
  81     "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
  82     "covariance matrix.",
  83     "All structures are fitted to the structure in the structure file.",
  84     "When this is not a run input file periodicity will not be taken into",
  85     "account. When the fit and analysis groups are identical and the analysis",
  86     "is non mass-weighted, the fit will also be non mass-weighted.",
  87     "[PAR]",
  88     "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
  89     "When the same atoms are used for the fit and the covariance analysis,",
  90     "the reference structure for the fit is written first with t=-1.",
  91     "The average (or reference when [TT]-ref[tt] is used) structure is",
  92     "written with t=0, the eigenvectors",
  93     "are written as frames with the eigenvector number as timestamp.",
  94     "[PAR]",
  95     "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
  96     "[PAR]",
  97     "Option [TT]-ascii[tt] writes the whole covariance matrix to",
  98     "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
  99     "[PAR]",
 100     "Option [TT]-xpm[tt] writes the whole covariance matrix to an [TT].xpm[tt] file.",
 101     "[PAR]",
 102     "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [TT].xpm[tt] file,",
 103     "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
 104     "written.",
 105     "[PAR]",
 106     "Note that the diagonalization of a matrix requires memory and time",
 107     "that will increase at least as fast as than the square of the number",
 108     "of atoms involved. It is easy to run out of memory, in which",
 109     "case this tool will probably exit with a 'Segmentation fault'. You",
 110     "should consider carefully whether a reduced set of atoms will meet",
 111     "your needs for lower costs."
 112   };
 113   static gmx_bool bFit=TRUE,bRef=FALSE,bM=FALSE,bPBC=TRUE;
 114   static int  end=-1;
 115   t_pargs pa[] = {
 116     { "-fit",  FALSE, etBOOL, {&bFit},
 117       "Fit to a reference structure"},
 118     { "-ref",  FALSE, etBOOL, {&bRef},
 119       "Use the deviation from the conformation in the structure file instead of from the average" },
 120     { "-mwa",  FALSE, etBOOL, {&bM},
 121       "Mass-weighted covariance analysis"},
 122     { "-last",  FALSE, etINT, {&end},
 123       "Last eigenvector to write away (-1 is till the last)" },
 124     { "-pbc",  FALSE,  etBOOL, {&bPBC},
 125       "Apply corrections for periodic boundary conditions" }
 126   };
 127   FILE       *out;
 128   t_trxstatus *status;
 129   t_trxstatus *trjout;
 130   t_topology top;
 131   int        ePBC;
 132   t_atoms    *atoms;
 133   rvec       *x,*xread,*xref,*xav,*xproj;
 134   matrix     box,zerobox;
 135   real       *sqrtm,*mat,*eigenvalues,sum,trace,inv_nframes;
 136   real       t,tstart,tend,**mat2;
 137   real       xj,*w_rls=NULL;
 138   real       min,max,*axis;
 139   int        ntopatoms,step;
 140   int        natoms,nat,count,nframes0,nframes,nlevels;
 141   gmx_large_int_t ndim,i,j,k,l;
 142   int        WriteXref;
 143   const char *fitfile,*trxfile,*ndxfile;
 144   const char *eigvalfile,*eigvecfile,*averfile,*logfile;
 145   const char *asciifile,*xpmfile,*xpmafile;
 146   char       str[STRLEN],*fitname,*ananame,*pcwd;
 147   int        d,dj,nfit;
 148   atom_id    *index,*ifit;
 149   gmx_bool       bDiffMass1,bDiffMass2;
 150   time_t     now;
 151   char       timebuf[STRLEN];
 152   t_rgb      rlo,rmi,rhi;
 153   real       *eigenvectors;
 154   output_env_t oenv;
 155   gmx_rmpbc_t  gpbc=NULL;
 156
 157   t_filenm fnm[] = {
 158     { efTRX, "-f",  NULL, ffREAD },
 159     { efTPS, NULL,  NULL, ffREAD },
 160     { efNDX, NULL,  NULL, ffOPTRD },
 161     { efXVG, NULL,  "eigenval", ffWRITE },
 162     { efTRN, "-v",  "eigenvec", ffWRITE },
 163     { efSTO, "-av", "average.pdb", ffWRITE },
 164     { efLOG, NULL,  "covar", ffWRITE },
 165     { efDAT, "-ascii","covar", ffOPTWR },
 166     { efXPM, "-xpm","covar", ffOPTWR },
 167     { efXPM, "-xpma","covara", ffOPTWR }
 168   };
 169 #define NFILE asize(fnm)
 170
 171   CopyRight(stderr,argv[0]);
 172   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 173                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 174
 175   clear_mat(zerobox);
 176
 177   fitfile    = ftp2fn(efTPS,NFILE,fnm);
 178   trxfile    = ftp2fn(efTRX,NFILE,fnm);
 179   ndxfile    = ftp2fn_null(efNDX,NFILE,fnm);
 180   eigvalfile = ftp2fn(efXVG,NFILE,fnm);
 181   eigvecfile = ftp2fn(efTRN,NFILE,fnm);
 182   averfile   = ftp2fn(efSTO,NFILE,fnm);
 183   logfile    = ftp2fn(efLOG,NFILE,fnm);
 184   asciifile  = opt2fn_null("-ascii",NFILE,fnm);
 185   xpmfile    = opt2fn_null("-xpm",NFILE,fnm);
 186   xpmafile   = opt2fn_null("-xpma",NFILE,fnm);
 187
 188   read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE);
 189   atoms=&top.atoms;
 190
 191   if (bFit) {
 192     printf("\nChoose a group for the least squares fit\n");
 193     get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname);
 194     if (nfit < 3)
 195       gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
 196   } else
 197     nfit=0;
 198   printf("\nChoose a group for the covariance analysis\n");
 199   get_index(atoms,ndxfile,1,&natoms,&index,&ananame);
 200
 201   bDiffMass1=FALSE;
 202   if (bFit) {
 203     snew(w_rls,atoms->nr);
 204     for(i=0; (i<nfit); i++) {
 205       w_rls[ifit[i]]=atoms->atom[ifit[i]].m;
 206       if (i)
 207         bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]]!=w_rls[ifit[i-1]]);
 208     }
 209   }
 210   bDiffMass2=FALSE;
 211   snew(sqrtm,natoms);
 212   for(i=0; (i<natoms); i++)
 213     if (bM) {
 214       sqrtm[i]=sqrt(atoms->atom[index[i]].m);
 215       if (i)
 216         bDiffMass2 = bDiffMass2 || (sqrtm[i]!=sqrtm[i-1]);
 217     }
 218     else
 219       sqrtm[i]=1.0;
 220
 221   if (bFit && bDiffMass1 && !bDiffMass2) {
 222     bDiffMass1 = natoms != nfit;
 223     i=0;
 224     for (i=0; (i<natoms) && !bDiffMass1; i++)
 225       bDiffMass1 = index[i] != ifit[i];
 226     if (!bDiffMass1) {
 227       fprintf(stderr,"\n"
 228               "Note: the fit and analysis group are identical,\n"
 229               "      while the fit is mass weighted and the analysis is not.\n"
 230               "      Making the fit non mass weighted.\n\n");
 231       for(i=0; (i<nfit); i++)
 232         w_rls[ifit[i]]=1.0;
 233     }
 234   }
 235
 236   /* Prepare reference frame */
 237   if (bPBC) {
 238     gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box);
 239     gmx_rmpbc(gpbc,atoms->nr,box,xref);
 240   }
 241   if (bFit)
 242     reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
 243
 244   snew(x,natoms);
 245   snew(xav,natoms);
 246   ndim=natoms*DIM;
 247   if (sqrt(GMX_LARGE_INT_MAX)<ndim) {
 248     gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
 249   }
 250   snew(mat,ndim*ndim);
 251
 252   fprintf(stderr,"Calculating the average structure ...\n");
 253   nframes0 = 0;
 254   nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
 255   if (nat != atoms->nr)
 256     fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);
 257   do {
 258     nframes0++;
 259     /* calculate x: a fitted struture of the selected atoms */
 260     if (bPBC)
 261       gmx_rmpbc(gpbc,nat,box,xread);
 262     if (bFit) {
 263       reset_x(nfit,ifit,nat,NULL,xread,w_rls);
 264       do_fit(nat,w_rls,xref,xread);
 265     }
 266     for (i=0; i<natoms; i++)
 267       rvec_inc(xav[i],xread[index[i]]);
 268   } while (read_next_x(oenv,status,&t,nat,xread,box));
 269   close_trj(status);
 270
 271   inv_nframes = 1.0/nframes0;
 272   for(i=0; i<natoms; i++)
 273     for(d=0; d<DIM; d++) {
 274       xav[i][d] *= inv_nframes;
 275       xread[index[i]][d] = xav[i][d];
 276     }
 277   write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
 278                          atoms,xread,NULL,epbcNONE,zerobox,natoms,index);
 279   sfree(xread);
 280
 281   fprintf(stderr,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim,(int)ndim);
 282   nframes=0;
 283   nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
 284   tstart = t;
 285   do {
 286     nframes++;
 287     tend = t;
 288     /* calculate x: a (fitted) structure of the selected atoms */
 289     if (bPBC)
 290       gmx_rmpbc(gpbc,nat,box,xread);
 291     if (bFit) {
 292       reset_x(nfit,ifit,nat,NULL,xread,w_rls);
 293       do_fit(nat,w_rls,xref,xread);
 294     }
 295     if (bRef)
 296       for (i=0; i<natoms; i++)
 297         rvec_sub(xread[index[i]],xref[index[i]],x[i]);
 298     else
 299       for (i=0; i<natoms; i++)
 300         rvec_sub(xread[index[i]],xav[i],x[i]);
 301
 302     for (j=0; j<natoms; j++) {
 303       for (dj=0; dj<DIM; dj++) {
 304         k=ndim*(DIM*j+dj);
 305         xj=x[j][dj];
 306         for (i=j; i<natoms; i++) {
 307           l=k+DIM*i;
 308           for(d=0; d<DIM; d++)
 309             mat[l+d] += x[i][d]*xj;
 310         }
 311       }
 312     }
 313   } while (read_next_x(oenv,status,&t,nat,xread,box) &&
 314            (bRef || nframes < nframes0));
 315   close_trj(status);
 316   gmx_rmpbc_done(gpbc);
 317
 318   fprintf(stderr,"Read %d frames\n",nframes);
 319
 320   if (bRef) {
 321     /* copy the reference structure to the ouput array x */
 322     snew(xproj,natoms);
 323     for (i=0; i<natoms; i++)
 324       copy_rvec(xref[index[i]],xproj[i]);
 325   } else {
 326     xproj = xav;
 327   }
 328
 329   /* correct the covariance matrix for the mass */
 330   inv_nframes = 1.0/nframes;
 331   for (j=0; j<natoms; j++)
 332     for (dj=0; dj<DIM; dj++)
 333       for (i=j; i<natoms; i++) {
 334         k = ndim*(DIM*j+dj)+DIM*i;
 335         for (d=0; d<DIM; d++)
 336           mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
 337       }
 338
 339   /* symmetrize the matrix */
 340   for (j=0; j<ndim; j++)
 341     for (i=j; i<ndim; i++)
 342       mat[ndim*i+j]=mat[ndim*j+i];
 343
 344   trace=0;
 345   for(i=0; i<ndim; i++)
 346     trace+=mat[i*ndim+i];
 347   fprintf(stderr,"\nTrace of the covariance matrix: %g (%snm^2)\n",
 348           trace,bM ? "u " : "");
 349
 350   if (asciifile) {
 351     out = ffopen(asciifile,"w");
 352     for (j=0; j<ndim; j++) {
 353       for (i=0; i<ndim; i+=3)
 354         fprintf(out,"%g %g %g\n",
 355                 mat[ndim*j+i],mat[ndim*j+i+1],mat[ndim*j+i+2]);
 356     }
 357     ffclose(out);
 358   }
 359
 360   if (xpmfile) {
 361     min = 0;
 362     max = 0;
 363     snew(mat2,ndim);
 364     for (j=0; j<ndim; j++) {
 365       mat2[j] = &(mat[ndim*j]);
 366       for (i=0; i<=j; i++) {
 367         if (mat2[j][i] < min)
 368           min = mat2[j][i];
 369         if (mat2[j][j] > max)
 370           max = mat2[j][i];
 371       }
 372     }
 373     snew(axis,ndim);
 374     for(i=0; i<ndim; i++)
 375       axis[i] = i+1;
 376     rlo.r = 0; rlo.g = 0; rlo.b = 1;
 377     rmi.r = 1; rmi.g = 1; rmi.b = 1;
 378     rhi.r = 1; rhi.g = 0; rhi.b = 0;
 379     out = ffopen(xpmfile,"w");
 380     nlevels = 80;
 381     write_xpm3(out,0,"Covariance",bM ? "u nm^2" : "nm^2",
 382                "dim","dim",ndim,ndim,axis,axis,
 383                mat2,min,0.0,max,rlo,rmi,rhi,&nlevels);
 384     ffclose(out);
 385     sfree(axis);
 386     sfree(mat2);
 387   }
 388
 389   if (xpmafile) {
 390     min = 0;
 391     max = 0;
 392     snew(mat2,ndim/DIM);
 393     for (i=0; i<ndim/DIM; i++)
 394       snew(mat2[i],ndim/DIM);
 395     for (j=0; j<ndim/DIM; j++) {
 396       for (i=0; i<=j; i++) {
 397         mat2[j][i] = 0;
 398         for(d=0; d<DIM; d++)
 399           mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
 400         if (mat2[j][i] < min)
 401           min = mat2[j][i];
 402         if (mat2[j][j] > max)
 403           max = mat2[j][i];
 404         mat2[i][j] = mat2[j][i];
 405       }
 406     }
 407     snew(axis,ndim/DIM);
 408     for(i=0; i<ndim/DIM; i++)
 409       axis[i] = i+1;
 410     rlo.r = 0; rlo.g = 0; rlo.b = 1;
 411     rmi.r = 1; rmi.g = 1; rmi.b = 1;
 412     rhi.r = 1; rhi.g = 0; rhi.b = 0;
 413     out = ffopen(xpmafile,"w");
 414     nlevels = 80;
 415     write_xpm3(out,0,"Covariance",bM ? "u nm^2" : "nm^2",
 416                "atom","atom",ndim/DIM,ndim/DIM,axis,axis,
 417                mat2,min,0.0,max,rlo,rmi,rhi,&nlevels);
 418     ffclose(out);
 419     sfree(axis);
 420     for (i=0; i<ndim/DIM; i++)
 421       sfree(mat2[i]);
 422     sfree(mat2);
 423   }
 424
 425
 426   /* call diagonalization routine */
 427
 428   snew(eigenvalues,ndim);
 429   snew(eigenvectors,ndim*ndim);
 430
 431   memcpy(eigenvectors,mat,ndim*ndim*sizeof(real));
 432   fprintf(stderr,"\nDiagonalizing ...\n");
 433   fflush(stderr);
 434   eigensolver(eigenvectors,ndim,0,ndim,eigenvalues,mat);
 435   sfree(eigenvectors);
 436
 437   /* now write the output */
 438
 439   sum=0;
 440   for(i=0; i<ndim; i++)
 441     sum+=eigenvalues[i];
 442   fprintf(stderr,"\nSum of the eigenvalues: %g (%snm^2)\n",
 443           sum,bM ? "u " : "");
 444   if (fabs(trace-sum)>0.01*trace)
 445     fprintf(stderr,"\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
 446
 447   fprintf(stderr,"\nWriting eigenvalues to %s\n",eigvalfile);
 448
 449   sprintf(str,"(%snm\\S2\\N)",bM ? "u " : "");
 450   out=xvgropen(eigvalfile,
 451                "Eigenvalues of the covariance matrix",
 452                "Eigenvector index",str,oenv);
 453   for (i=0; (i<ndim); i++)
 454     fprintf (out,"%10d %g\n",(int)i+1,eigenvalues[ndim-1-i]);
 455   ffclose(out);
 456
 457   if (end==-1) {
 458     if (nframes-1 < ndim)
 459       end=nframes-1;
 460     else
 461       end=ndim;
 462   }
 463   if (bFit) {
 464     /* misuse lambda: 0/1 mass weighted analysis no/yes */
 465     if (nfit==natoms) {
 466       WriteXref = eWXR_YES;
 467       for(i=0; i<nfit; i++)
 468         copy_rvec(xref[ifit[i]],x[i]);
 469     } else
 470       WriteXref = eWXR_NO;
 471   } else {
 472     /* misuse lambda: -1 for no fit */
 473     WriteXref = eWXR_NOFIT;
 474   }
 475
 476   write_eigenvectors(eigvecfile,natoms,mat,TRUE,1,end,
 477                      WriteXref,x,bDiffMass1,xproj,bM,eigenvalues);
 478
 479   out = ffopen(logfile,"w");
 480
 481   time(&now);
 482   gmx_ctime_r(&now,timebuf,STRLEN);
 483   fprintf(out,"Covariance analysis log, written %s\n",timebuf);
 484
 485   fprintf(out,"Program: %s\n",argv[0]);
 486   gmx_getcwd(str, STRLEN);
 487
 488   fprintf(out,"Working directory: %s\n\n",str);
 489
 490   fprintf(out,"Read %d frames from %s (time %g to %g %s)\n",nframes,trxfile,
 491           output_env_conv_time(oenv,tstart),output_env_conv_time(oenv,tend),output_env_get_time_unit(oenv));
 492   if (bFit)
 493     fprintf(out,"Read reference structure for fit from %s\n",fitfile);
 494   if (ndxfile)
 495     fprintf(out,"Read index groups from %s\n",ndxfile);
 496   fprintf(out,"\n");
 497
 498   fprintf(out,"Analysis group is '%s' (%d atoms)\n",ananame,natoms);
 499   if (bFit)
 500     fprintf(out,"Fit group is '%s' (%d atoms)\n",fitname,nfit);
 501   else
 502     fprintf(out,"No fit was used\n");
 503   fprintf(out,"Analysis is %smass weighted\n", bDiffMass2 ? "":"non-");
 504   if (bFit)
 505     fprintf(out,"Fit is %smass weighted\n", bDiffMass1 ? "":"non-");
 506   fprintf(out,"Diagonalized the %dx%d covariance matrix\n",(int)ndim,(int)ndim);
 507   fprintf(out,"Trace of the covariance matrix before diagonalizing: %g\n",
 508           trace);
 509   fprintf(out,"Trace of the covariance matrix after diagonalizing: %g\n\n",
 510           sum);
 511
 512   fprintf(out,"Wrote %d eigenvalues to %s\n",(int)ndim,eigvalfile);
 513   if (WriteXref == eWXR_YES)
 514     fprintf(out,"Wrote reference structure to %s\n",eigvecfile);
 515   fprintf(out,"Wrote average structure to %s and %s\n",averfile,eigvecfile);
 516   fprintf(out,"Wrote eigenvectors %d to %d to %s\n",1,end,eigvecfile);
 517
 518   ffclose(out);
 519
 520   fprintf(stderr,"Wrote the log to %s\n",logfile);
 521
 522   thanx(stderr);
 523
 524   return 0;
 525 }