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