?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / g_anaeig.c
blob75fe2f8b97a4094981aaa89af3d2326fb83d6689
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
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
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)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_g_anaeig_c = "$Id$";
30
31 #include <math.h>
32 #include <string.h>
33 #include "statutil.h"
34 #include "sysstuff.h"
35 #include "typedefs.h"
36 #include "smalloc.h"
37 #include "macros.h"
38 #include "fatal.h"
39 #include "vec.h"
40 #include "pbc.h"
41 #include "copyrite.h"
42 #include "futil.h"
43 #include "statutil.h"
44 #include "rdgroup.h"
45 #include "pdbio.h"
46 #include "confio.h"
47 #include "tpxio.h"
48 #include "matio.h"
49 #include "mshift.h"
50 #include "xvgr.h"
51 #include "do_fit.h"
52 #include "rmpbc.h"
53 #include "txtdump.h"
54 #include "eigio.h"
55
56 char *proj_unit;
57
58 static real tick_spacing(real range,int minticks)
59 {
60 real sp;
61
62 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
63 while (range/sp < minticks-1)
64 sp = sp/2;
65
66 return sp;
67 }
68
69 static void write_xvgr_graphs(char *file,int ngraphs,
70 char *title,char *xlabel,char **ylabel,
71 int n,real *x, real **y,bool bZero)
72 {
73 FILE *out;
74 int g,i;
75 real min,max,xsp,ysp;
76
77 out=ffopen(file,"w");
78 for(g=0; g<ngraphs; g++) {
79 min=y[g][0];
80 max=y[g][0];
81 for(i=0; i<n; i++) {
82 if (y[g][i]<min) min=y[g][i];
83 if (y[g][i]>max) max=y[g][i];
84 }
85 if (bZero)
86 min=0;
87 else
88 min=min-0.1*(max-min);
89 max=max+0.1*(max-min);
90 xsp=tick_spacing(x[n-1]-x[0],4);
91 ysp=tick_spacing(max-min,3);
92 fprintf(out,"@ with g%d\n@ g%d on\n",ngraphs-1-g,ngraphs-1-g);
93 fprintf(out,"@ g%d autoscale type AUTO\n",ngraphs-1-g);
94 if (g==0)
95 fprintf(out,"@ title \"%s\"\n",title);
96 if (g==ngraphs-1)
97 fprintf(out,"@ xaxis label \"%s\"\n",xlabel);
98 else
99 fprintf(out,"@ xaxis ticklabel off\n");
100 fprintf(out,"@ world xmin %g\n",x[0]);
101 fprintf(out,"@ world xmax %g\n",x[n-1]);
102 fprintf(out,"@ world ymin %g\n",min);
103 fprintf(out,"@ world ymax %g\n",max);
104 fprintf(out,"@ view xmin 0.15\n");
105 fprintf(out,"@ view xmax 0.85\n");
106 fprintf(out,"@ view ymin %g\n",0.15+(ngraphs-1-g)*0.7/ngraphs);
107 fprintf(out,"@ view ymax %g\n",0.15+(ngraphs-g)*0.7/ngraphs);
108 fprintf(out,"@ yaxis label \"%s\"\n",ylabel[g]);
109 fprintf(out,"@ xaxis tick major %g\n",xsp);
110 fprintf(out,"@ xaxis tick minor %g\n",xsp/2);
111 fprintf(out,"@ xaxis ticklabel start type spec\n");
112 fprintf(out,"@ xaxis ticklabel start %g\n",ceil(min/xsp)*xsp);
113 fprintf(out,"@ yaxis tick major %g\n",ysp);
114 fprintf(out,"@ yaxis tick minor %g\n",ysp/2);
115 fprintf(out,"@ yaxis ticklabel start type spec\n");
116 fprintf(out,"@ yaxis ticklabel start %g\n",ceil(min/ysp)*ysp);
117 if ((min<0) && (max>0)) {
118 fprintf(out,"@ zeroxaxis bar on\n");
119 fprintf(out,"@ zeroxaxis bar linestyle 3\n");
120 }
121 for(i=0; i<n; i++)
122 fprintf(out,"%10.4f %10.5f\n",x[i],y[g][i]);
123 fprintf(out,"&\n");
124 }
125 fclose(out);
126 }
127
128 static void inprod_matrix(char *matfile,int natoms,
129 int nvec1,int *eignr1,rvec **eigvec1,
130 int nvec2,int *eignr2,rvec **eigvec2,
131 bool bSelect,int noutvec,int *outvec)
132 {
133 FILE *out;
134 real **mat;
135 int i,x1,x,y,nlevels;
136 int n1;
137 real inp,*t_x,*t_y,max;
138 t_rgb rlo,rhi;
139
140 if (bSelect)
141 n1 = noutvec;
142 else
143 n1 = nvec1;
144
145 fprintf(stderr,"Calculating inner-product matrix of %dx%d eigenvectors\n",
146 n1,nvec2);
147
148 snew(mat,n1);
149 for(x=0; x<n1; x++)
150 snew(mat[x],nvec2);
151
152 snew(t_x,n1);
153 max = 0;
154 for(x1=0; x1<n1; x1++) {
155 if (bSelect)
156 x = outvec[x1];
157 else
158 x = x1;
159 t_x[x1] = eignr1[x]+1;
160 fprintf(stderr," %d",eignr1[x]+1);
161 for(y=0; y<nvec2; y++) {
162 inp = 0;
163 for(i=0; i<natoms; i++)
164 inp += iprod(eigvec1[x][i],eigvec2[y][i]);
165 mat[x1][y] = fabs(inp);
166 if (mat[x1][y]>max)
167 max = mat[x1][y];
168 }
169 }
170 fprintf(stderr,"\n");
171 snew(t_y,nvec2);
172 for(i=0; i<nvec2; i++)
173 t_y[i] = eignr2[i]+1;
174 rlo.r = 1; rlo.g = 1; rlo.b = 1;
175 rhi.r = 0; rhi.g = 0; rhi.b = 0;
176 nlevels = 41;
177 out = ffopen(matfile,"w");
178 write_xpm(out,"Eigenvector inner-products","in.prod.","run 1","run 2",
179 n1,nvec2,t_x,t_y,mat,0.0,max,rlo,rhi,&nlevels);
180 fclose(out);
181 }
182
183 static void overlap(char *outfile,int natoms,
184 rvec **eigvec1,
185 int nvec2,int *eignr2,rvec **eigvec2,
186 int noutvec,int *outvec)
187 {
188 FILE *out;
189 int i,v,vec,x;
190 real overlap,inp;
191
192 fprintf(stderr,"Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
193 for(i=0; i<noutvec; i++)
194 fprintf(stderr,"%d ",outvec[i]+1);
195 fprintf(stderr,"\n");
196
197 out=xvgropen(outfile,"Subspace overlap",
198 "Eigenvectors of trajectory 2","Overlap");
199 fprintf(out,"@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
200 noutvec);
201 overlap=0;
202 for(x=0; x<nvec2; x++) {
203 for(v=0; v<noutvec; v++) {
204 vec=outvec[v];
205 inp=0;
206 for(i=0; i<natoms; i++)
207 inp+=iprod(eigvec1[vec][i],eigvec2[x][i]);
208 overlap+=sqr(inp);
209 }
210 fprintf(out,"%5d %5.3f\n",eignr2[x]+1,overlap/noutvec);
211 }
212
213 fclose(out);
214 }
215
216 static void project(char *trajfile,t_topology *top,matrix topbox,rvec *xtop,
217 char *projfile,char *twodplotfile,char *threedplotfile,
218 char *filterfile,int skip,
219 char *extremefile,bool bExtrAll,real extreme,int nextr,
220 t_atoms *atoms,int natoms,atom_id *index,
221 bool bFit,rvec *xref,int nfit,atom_id *ifit,real *w_rls,
222 real *sqrtm,rvec *xav,
223 int *eignr,rvec **eigvec,
224 int noutvec,int *outvec)
225 {
226 FILE *xvgrout=NULL;
227 int status,out=0,nat,i,j,d,v,vec,nfr,nframes=0,snew_size,frame;
228 int noutvec_extr,*imin,*imax;
229 atom_id *all_at;
230 matrix box;
231 rvec *xread,*x;
232 real t,inp,**inprod=NULL,min=0,max=0;
233 char str[STRLEN],str2[STRLEN],**ylabel,*c;
234
235 snew(x,natoms);
236
237 if (bExtrAll)
238 noutvec_extr=noutvec;
239 else
240 noutvec_extr=1;
241
242
243 if (trajfile) {
244 snew(inprod,noutvec+1);
245
246 if (filterfile) {
247 fprintf(stderr,"Writing a filtered trajectory to %s using eigenvectors\n",
248 filterfile);
249 for(i=0; i<noutvec; i++)
250 fprintf(stderr,"%d ",outvec[i]+1);
251 fprintf(stderr,"\n");
252 out=open_trx(filterfile,"w");
253 }
254 snew_size=0;
255 nfr=0;
256 nframes=0;
257 nat=read_first_x(&status,trajfile,&t,&xread,box);
258 if (nat>atoms->nr)
259 fatal_error(0,"the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)",nat,atoms->nr);
260 snew(all_at,nat);
261 for(i=0; (i<nat); i++)
262 all_at[i]=i;
263 do {
264 if (nfr % skip == 0) {
265 if (top)
266 rm_pbc(&(top->idef),nat,box,xread,xread);
267 if (nframes>=snew_size) {
268 snew_size+=100;
269 for(i=0; i<noutvec+1; i++)
270 srenew(inprod[i],snew_size);
271 }
272 inprod[noutvec][nframes]=t;
273 /* calculate x: a fitted struture of the selected atoms */
274 if (bFit && (xref==NULL)) {
275 reset_x(nfit,ifit,nat,all_at,xread,w_rls);
276 do_fit(atoms->nr,w_rls,xtop,xread);
277 }
278 for (i=0; i<natoms; i++)
279 copy_rvec(xread[index[i]],x[i]);
280 if (bFit && xref) {
281 reset_x(natoms,all_at,natoms,all_at,x,w_rls);
282 do_fit(natoms,w_rls,xref,x);
283 }
284
285 for(v=0; v<noutvec; v++) {
286 vec=outvec[v];
287 /* calculate (mass-weighted) projection */
288 inp=0;
289 for (i=0; i<natoms; i++) {
290 inp+=(eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
291 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
292 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
293 }
294 inprod[v][nframes]=inp;
295 }
296 if (filterfile) {
297 for(i=0; i<natoms; i++)
298 for(d=0; d<DIM; d++) {
299 /* misuse xread for output */
300 xread[index[i]][d] = xav[i][d];
301 for(v=0; v<noutvec; v++)
302 xread[index[i]][d] +=
303 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
304 }
305 write_trx(out,natoms,index,atoms,0,t,box,xread,NULL);
306 }
307 nframes++;
308 }
309 nfr++;
310 } while (read_next_x(status,&t,nat,xread,box));
311 close_trj(status);
312 sfree(x);
313 if (filterfile)
314 close_trx(out);
315 }
316 else
317 snew(xread,atoms->nr);
318
319
320 if (projfile) {
321 snew(ylabel,noutvec);
322 for(v=0; v<noutvec; v++) {
323 sprintf(str,"vec %d",eignr[outvec[v]]+1);
324 ylabel[v]=strdup(str);
325 }
326 sprintf(str,"projection on eigenvectors (%s)",proj_unit);
327 write_xvgr_graphs(projfile,noutvec,
328 str,"Time (ps)",
329 ylabel,nframes,inprod[noutvec],inprod,FALSE);
330 }
331
332 if (twodplotfile) {
333 sprintf(str,"projection on eigenvector %d (%s)",
334 eignr[outvec[0]]+1,proj_unit);
335 sprintf(str2,"projection on eigenvector %d (%s)",
336 eignr[outvec[noutvec-1]]+1,proj_unit);
337 xvgrout=xvgropen(twodplotfile,"2D projection of trajectory",str,str2);
338 for(i=0; i<nframes; i++)
339 fprintf(xvgrout,"%10.5f %10.5f\n",inprod[0][i],inprod[noutvec-1][i]);
340 fclose(xvgrout);
341 }
342
343 if (threedplotfile) {
344 t_atoms atoms;
345 rvec *x;
346 matrix box;
347 char *resnm,*atnm;
348
349 if (noutvec < 3)
350 fatal_error(0,"You have selected less than 3 eigenvectors");
351
352 sprintf(str,"3D proj. of traj. on eigenv. %d, %d and %d",
353 eignr[outvec[0]]+1,eignr[outvec[1]]+1,eignr[outvec[2]]+1);
354 init_t_atoms(&atoms,nframes,FALSE);
355 snew(x,nframes);
356 atnm=strdup("CA");
357 resnm=strdup("GLY");
358 for(i=0; i<nframes; i++) {
359 atoms.resname[i]=&resnm;
360 atoms.atomname[i]=&atnm;
361 atoms.atom[i].resnr=i;
362 x[i][XX]=inprod[0][i];
363 x[i][YY]=inprod[1][i];
364 x[i][ZZ]=inprod[2][i];
365 }
366 clear_mat(box);
367 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
368 write_sto_conf(threedplotfile,str,&atoms,x,NULL,box);
369 free_t_atoms(&atoms);
370 fclose(xvgrout);
371 }
372
373 if (extremefile) {
374 if (extreme==0) {
375 fprintf(stderr,"%11s %17s %17s\n","eigenvector","Minimum","Maximum");
376 fprintf(stderr,
377 "%11s %10s %10s %10s %10s\n","","value","time","value","time");
378 snew(imin,noutvec_extr);
379 snew(imax,noutvec_extr);
380 for(v=0; v<noutvec_extr; v++) {
381 for(i=0; i<nframes; i++) {
382 if (inprod[v][i]<inprod[v][imin[v]])
383 imin[v]=i;
384 if (inprod[v][i]>inprod[v][imax[v]])
385 imax[v]=i;
386 }
387 min=inprod[v][imin[v]];
388 max=inprod[v][imax[v]];
389 fprintf(stderr,"%7d %10.6f %10.1f %10.6f %10.1f\n",
390 eignr[outvec[v]]+1,
391 min,inprod[noutvec][imin[v]],max,inprod[noutvec][imax[v]]);
392 }
393 }
394 else {
395 min=-extreme;
396 max=+extreme;
397 }
398 /* build format string for filename: */
399 strcpy(str,extremefile);/* copy filename */
400 c=strrchr(str,'.'); /* find where extention begins */
401 strcpy(str2,c); /* get extention */
402 sprintf(c,"%%d%s",str2); /* append '%s' and extention to filename */
403 for(v=0; v<noutvec_extr; v++) {
404 /* make filename using format string */
405 if (noutvec_extr==1)
406 strcpy(str2,extremefile);
407 else
408 sprintf(str2,str,eignr[outvec[v]]+1);
409 fprintf(stderr,"Writing %d frames along eigenvector %d to %s\n",
410 nextr,outvec[v]+1,str2);
411 out=open_trx(str2,"w");
412 for(frame=0; frame<nextr; frame++) {
413 if ((extreme==0) && (nextr<=3))
414 for(i=0; i<natoms; i++)
415 atoms->atom[index[i]].chain='A'+frame;
416 for(i=0; i<natoms; i++)
417 for(d=0; d<DIM; d++)
418 xread[index[i]][d] =
419 (xav[i][d] + (min*(nextr-frame-1)+max*frame)/(nextr-1)
420 *eigvec[outvec[v]][i][d]/sqrtm[i]);
421 write_trx(out,natoms,index,atoms,0,frame,topbox,xread,NULL);
422 }
423 close_trx(out);
424 }
425 }
426 fprintf(stderr,"\n");
427 }
428
429 static void components(char *outfile,int natoms,real *sqrtm,
430 int *eignr,rvec **eigvec,
431 int noutvec,int *outvec)
432 {
433 int g,v,i;
434 real *x,**y;
435 char str[STRLEN],**ylabel;
436
437 fprintf(stderr,"Writing atom displacements to %s\n",outfile);
438
439 snew(ylabel,noutvec);
440 snew(y,noutvec);
441 snew(x,natoms);
442 for(i=0; i<natoms; i++)
443 x[i]=i+1;
444 for(g=0; g<noutvec; g++) {
445 v=outvec[g];
446 sprintf(str,"vec %d",eignr[v]+1);
447 ylabel[g]=strdup(str);
448 snew(y[g],natoms);
449 for(i=0; i<natoms; i++)
450 y[g][i] = norm(eigvec[v][i])/sqrtm[i];
451 }
452 write_xvgr_graphs(outfile,noutvec,"Atom displacements","Atom number",
453 ylabel,natoms,x,y,TRUE);
454 fprintf(stderr,"\n");
455 }
456
457 int main(int argc,char *argv[])
458 {
459 static char *desc[] = {
460 "[TT]g_anaeig[tt] analyzes eigenvectors.",
461 "The eigenvectors can be of a covariance matrix ([TT]g_covar[tt])",
462 "or of a Normal Modes anaysis ([TT]g_nmeig[tt]).[PAR]",
463 "When a trajectory is projected on eigenvectors,",
464 "all structures are fitted to the structure in the eigenvector file,",
465 "if present, otherwise to the structure in the structure file.",
466 "When no run input file is supplied, periodicity will not be taken into",
467 "account.",
468 "Most analyses are done on eigenvectors [TT]-first[tt] to [TT]-last[tt],",
469 "but when [TT]-first[tt] is set to -1 you will be prompted for a",
470 "selection.[PAR]",
471 "[TT]-disp[tt]: plot all atom displacements of eigenvectors",
472 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
473 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
474 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
475 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
476 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
477 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
478 "three selected eigenvectors.[PAR]",
479 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
480 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
481 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
482 "on the average structure and interpolate [TT]-nframes[tt] frames between",
483 "them, or set your own extremes with [TT]-max[tt]. The eigenvector",
484 "[TT]-first[tt] will be written unless [TT]-first[tt] and [TT]-last[tt]",
485 "have been set explicitly, in which case all eigenvectors will be written",
486 "to separate files.",
487 "Chain identifiers will be added when",
488 "writing a [TT].pdb[tt] file with two or three structures",
489 "(you can use [TT]rasmol -nmrpdb[tt] to view such a pdb file).[PAR]",
490 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
491 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
492 "in file [TT]-v[tt].[PAR]",
493 "[TT]-inpr[tt]: calculate a matrix of inner-products between eigenvectors",
494 "in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors of the first file",
495 "will be used unless [TT]-first[tt] and [TT]-last[tt] have been set",
496 "explicitly."
497 };
498 static int first=1,last=8,skip=1,nextr=2;
499 static real max=0.0;
500 t_pargs pa[] = {
501 { "-first", FALSE, etINT, {&first},
502 "First eigenvector for analysis (-1 is select)" },
503 { "-last", FALSE, etINT, {&last},
504 "Last eigenvector for analysis (-1 is till the last)" },
505 { "-skip", FALSE, etINT, {&skip},
506 "Only analyse every nr-th frame" },
507 { "-max", FALSE, etREAL, {&max},
508 "Maximum for projection of the eigenvector on the average structure, max=0 gives the extremes" },
509 { "-nframes", FALSE, etINT, {&nextr},
510 "Number of frames for the extremes output" }
511 };
512 #define NPA asize(pa)
513
514 FILE *out;
515 int status,trjout;
516 t_topology top;
517 t_atoms *atoms=NULL;
518 rvec *xtop,*xref1,*xref2;
519 bool bDMR1,bDMA1,bDMR2,bDMA2;
520 int nvec1,nvec2,*eignr1=NULL,*eignr2=NULL;
521 rvec *x,*xread,*xav1,*xav2,**eigvec1=NULL,**eigvec2=NULL;
522 matrix topbox;
523 real xid,totmass,*sqrtm,*w_rls,t,lambda;
524 int natoms,step;
525 char *grpname,*indexfile,title[STRLEN];
526 int i,j,d;
527 int nout,*iout,noutvec,*outvec,nfit;
528 atom_id *index,*ifit;
529 char *Vec2File,*topfile,*CompFile;
530 char *ProjOnVecFile,*TwoDPlotFile,*ThreeDPlotFile;
531 char *FilterFile,*ExtremeFile;
532 char *OverlapFile,*InpMatFile;
533 bool bFit1,bFit2,bM,bIndex,bTPS,bTop,bVec2,bProj;
534 bool bFirstToLast,bFirstLastSet,bTraj;
535
536 t_filenm fnm[] = {
537 { efTRN, "-v", "eigenvec", ffREAD },
538 { efTRN, "-v2", "eigenvec2", ffOPTRD },
539 { efTRX, "-f", NULL, ffOPTRD },
540 { efTPS, NULL, NULL, ffOPTRD },
541 { efNDX, NULL, NULL, ffOPTRD },
542 { efXVG, "-disp", "eigdisp", ffOPTWR },
543 { efXVG, "-proj", "proj", ffOPTWR },
544 { efXVG, "-2d", "2dproj", ffOPTWR },
545 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
546 { efTRX, "-filt", "filtered", ffOPTWR },
547 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
548 { efXVG, "-over", "overlap", ffOPTWR },
549 { efXPM, "-inpr", "inprod", ffOPTWR }
550 };
551 #define NFILE asize(fnm)
552
553 CopyRight(stderr,argv[0]);
554 parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,
555 NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
556
557 indexfile=ftp2fn_null(efNDX,NFILE,fnm);
558
559 Vec2File = opt2fn_null("-v2",NFILE,fnm);
560 topfile = ftp2fn(efTPS,NFILE,fnm);
561 CompFile = opt2fn_null("-disp",NFILE,fnm);
562 ProjOnVecFile = opt2fn_null("-proj",NFILE,fnm);
563 TwoDPlotFile = opt2fn_null("-2d",NFILE,fnm);
564 ThreeDPlotFile = opt2fn_null("-3d",NFILE,fnm);
565 FilterFile = opt2fn_null("-filt",NFILE,fnm);
566 ExtremeFile = opt2fn_null("-extr",NFILE,fnm);
567 OverlapFile = opt2fn_null("-over",NFILE,fnm);
568 InpMatFile = ftp2fn_null(efXPM,NFILE,fnm);
569 bTop = fn2bTPX(topfile);
570 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
571 || FilterFile || ExtremeFile;
572 bFirstLastSet =
573 opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
574 bFirstToLast = CompFile || ProjOnVecFile || FilterFile || OverlapFile ||
575 ((ExtremeFile || InpMatFile) && bFirstLastSet);
576 bVec2 = Vec2File || OverlapFile || InpMatFile;
577 bM = CompFile || bProj;
578 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
579 || TwoDPlotFile || ThreeDPlotFile;
580 bIndex = bM || bProj;
581 bTPS = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
582 FilterFile || (bIndex && indexfile);
583
584 read_eigenvectors(opt2fn("-v",NFILE,fnm),&natoms,&bFit1,
585 &xref1,&bDMR1,&xav1,&bDMA1,&nvec1,&eignr1,&eigvec1);
586 if (bVec2) {
587 read_eigenvectors(Vec2File,&i,&bFit2,
588 &xref2,&bDMR2,&xav2,&bDMA2,&nvec2,&eignr2,&eigvec2);
589 if (i!=natoms)
590 fatal_error(0,"Dimensions in the eigenvector files don't match");
591 }
592
593 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
594 bM=FALSE;
595 if ((xref1==NULL) && (bM || bTraj))
596 bTPS=TRUE;
597
598 xtop=NULL;
599 nfit=0;
600 ifit=NULL;
601 w_rls=NULL;
602 if (!bTPS)
603 bTop=FALSE;
604 else {
605 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
606 title,&top,&xtop,NULL,topbox,bM);
607 atoms=&top.atoms;
608 rm_pbc(&(top.idef),atoms->nr,topbox,xtop,xtop);
609 /* Fitting is only needed when we need to read a trajectory */
610 if (bTraj) {
611 if ((xref1==NULL) || (bM && bDMR1)) {
612 if (bFit1) {
613 printf("\nNote: the structure in %s should be the same\n"
614 " as the one used for the fit in g_covar\n",topfile);
615 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
616 get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
617 snew(w_rls,atoms->nr);
618 for(i=0; (i<nfit); i++)
619 if (bM && bDMR1)
620 w_rls[ifit[i]]=atoms->atom[ifit[i]].m;
621 else
622 w_rls[ifit[i]]=1.0;
623 }
624 else {
625 /* make the fit index in xref instead of xtop */
626 nfit=natoms;
627 snew(ifit,natoms);
628 snew(w_rls,nfit);
629 for(i=0; (i<nfit); i++) {
630 ifit[i]=i;
631 w_rls[i]=atoms->atom[ifit[i]].m;
632 }
633 }
634 }
635 else {
636 /* make the fit non mass weighted on xref */
637 nfit=natoms;
638 snew(ifit,nfit);
639 snew(w_rls,nfit);
640 for(i=0; i<nfit; i++) {
641 ifit[i]=i;
642 w_rls[i]=1.0;
643 }
644 }
645 }
646 }
647
648 if (bIndex) {
649 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
650 get_index(atoms,indexfile,1,&i,&index,&grpname);
651 if (i!=natoms)
652 fatal_error(0,"you selected a group with %d elements instead of %d",
653 i,natoms);
654 printf("\n");
655 }
656
657 snew(sqrtm,natoms);
658 if (bM && bDMA1) {
659 proj_unit="amu\\S1/2\\Nnm";
660 for(i=0; (i<natoms); i++)
661 sqrtm[i]=sqrt(atoms->atom[index[i]].m);
662 } else {
663 proj_unit="nm";
664 for(i=0; (i<natoms); i++)
665 sqrtm[i]=1.0;
666 }
667
668 if (bVec2) {
669 t=0;
670 totmass=0;
671 for(i=0; (i<natoms); i++)
672 for(d=0;(d<DIM);d++) {
673 t+=sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
674 totmass+=sqr(sqrtm[i]);
675 }
676 fprintf(stderr,"RMSD (without fit) between the two average structures:"
677 " %.3f (nm)\n\n",sqrt(t/totmass));
678 }
679
680 if (last==-1)
681 last=natoms*DIM;
682 if (first>-1) {
683 if (bFirstToLast) {
684 /* make an index from first to last */
685 nout=last-first+1;
686 snew(iout,nout);
687 for(i=0; i<nout; i++)
688 iout[i]=first-1+i;
689 } else if (ThreeDPlotFile) {
690 /* make an index of first, first+1 and last */
691 nout=3;
692 snew(iout,nout);
693 iout[0]=first-1;
694 iout[1]=first;
695 iout[2]=last-1;
696 } else {
697 /* make an index of first and last */
698 nout=2;
699 snew(iout,nout);
700 iout[0]=first-1;
701 iout[1]=last-1;
702 }
703 }
704 else {
705 printf("Select eigenvectors for output, end your selection with 0\n");
706 nout=-1;
707 iout=NULL;
708 do {
709 nout++;
710 srenew(iout,nout+1);
711 scanf("%d",&iout[nout]);
712 iout[nout]--;
713 } while (iout[nout]>=0);
714 printf("\n");
715 }
716 /* make an index of the eigenvectors which are present */
717 snew(outvec,nout);
718 noutvec=0;
719 for(i=0; i<nout; i++) {
720 j=0;
721 while ((j<nvec1) && (eignr1[j]!=iout[i]))
722 j++;
723 if ((j<nvec1) && (eignr1[j]==iout[i])) {
724 outvec[noutvec]=j;
725 noutvec++;
726 }
727 }
728 fprintf(stderr,"%d eigenvectors selected for output",noutvec);
729 if (noutvec <= 100) {
730 fprintf(stderr,":");
731 for(j=0; j<noutvec; j++)
732 fprintf(stderr," %d",eignr1[outvec[j]]+1);
733 }
734 fprintf(stderr,"\n");
735
736 if (CompFile)
737 components(CompFile,natoms,sqrtm,eignr1,eigvec1,noutvec,outvec);
738
739 if (bProj)
740 project(bTraj ? opt2fn("-f",NFILE,fnm) : NULL,
741 bTop ? &top : NULL,topbox,xtop,
742 ProjOnVecFile,TwoDPlotFile,ThreeDPlotFile,FilterFile,skip,
743 ExtremeFile,bFirstLastSet,max,nextr,atoms,natoms,index,
744 bFit1,xref1,nfit,ifit,w_rls,
745 sqrtm,xav1,eignr1,eigvec1,noutvec,outvec);
746
747 if (OverlapFile)
748 overlap(OverlapFile,natoms,
749 eigvec1,nvec2,eignr2,eigvec2,noutvec,outvec);
750
751 if (InpMatFile)
752 inprod_matrix(InpMatFile,natoms,
753 nvec1,eignr1,eigvec1,nvec2,eignr2,eigvec2,
754 bFirstLastSet,noutvec,outvec);
755
756 if (!CompFile && !bProj && !OverlapFile && !InpMatFile)
757 fprintf(stderr,"\nIf you want some output,"
758 " set one (or two or ...) of the output file options\n");
759
760 return 0;
761 }
762
got merge a long time ago