changed reading hint
[gromacs/adressmacs.git] / src / tools / g_covar.c
blobe7de474a09bcfc7ec222dcf0418c8fe3279054c3
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
29 static char *SRCID_g_covar_c = "$Id$";
31 #include <math.h>
32 #include <string.h>
33 #include <unistd.h>
34 #include <time.h>
35 #include "statutil.h"
36 #include "sysstuff.h"
37 #include "typedefs.h"
38 #include "smalloc.h"
39 #include "macros.h"
40 #include "vec.h"
41 #include "pbc.h"
42 #include "copyrite.h"
43 #include "futil.h"
44 #include "statutil.h"
45 #include "rdgroup.h"
46 #include "confio.h"
47 #include "trnio.h"
48 #include "mshift.h"
49 #include "xvgr.h"
50 #include "do_fit.h"
51 #include "rmpbc.h"
52 #include "txtdump.h"
53 #include "eigio.h"
54 #include "ql77.h"
56 int main(int argc,char *argv[])
58 static char *desc[] = {
59 "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
60 "covariance matrix.",
61 "All structures are fitted to the structure in the structure file.",
62 "When this is not a run input file periodicity will not be taken into",
63 "account. When the fit and analysis groups are identical and the analysis",
64 "is non mass-weighted, the fit will also be non mass-weighted.[PAR]",
65 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
66 "When the same atoms are used for the fit and the covariance analysis,",
67 "the reference structure is written first with t=-1.",
68 "The average structure is written with t=0, the eigenvectors",
69 "are written as frames with the eigenvector number as timestamp.",
70 "The eigenvectors can be analyzed with [TT]g_anaeig[tt]."
72 static bool bFit=TRUE,bM=FALSE;
73 static int end=-1;
74 t_pargs pa[] = {
75 { "-fit", FALSE, etBOOL, {&bFit},
76 "Fit to a reference structure"},
77 { "-mwa", FALSE, etBOOL, {&bM},
78 "Mass-weighted covariance analysis"},
79 { "-last", FALSE, etINT, {&end},
80 "Last eigenvector to write away (-1 is till the last)" }
82 FILE *out;
83 int status,trjout;
84 t_topology top;
85 t_atoms *atoms;
86 rvec *x,*xread,*xref,*xav;
87 matrix box,zerobox;
88 real t,tstart,tend,*mat,dev,trace,sum,*eigval,inv_nframes;
89 real xj,*sqrtm,*w_rls=NULL;
90 int ntopatoms,step;
91 int natoms,nat,ndim,count,nframes;
92 int WriteXref;
93 char *fitfile,*trxfile,*ndxfile;
94 char *eigvalfile,*eigvecfile,*averfile,*logfile;
95 char str[STRLEN],*fitname,*ananame;
96 int i,j,k,l,d,dj,nfit;
97 atom_id *index,*all_at,*ifit;
98 bool bDiffMass1,bDiffMass2;
99 time_t now;
100 t_filenm fnm[] = {
101 { efTRX, "-f", NULL, ffREAD },
102 { efTPS, NULL, NULL, ffREAD },
103 { efNDX, NULL, NULL, ffOPTRD },
104 { efXVG, NULL, "eigenval", ffWRITE },
105 { efTRN, "-v", "eigenvec", ffWRITE },
106 { efSTO, "-av", "average.pdb", ffWRITE },
107 { efLOG, NULL, "covar", ffWRITE }
109 #define NFILE asize(fnm)
111 CopyRight(stderr,argv[0]);
112 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_SET_NPRI,TRUE,
113 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
115 clear_mat(zerobox);
117 fitfile = ftp2fn(efTPS,NFILE,fnm);
118 trxfile = ftp2fn(efTRX,NFILE,fnm);
119 ndxfile = ftp2fn_null(efNDX,NFILE,fnm);
120 eigvalfile = ftp2fn(efXVG,NFILE,fnm);
121 eigvecfile = ftp2fn(efTRN,NFILE,fnm);
122 averfile = ftp2fn(efSTO,NFILE,fnm);
123 logfile = ftp2fn(efLOG,NFILE,fnm);
125 read_tps_conf(fitfile,str,&top,&xref,NULL,box,TRUE);
126 atoms=&top.atoms;
128 if (bFit) {
129 printf("\nChoose a group for the least squares fit\n");
130 get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname);
131 if (nfit < 3)
132 fatal_error(0,"Need >= 3 points to fit!\n");
133 } else
134 nfit=0;
135 printf("\nChoose a group for the covariance analysis\n");
136 get_index(atoms,ndxfile,1,&natoms,&index,&ananame);
138 bDiffMass1=FALSE;
139 if (bFit) {
140 snew(w_rls,atoms->nr);
141 for(i=0; (i<nfit); i++) {
142 w_rls[ifit[i]]=atoms->atom[ifit[i]].m;
143 if (i)
144 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]]!=w_rls[ifit[i-1]]);
147 bDiffMass2=FALSE;
148 snew(sqrtm,natoms);
149 for(i=0; (i<natoms); i++)
150 if (bM) {
151 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
152 if (i)
153 bDiffMass2 = bDiffMass2 || (sqrtm[i]!=sqrtm[i-1]);
155 else
156 sqrtm[i]=1.0;
158 if (bFit && bDiffMass1 && !bDiffMass2) {
159 bDiffMass1 = natoms != nfit;
160 i=0;
161 for (i=0; (i<natoms) && !bDiffMass1; i++)
162 bDiffMass1 = index[i] != ifit[i];
163 if (!bDiffMass1) {
164 fprintf(stderr,"\n"
165 "Note: the fit and analysis group are identical,\n"
166 " while the fit is mass weighted and the analysis is not.\n"
167 " Making the fit non mass weighted.\n\n");
168 for(i=0; (i<nfit); i++)
169 w_rls[ifit[i]]=1.0;
173 snew(all_at,atoms->nr);
174 for(i=0; (i<atoms->nr); i++)
175 all_at[i]=i;
177 /* Prepare reference frame */
178 rm_pbc(&(top.idef),atoms->nr,box,xref,xref);
179 if (bFit)
180 reset_x(nfit,ifit,atoms->nr,all_at,xref,w_rls);
182 snew(x,natoms);
183 snew(xav,natoms);
184 ndim=natoms*DIM;
185 snew(mat,ndim*ndim);
187 fprintf(stderr,"Constructing covariance matrix (%dx%d)...\n",ndim,ndim);
189 nframes=0;
190 nat=read_first_x(&status,trxfile,&t,&xread,box);
191 tstart = t;
192 do {
193 nframes++;
194 tend = t;
195 /* calculate x: a fitted struture of the selected atoms */
196 rm_pbc(&(top.idef),nat,box,xread,xread);
197 if (bFit) {
198 reset_x(nfit,ifit,nat,all_at,xread,w_rls);
199 do_fit(nat,w_rls,xref,xread);
201 for (i=0; i<natoms; i++)
202 copy_rvec(xread[index[i]],x[i]);
204 for (j=0; j<natoms; j++) {
205 /* calculate average structure */
206 rvec_inc(xav[j],x[j]);
207 /* calculate cross product matrix */
208 for (dj=0; dj<DIM; dj++) {
209 k=ndim*(DIM*j+dj);
210 xj=x[j][dj];
211 for (i=j; i<natoms; i++) {
212 l=k+DIM*i;
213 for(d=0; d<DIM; d++)
214 mat[l+d] += x[i][d]*xj;
218 } while (read_next_x(status,&t,nat,xread,box));
219 close_trj(status);
221 /* calculate the mass-weighted covariance matrix */
222 inv_nframes=1.0/nframes;
223 for (i=0; i<natoms; i++) {
224 svmul(inv_nframes,xav[i],xav[i]);
225 copy_rvec(xav[i],xread[index[i]]);
228 write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
229 atoms,xread,NULL,NULL,natoms,index);
231 for (j=0; j<natoms; j++)
232 for (dj=0; dj<DIM; dj++)
233 for (i=j; i<natoms; i++) {
234 k=ndim*(DIM*j+dj)+DIM*i;
235 for (d=0; d<DIM; d++) {
236 mat[k+d]=(mat[k+d]*inv_nframes-xav[i][d]*xav[j][dj])
237 *sqrtm[i]*sqrtm[j];
240 for (j=0; j<ndim; j++)
241 for (i=j; i<ndim; i++)
242 mat[ndim*i+j]=mat[ndim*j+i];
244 trace=0;
245 for(i=0; i<ndim; i++)
246 trace+=mat[i*ndim+i];
247 fprintf(stderr,"\nTrace of the covariance matrix: %g (%snm^2)\n",
248 trace,bM ? "amu " : "");
250 if (debug) {
251 fprintf(stderr,"Dumping the covariance matrix to covmat.dat\n");
252 out = ffopen("covmat.dat","w");
253 for (j=0; j<ndim; j++) {
254 for (i=0; i<ndim; i++)
255 fprintf(out," %g",mat[ndim*j+i]);
256 fprintf(out,"\n");
260 /* call diagonalization routine */
262 fprintf(stderr,"\nDiagonalizing...\n");
263 fflush(stderr);
265 snew(eigval,ndim);
266 ql77 (ndim,mat,eigval);
268 /* now write the output */
270 sum=0;
271 for(i=0; i<ndim; i++)
272 sum+=eigval[i];
273 fprintf(stderr,"\nSum of the eigenvalues: %g (%snm^2)\n",
274 sum,bM ? "amu " : "");
275 if (fabs(trace-sum)>0.01*trace)
276 fprintf(stderr,"\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
278 fprintf(stderr,"\nWriting eigenvalues to %s\n",eigvalfile);
280 sprintf(str,"(%snm\\S2\\N)",bM ? "amu " : "");
281 out=xvgropen(eigvalfile,
282 "Eigenvalues of the covariance matrix",
283 "Eigenvector index",str);
284 for (i=0; (i<ndim); i++)
285 fprintf (out,"%10d %g\n",i+1,eigval[ndim-1-i]);
286 fclose(out);
288 if (end==-1) {
289 if (nframes-1 < ndim)
290 end=nframes-1;
291 else
292 end=ndim;
294 if (bFit) {
295 /* misuse lambda: 0/1 mass weighted analysis no/yes */
296 if (nfit==natoms) {
297 WriteXref = eWXR_YES;
298 for(i=0; i<nfit; i++)
299 copy_rvec(xref[ifit[i]],x[i]);
300 } else
301 WriteXref = eWXR_NO;
302 } else {
303 /* misuse lambda: -1 for no fit */
304 WriteXref = eWXR_NOFIT;
307 write_eigenvectors(eigvecfile,natoms,mat,TRUE,1,end,
308 WriteXref,x,bDiffMass1,xav,bM);
310 out = ffopen(logfile,"w");
312 now = time(NULL);
313 fprintf(out,"Covariance analysis log, written %s\n",
314 ctime(&now));
315 fprintf(out,"Program: %s\n",argv[0]);
316 getcwd(str,STRLEN);
317 fprintf(out,"Working directory: %s\n\n",str);
319 fprintf(out,"Read %d frames from %s (time %g to %g)\n",nframes,trxfile,
320 tstart,tend);
321 if (bFit)
322 fprintf(out,"Read reference structure for fit from %s\n",fitfile);
323 if (ndxfile)
324 fprintf(out,"Read index groups from %s\n",ndxfile);
325 fprintf(out,"\n");
327 fprintf(out,"Analysis group is '%s' (%d atoms)\n",ananame,natoms);
328 if (bFit)
329 fprintf(out,"Fit group is '%s' (%d atoms)\n",fitname,nfit);
330 else
331 fprintf(out,"No fit was used\n");
332 fprintf(out,"Analysis is %smass weighted\n", bDiffMass2 ? "":"non-");
333 if (bFit)
334 fprintf(out,"Fit is %smass weighted\n", bDiffMass1 ? "":"non-");
335 fprintf(out,"Diagonalized the %dx%d covariance matrix\n",ndim,ndim);
336 fprintf(out,"Trace of the covariance matrix before diagonalizing: %g\n",
337 trace);
338 fprintf(out,"Trace of the covariance matrix after diagonalizing: %g\n\n",
339 sum);
341 fprintf(out,"Wrote %d eigenvalues to %s\n",ndim,eigvalfile);
342 if (WriteXref == eWXR_YES)
343 fprintf(out,"Wrote reference structure to %s\n",eigvecfile);
344 fprintf(out,"Wrote average structure to %s and %s\n",averfile,eigvecfile);
345 fprintf(out,"Wrote eigenvectors %d to %d to %s\n",1,end,eigvecfile);
347 fclose(out);
349 fprintf(stderr,"Wrote the log to %s\n",logfile);
351 thanx(stdout);
353 return 0;