3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
42 #ifdef HAVE_SYS_TIME_H
47 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
72 #include "eigensolver.h"
77 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
79 gmx_ctime_r(const time_t *clock,char *buf, int n);
82 int gmx_covar(int argc,char *argv[])
84 const char *desc[] = {
85 "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
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.",
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.",
99 "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
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, ...",
104 "Option [TT]-xpm[tt] writes the whole covariance matrix to an xpm file.",
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",
110 static bool bFit=TRUE,bRef=FALSE,bM=FALSE,bPBC=TRUE;
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" }
130 rvec *x,*xread,*xref,*xav,*xproj;
132 real *sqrtm,*mat,*eigval,sum,trace,inv_nframes;
133 real t,tstart,tend,**mat2;
137 int natoms,nat,count,nframes0,nframes,nlevels;
138 gmx_large_int_t ndim,i,j,k,l;
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;
145 atom_id *index,*ifit;
146 bool bDiffMass1,bDiffMass2;
148 char timebuf[STRLEN];
152 gmx_rmpbc_t gpbc=NULL;
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 }
166 #define NFILE asize(fnm)
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);
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);
185 read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE);
189 printf("\nChoose a group for the least squares fit\n");
190 get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname);
192 gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
195 printf("\nChoose a group for the covariance analysis\n");
196 get_index(atoms,ndxfile,1,&natoms,&index,&ananame);
200 snew(w_rls,atoms->nr);
201 for(i=0; (i<nfit); i++) {
202 w_rls[ifit[i]]=atoms->atom[ifit[i]].m;
204 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]]!=w_rls[ifit[i-1]]);
209 for(i=0; (i<natoms); i++)
211 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
213 bDiffMass2 = bDiffMass2 || (sqrtm[i]!=sqrtm[i-1]);
218 if (bFit && bDiffMass1 && !bDiffMass2) {
219 bDiffMass1 = natoms != nfit;
221 for (i=0; (i<natoms) && !bDiffMass1; i++)
222 bDiffMass1 = index[i] != ifit[i];
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++)
233 /* Prepare reference frame */
235 gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box);
236 gmx_rmpbc(gpbc,atoms->nr,box,xref);
239 reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
244 if (sqrt(GMX_LARGE_INT_MAX)<ndim) {
245 gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
249 fprintf(stderr,"Calculating the average structure ...\n");
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);
256 /* calculate x: a fitted struture of the selected atoms */
258 gmx_rmpbc(gpbc,nat,box,xread);
260 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
261 do_fit(nat,w_rls,xref,xread);
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));
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];
274 write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
275 atoms,xread,NULL,epbcNONE,zerobox,natoms,index);
278 fprintf(stderr,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim,(int)ndim);
280 nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
285 /* calculate x: a (fitted) structure of the selected atoms */
287 gmx_rmpbc(gpbc,nat,box,xread);
289 reset_x(nfit,ifit,nat,NULL,xread,w_rls);
290 do_fit(nat,w_rls,xref,xread);
293 for (i=0; i<natoms; i++)
294 rvec_sub(xread[index[i]],xref[index[i]],x[i]);
296 for (i=0; i<natoms; i++)
297 rvec_sub(xread[index[i]],xav[i],x[i]);
299 for (j=0; j<natoms; j++) {
300 for (dj=0; dj<DIM; dj++) {
303 for (i=j; i<natoms; i++) {
306 mat[l+d] += x[i][d]*xj;
310 } while (read_next_x(oenv,status,&t,nat,xread,box) &&
311 (bRef || nframes < nframes0));
313 gmx_rmpbc_done(gpbc);
315 fprintf(stderr,"Read %d frames\n",nframes);
318 /* copy the reference structure to the ouput array x */
320 for (i=0; i<natoms; i++)
321 copy_rvec(xref[index[i]],xproj[i]);
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];
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];
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 " : "");
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]);
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)
366 if (mat2[j][j] > max)
371 for(i=0; i<ndim; i++)
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");
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);
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++) {
396 mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
397 if (mat2[j][i] < min)
399 if (mat2[j][j] > max)
401 mat2[i][j] = mat2[j][i];
405 for(i=0; i<ndim/DIM; i++)
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");
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);
417 for (i=0; i<ndim/DIM; i++)
423 /* call diagonalization routine */
425 fprintf(stderr,"\nDiagonalizing ...\n");
430 memcpy(tmp,mat,ndim*ndim*sizeof(real));
431 eigensolver(tmp,ndim,0,ndim,eigval,mat);
434 /* now write the output */
437 for(i=0; i<ndim; i++)
439 fprintf(stderr,"\nSum of the eigenvalues: %g (%snm^2)\n",
441 if (fabs(trace-sum)>0.01*trace)
442 fprintf(stderr,"\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
444 fprintf(stderr,"\nWriting eigenvalues to %s\n",eigvalfile);
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]);
455 if (nframes-1 < ndim)
461 /* misuse lambda: 0/1 mass weighted analysis no/yes */
463 WriteXref = eWXR_YES;
464 for(i=0; i<nfit; i++)
465 copy_rvec(xref[ifit[i]],x[i]);
469 /* misuse lambda: -1 for no fit */
470 WriteXref = eWXR_NOFIT;
473 write_eigenvectors(eigvecfile,natoms,mat,TRUE,1,end,
474 WriteXref,x,bDiffMass1,xproj,bM,eigval);
476 out = ffopen(logfile,"w");
479 gmx_ctime_r(&now,timebuf,STRLEN);
480 fprintf(out,"Covariance analysis log, written %s\n",timebuf);
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);
486 pcwd=getcwd(str,STRLEN);
490 gmx_fatal(FARGS,"Current working directory is undefined");
493 fprintf(out,"Working directory: %s\n\n",str);
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));
498 fprintf(out,"Read reference structure for fit from %s\n",fitfile);
500 fprintf(out,"Read index groups from %s\n",ndxfile);
503 fprintf(out,"Analysis group is '%s' (%d atoms)\n",ananame,natoms);
505 fprintf(out,"Fit group is '%s' (%d atoms)\n",fitname,nfit);
507 fprintf(out,"No fit was used\n");
508 fprintf(out,"Analysis is %smass weighted\n", bDiffMass2 ? "":"non-");
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",
514 fprintf(out,"Trace of the covariance matrix after diagonalizing: %g\n\n",
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);
525 fprintf(stderr,"Wrote the log to %s\n",logfile);