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