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
38 #include "gmx_header_config.h"
43 #ifdef HAVE_SYS_TIME_H
48 #ifdef GMX_NATIVE_WINDOWS
73 #include "eigensolver.h"
78 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
80 gmx_ctime_r(const time_t *clock
,char *buf
, int n
);
83 int gmx_covar(int argc
,char *argv
[])
85 const char *desc
[] = {
86 "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
88 "All structures are fitted to the structure in the structure file.",
89 "When this is not a run input file periodicity will not be taken into",
90 "account. When the fit and analysis groups are identical and the analysis",
91 "is non mass-weighted, the fit will also be non mass-weighted.",
93 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
94 "When the same atoms are used for the fit and the covariance analysis,",
95 "the reference structure for the fit is written first with t=-1.",
96 "The average (or reference when [TT]-ref[tt] is used) structure is",
97 "written with t=0, the eigenvectors",
98 "are written as frames with the eigenvector number as timestamp.",
100 "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
102 "Option [TT]-ascii[tt] writes the whole covariance matrix to",
103 "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
105 "Option [TT]-xpm[tt] writes the whole covariance matrix to an [TT].xpm[tt] file.",
107 "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [TT].xpm[tt] file,",
108 "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
111 "Note that the diagonalization of a matrix requires memory and time",
112 "that will increase at least as fast as than the square of the number",
113 "of atoms involved. It is easy to run out of memory, in which",
114 "case this tool will probably exit with a 'Segmentation fault'. You",
115 "should consider carefully whether a reduced set of atoms will meet",
116 "your needs for lower costs."
118 static gmx_bool bFit
=TRUE
,bRef
=FALSE
,bM
=FALSE
,bPBC
=TRUE
;
121 { "-fit", FALSE
, etBOOL
, {&bFit
},
122 "Fit to a reference structure"},
123 { "-ref", FALSE
, etBOOL
, {&bRef
},
124 "Use the deviation from the conformation in the structure file instead of from the average" },
125 { "-mwa", FALSE
, etBOOL
, {&bM
},
126 "Mass-weighted covariance analysis"},
127 { "-last", FALSE
, etINT
, {&end
},
128 "Last eigenvector to write away (-1 is till the last)" },
129 { "-pbc", FALSE
, etBOOL
, {&bPBC
},
130 "Apply corrections for periodic boundary conditions" }
138 rvec
*x
,*xread
,*xref
,*xav
,*xproj
;
140 real
*sqrtm
,*mat
,*eigenvalues
,sum
,trace
,inv_nframes
;
141 real t
,tstart
,tend
,**mat2
;
145 int natoms
,nat
,count
,nframes0
,nframes
,nlevels
;
146 gmx_large_int_t ndim
,i
,j
,k
,l
;
148 const char *fitfile
,*trxfile
,*ndxfile
;
149 const char *eigvalfile
,*eigvecfile
,*averfile
,*logfile
;
150 const char *asciifile
,*xpmfile
,*xpmafile
;
151 char str
[STRLEN
],*fitname
,*ananame
,*pcwd
;
153 atom_id
*index
,*ifit
;
154 gmx_bool bDiffMass1
,bDiffMass2
;
156 char timebuf
[STRLEN
];
160 gmx_rmpbc_t gpbc
=NULL
;
163 { efTRX
, "-f", NULL
, ffREAD
},
164 { efTPS
, NULL
, NULL
, ffREAD
},
165 { efNDX
, NULL
, NULL
, ffOPTRD
},
166 { efXVG
, NULL
, "eigenval", ffWRITE
},
167 { efTRN
, "-v", "eigenvec", ffWRITE
},
168 { efSTO
, "-av", "average.pdb", ffWRITE
},
169 { efLOG
, NULL
, "covar", ffWRITE
},
170 { efDAT
, "-ascii","covar", ffOPTWR
},
171 { efXPM
, "-xpm","covar", ffOPTWR
},
172 { efXPM
, "-xpma","covara", ffOPTWR
}
174 #define NFILE asize(fnm)
176 CopyRight(stderr
,argv
[0]);
177 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_BE_NICE
,
178 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
182 fitfile
= ftp2fn(efTPS
,NFILE
,fnm
);
183 trxfile
= ftp2fn(efTRX
,NFILE
,fnm
);
184 ndxfile
= ftp2fn_null(efNDX
,NFILE
,fnm
);
185 eigvalfile
= ftp2fn(efXVG
,NFILE
,fnm
);
186 eigvecfile
= ftp2fn(efTRN
,NFILE
,fnm
);
187 averfile
= ftp2fn(efSTO
,NFILE
,fnm
);
188 logfile
= ftp2fn(efLOG
,NFILE
,fnm
);
189 asciifile
= opt2fn_null("-ascii",NFILE
,fnm
);
190 xpmfile
= opt2fn_null("-xpm",NFILE
,fnm
);
191 xpmafile
= opt2fn_null("-xpma",NFILE
,fnm
);
193 read_tps_conf(fitfile
,str
,&top
,&ePBC
,&xref
,NULL
,box
,TRUE
);
197 printf("\nChoose a group for the least squares fit\n");
198 get_index(atoms
,ndxfile
,1,&nfit
,&ifit
,&fitname
);
200 gmx_fatal(FARGS
,"Need >= 3 points to fit!\n");
203 printf("\nChoose a group for the covariance analysis\n");
204 get_index(atoms
,ndxfile
,1,&natoms
,&index
,&ananame
);
208 snew(w_rls
,atoms
->nr
);
209 for(i
=0; (i
<nfit
); i
++) {
210 w_rls
[ifit
[i
]]=atoms
->atom
[ifit
[i
]].m
;
212 bDiffMass1
= bDiffMass1
|| (w_rls
[ifit
[i
]]!=w_rls
[ifit
[i
-1]]);
217 for(i
=0; (i
<natoms
); i
++)
219 sqrtm
[i
]=sqrt(atoms
->atom
[index
[i
]].m
);
221 bDiffMass2
= bDiffMass2
|| (sqrtm
[i
]!=sqrtm
[i
-1]);
226 if (bFit
&& bDiffMass1
&& !bDiffMass2
) {
227 bDiffMass1
= natoms
!= nfit
;
229 for (i
=0; (i
<natoms
) && !bDiffMass1
; i
++)
230 bDiffMass1
= index
[i
] != ifit
[i
];
233 "Note: the fit and analysis group are identical,\n"
234 " while the fit is mass weighted and the analysis is not.\n"
235 " Making the fit non mass weighted.\n\n");
236 for(i
=0; (i
<nfit
); i
++)
241 /* Prepare reference frame */
243 gpbc
= gmx_rmpbc_init(&top
.idef
,ePBC
,atoms
->nr
,box
);
244 gmx_rmpbc(gpbc
,atoms
->nr
,box
,xref
);
247 reset_x(nfit
,ifit
,atoms
->nr
,NULL
,xref
,w_rls
);
252 if (sqrt(GMX_LARGE_INT_MAX
)<ndim
) {
253 gmx_fatal(FARGS
,"Number of degrees of freedoms to large for matrix.\n");
257 fprintf(stderr
,"Calculating the average structure ...\n");
259 nat
=read_first_x(oenv
,&status
,trxfile
,&t
,&xread
,box
);
260 if (nat
!= atoms
->nr
)
261 fprintf(stderr
,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms
,nat
);
264 /* calculate x: a fitted struture of the selected atoms */
266 gmx_rmpbc(gpbc
,nat
,box
,xread
);
268 reset_x(nfit
,ifit
,nat
,NULL
,xread
,w_rls
);
269 do_fit(nat
,w_rls
,xref
,xread
);
271 for (i
=0; i
<natoms
; i
++)
272 rvec_inc(xav
[i
],xread
[index
[i
]]);
273 } while (read_next_x(oenv
,status
,&t
,nat
,xread
,box
));
276 inv_nframes
= 1.0/nframes0
;
277 for(i
=0; i
<natoms
; i
++)
278 for(d
=0; d
<DIM
; d
++) {
279 xav
[i
][d
] *= inv_nframes
;
280 xread
[index
[i
]][d
] = xav
[i
][d
];
282 write_sto_conf_indexed(opt2fn("-av",NFILE
,fnm
),"Average structure",
283 atoms
,xread
,NULL
,epbcNONE
,zerobox
,natoms
,index
);
286 fprintf(stderr
,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim
,(int)ndim
);
288 nat
=read_first_x(oenv
,&status
,trxfile
,&t
,&xread
,box
);
293 /* calculate x: a (fitted) structure of the selected atoms */
295 gmx_rmpbc(gpbc
,nat
,box
,xread
);
297 reset_x(nfit
,ifit
,nat
,NULL
,xread
,w_rls
);
298 do_fit(nat
,w_rls
,xref
,xread
);
301 for (i
=0; i
<natoms
; i
++)
302 rvec_sub(xread
[index
[i
]],xref
[index
[i
]],x
[i
]);
304 for (i
=0; i
<natoms
; i
++)
305 rvec_sub(xread
[index
[i
]],xav
[i
],x
[i
]);
307 for (j
=0; j
<natoms
; j
++) {
308 for (dj
=0; dj
<DIM
; dj
++) {
311 for (i
=j
; i
<natoms
; i
++) {
314 mat
[l
+d
] += x
[i
][d
]*xj
;
318 } while (read_next_x(oenv
,status
,&t
,nat
,xread
,box
) &&
319 (bRef
|| nframes
< nframes0
));
321 gmx_rmpbc_done(gpbc
);
323 fprintf(stderr
,"Read %d frames\n",nframes
);
326 /* copy the reference structure to the ouput array x */
328 for (i
=0; i
<natoms
; i
++)
329 copy_rvec(xref
[index
[i
]],xproj
[i
]);
334 /* correct the covariance matrix for the mass */
335 inv_nframes
= 1.0/nframes
;
336 for (j
=0; j
<natoms
; j
++)
337 for (dj
=0; dj
<DIM
; dj
++)
338 for (i
=j
; i
<natoms
; i
++) {
339 k
= ndim
*(DIM
*j
+dj
)+DIM
*i
;
340 for (d
=0; d
<DIM
; d
++)
341 mat
[k
+d
] = mat
[k
+d
]*inv_nframes
*sqrtm
[i
]*sqrtm
[j
];
344 /* symmetrize the matrix */
345 for (j
=0; j
<ndim
; j
++)
346 for (i
=j
; i
<ndim
; i
++)
347 mat
[ndim
*i
+j
]=mat
[ndim
*j
+i
];
350 for(i
=0; i
<ndim
; i
++)
351 trace
+=mat
[i
*ndim
+i
];
352 fprintf(stderr
,"\nTrace of the covariance matrix: %g (%snm^2)\n",
353 trace
,bM
? "u " : "");
356 out
= ffopen(asciifile
,"w");
357 for (j
=0; j
<ndim
; j
++) {
358 for (i
=0; i
<ndim
; i
+=3)
359 fprintf(out
,"%g %g %g\n",
360 mat
[ndim
*j
+i
],mat
[ndim
*j
+i
+1],mat
[ndim
*j
+i
+2]);
369 for (j
=0; j
<ndim
; j
++) {
370 mat2
[j
] = &(mat
[ndim
*j
]);
371 for (i
=0; i
<=j
; i
++) {
372 if (mat2
[j
][i
] < min
)
374 if (mat2
[j
][j
] > max
)
379 for(i
=0; i
<ndim
; i
++)
381 rlo
.r
= 0; rlo
.g
= 0; rlo
.b
= 1;
382 rmi
.r
= 1; rmi
.g
= 1; rmi
.b
= 1;
383 rhi
.r
= 1; rhi
.g
= 0; rhi
.b
= 0;
384 out
= ffopen(xpmfile
,"w");
386 write_xpm3(out
,0,"Covariance",bM
? "u nm^2" : "nm^2",
387 "dim","dim",ndim
,ndim
,axis
,axis
,
388 mat2
,min
,0.0,max
,rlo
,rmi
,rhi
,&nlevels
);
398 for (i
=0; i
<ndim
/DIM
; i
++)
399 snew(mat2
[i
],ndim
/DIM
);
400 for (j
=0; j
<ndim
/DIM
; j
++) {
401 for (i
=0; i
<=j
; i
++) {
404 mat2
[j
][i
] += mat
[ndim
*(DIM
*j
+d
)+DIM
*i
+d
];
405 if (mat2
[j
][i
] < min
)
407 if (mat2
[j
][j
] > max
)
409 mat2
[i
][j
] = mat2
[j
][i
];
413 for(i
=0; i
<ndim
/DIM
; i
++)
415 rlo
.r
= 0; rlo
.g
= 0; rlo
.b
= 1;
416 rmi
.r
= 1; rmi
.g
= 1; rmi
.b
= 1;
417 rhi
.r
= 1; rhi
.g
= 0; rhi
.b
= 0;
418 out
= ffopen(xpmafile
,"w");
420 write_xpm3(out
,0,"Covariance",bM
? "u nm^2" : "nm^2",
421 "atom","atom",ndim
/DIM
,ndim
/DIM
,axis
,axis
,
422 mat2
,min
,0.0,max
,rlo
,rmi
,rhi
,&nlevels
);
425 for (i
=0; i
<ndim
/DIM
; i
++)
431 /* call diagonalization routine */
433 snew(eigenvalues
,ndim
);
434 snew(eigenvectors
,ndim
*ndim
);
436 memcpy(eigenvectors
,mat
,ndim
*ndim
*sizeof(real
));
437 fprintf(stderr
,"\nDiagonalizing ...\n");
439 eigensolver(eigenvectors
,ndim
,0,ndim
,eigenvalues
,mat
);
442 /* now write the output */
445 for(i
=0; i
<ndim
; i
++)
447 fprintf(stderr
,"\nSum of the eigenvalues: %g (%snm^2)\n",
449 if (fabs(trace
-sum
)>0.01*trace
)
450 fprintf(stderr
,"\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
452 fprintf(stderr
,"\nWriting eigenvalues to %s\n",eigvalfile
);
454 sprintf(str
,"(%snm\\S2\\N)",bM
? "u " : "");
455 out
=xvgropen(eigvalfile
,
456 "Eigenvalues of the covariance matrix",
457 "Eigenvector index",str
,oenv
);
458 for (i
=0; (i
<ndim
); i
++)
459 fprintf (out
,"%10d %g\n",(int)i
+1,eigenvalues
[ndim
-1-i
]);
463 if (nframes
-1 < ndim
)
469 /* misuse lambda: 0/1 mass weighted analysis no/yes */
471 WriteXref
= eWXR_YES
;
472 for(i
=0; i
<nfit
; i
++)
473 copy_rvec(xref
[ifit
[i
]],x
[i
]);
477 /* misuse lambda: -1 for no fit */
478 WriteXref
= eWXR_NOFIT
;
481 write_eigenvectors(eigvecfile
,natoms
,mat
,TRUE
,1,end
,
482 WriteXref
,x
,bDiffMass1
,xproj
,bM
,eigenvalues
);
484 out
= ffopen(logfile
,"w");
487 gmx_ctime_r(&now
,timebuf
,STRLEN
);
488 fprintf(out
,"Covariance analysis log, written %s\n",timebuf
);
490 fprintf(out
,"Program: %s\n",argv
[0]);
491 #ifdef GMX_NATIVE_WINDOWS
492 pcwd
=_getcwd(str
,STRLEN
);
494 pcwd
=getcwd(str
,STRLEN
);
498 gmx_fatal(FARGS
,"Current working directory is undefined");
501 fprintf(out
,"Working directory: %s\n\n",str
);
503 fprintf(out
,"Read %d frames from %s (time %g to %g %s)\n",nframes
,trxfile
,
504 output_env_conv_time(oenv
,tstart
),output_env_conv_time(oenv
,tend
),output_env_get_time_unit(oenv
));
506 fprintf(out
,"Read reference structure for fit from %s\n",fitfile
);
508 fprintf(out
,"Read index groups from %s\n",ndxfile
);
511 fprintf(out
,"Analysis group is '%s' (%d atoms)\n",ananame
,natoms
);
513 fprintf(out
,"Fit group is '%s' (%d atoms)\n",fitname
,nfit
);
515 fprintf(out
,"No fit was used\n");
516 fprintf(out
,"Analysis is %smass weighted\n", bDiffMass2
? "":"non-");
518 fprintf(out
,"Fit is %smass weighted\n", bDiffMass1
? "":"non-");
519 fprintf(out
,"Diagonalized the %dx%d covariance matrix\n",(int)ndim
,(int)ndim
);
520 fprintf(out
,"Trace of the covariance matrix before diagonalizing: %g\n",
522 fprintf(out
,"Trace of the covariance matrix after diagonalizing: %g\n\n",
525 fprintf(out
,"Wrote %d eigenvalues to %s\n",(int)ndim
,eigvalfile
);
526 if (WriteXref
== eWXR_YES
)
527 fprintf(out
,"Wrote reference structure to %s\n",eigvecfile
);
528 fprintf(out
,"Wrote average structure to %s and %s\n",averfile
,eigvecfile
);
529 fprintf(out
,"Wrote eigenvectors %d to %d to %s\n",1,end
,eigvecfile
);
533 fprintf(stderr
,"Wrote the log to %s\n",logfile
);