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