?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / g_msd.c
bloba929d3c854a1fee36d453333e8f43cba7840d612
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_msd_c = "$Id$";
30
31
32 #include <string.h>
33 #include <ctype.h>
34 #include <math.h>
35 #include "sysstuff.h"
36 #include "smalloc.h"
37 #include "macros.h"
38 #include "statutil.h"
39 #include "maths.h"
40 #include "futil.h"
41 #include "rdgroup.h"
42 #include "copyrite.h"
43 #include "typedefs.h"
44 #include "xvgr.h"
45 #include "gstat.h"
46 #include "tpxio.h"
47
48 #define FACTOR 1000.0 /* Convert nm^2/ps to 10e-5 cm^2/s */
49 /* NORMAL = total diffusion coefficient (default). X,Y,Z is diffusion
50 coefficient in X,Y,Z direction. LATERAL is diffusion coefficient in
51 plane perpendicular to axis
52 */
53 enum { NOT_USED, NORMAL, X, Y, Z, LATERAL };
54
55 typedef struct {
56 real t0,delta_t,dim_factor;
57 real **data,*time,*data_x,*data_y,*data_z,*data_xy,*mass;
58 rvec **x0;
59 t_block *mols;
60 t_lsq **lsq;
61 int type,axis,natoms,nrestart,nnx,nframes,nlast,ngrp;
62 int *n_offs;
63 int **ndata;
64 } t_corr;
65
66 typedef real t_calc_func(t_corr *,int,atom_id[],int,rvec[]);
67 typedef void t_prep_data_func(t_corr *this,int gnx,atom_id index[],
68 rvec xcur[],rvec xprev[],matrix box);
69
70 static real thistime(t_corr *this)
71 {
72 return this->time[this->nframes];
73 }
74
75 static bool in_data(t_corr *this,int nx00)
76 {
77 return this->nframes-this->n_offs[nx00];
78 }
79
80 t_corr *init_corr(int nrgrp,int type,int axis,real dim_factor,
81 bool bMass,bool bMol,t_topology *top)
82 {
83 t_corr *this;
84 t_atoms *atoms;
85 int i;
86
87 snew(this,1);
88 this->type = type;
89 this->ngrp = nrgrp;
90 this->nframes = 0;
91 this->nlast = 0;
92
93 snew(this->ndata,nrgrp);
94 snew(this->data,nrgrp);
95 for(i=0; (i<nrgrp); i++) {
96 this->ndata[i] = NULL;
97 this->data[i] = NULL;
98 }
99 this->time = NULL;
100 this->lsq = NULL;
101 if (bMass) {
102 atoms=&top->atoms;
103 snew(this->mass,atoms->nr);
104 for(i=0; (i<atoms->nr); i++) {
105 this->mass[i]=atoms->atom[i].m;
106 }
107 }
108 else if (bMol) {
109 this->mols = &(top->blocks[ebMOLS]);
110 }
111
112 return this;
113 }
114
115 static void done_corr(t_corr *this)
116 {
117 int i;
118
119 sfree(this->n_offs);
120 for(i=0; (i<this->nrestart); i++)
121 sfree(this->x0[i]);
122 sfree(this->x0);
123 }
124
125 static void init_restart(t_corr *this,int nrestart,real dt)
126 {
127 int i;
128
129 this->nrestart = max(1,nrestart);
130 this->delta_t = (nrestart > 1) ? dt : 0;
131
132 printf("\nThe number of starting points you requested takes %lu"
133 " bytes memory\n\n",
134 (unsigned long) this->natoms*this->nrestart*sizeof(rvec));
135
136 snew(this->x0,this->nrestart);
137 for(i=0; (i<this->nrestart); i++)
138 snew(this->x0[i],this->natoms);
139 snew(this->n_offs,this->nrestart);
140 }
141
142 void lsq_y_ax_b2(int n, real x[], real y[], real *a, real *b,
143 real *sa,real *sb)
144 {
145 int i;
146 double yx,xx,sx,sy,sa2;
147
148 yx=xx=sx=sy=0.0;
149 for (i=0; (i < n); i++) {
150 yx+=y[i]*x[i];
151 xx+=x[i]*x[i];
152 sx+=x[i];
153 sy+=y[i];
154 }
155 *a=(n*yx-sy*sx)/(n*xx-sx*sx);
156 *b=(sy-(*a)*sx)/n;
157 sa2=0.0;
158 for(i=0; (i<n); i++) {
159 real tmp=(y[i]-(*a)*x[i]-(*b));
160 sa2+=tmp*tmp;
161 }
162 *sa=sqrt((sa2/(n*(n-2)))*(1.0/(xx/n-(sx/n)*(sx/n))));
163 *sb=(*sa)*sqrt(xx/n);
164 }
165
166 static void subtitle(FILE *out,t_corr *this)
167 {
168 int i;
169 real aa,bb,da,db,D;
170
171 for(i=0; (i<this->ngrp); i++) {
172 lsq_y_ax_b2(this->nframes,this->time,this->data[i],&aa,&bb,&da,&db);
173 D=aa*FACTOR/this->dim_factor;
174 fprintf(out,"# D[%d] = %.4f (1e-5 cm^2 s^-1)\n",i,D);
175 }
176 }
177
178 static void corr_print(t_corr *this,char *fn,char *title,char *yaxis,
179 bool bXvgr,bool bDiff)
180 {
181 FILE *out;
182 int i,j,imin;
183 real t,aa,bb,da,db;
184
185 /* real err_fac; */
186
187 /* err_fac=sqrt(2.0/3.0*sqrt(M_PI)); */
188 for(j=0; (j<this->ngrp); j++)
189 for(i=0; (i<this->nframes); i++)
190 this->data[j][i]/=this->ndata[j][i];
191 out=xvgropen(fn,title,"Time (ps)",yaxis);
192 if (bXvgr)
193 subtitle(out,this);
194
195 lsq_y_ax_b2(this->nframes,this->time,this->data[0],&aa,&bb,&da,&db);
196 if (bDiff)
197 imin = 1;
198 else
199 imin = 0;
200 for(i=imin; (i<this->nframes); i++) {
201 t = this->time[i]-this->time[0];
202 fprintf(out,"%10g",t);
203 for(j=0; (j<this->ngrp); j++)
204 fprintf(out," %10g",
205 bDiff ? 1000*this->data[j][i]/(6*t) : this->data[j][i]);
206 fprintf(out,"\n");
207 }
208 fclose(out);
209 }
210
211 /* called from corr_loop, to do the main calculations */
212 static void calc_corr(t_corr *this,int nr,int nx,atom_id index[],rvec xc[],
213 t_calc_func *calc1)
214 {
215 int nx0;
216 real g;
217
218 /* Check for new starting point */
219 if (this->nlast < this->nrestart) {
220 if ((thistime(this) >= (this->t0+this->nlast*this->delta_t)) && (nr==0)) {
221 printf("New starting point\n");
222 memcpy(this->x0[this->nlast],xc,this->natoms*sizeof(xc[0]));
223 this->n_offs[this->nlast]=this->nframes;
224 this->nlast++;
225 }
226 }
227
228 /* nx0 appears to be the number of new starting points,
229 * so for all starting points, call calc1.
230 */
231 for(nx0=0; (nx0<this->nlast); nx0++) {
232 g = calc1(this,nx,index,nx0,xc);
233 #ifdef DEBUG2
234 printf("g[%d]=%g\n",nx0,g);
235 #endif
236 this->data [nr][in_data(this,nx0)]+=g;
237 this->ndata[nr][in_data(this,nx0)]++;
238 }
239 }
240
241 static real calc1_norm(t_corr *this,int nx,atom_id index[],int nx0,rvec xc[])
242 {
243 int i,ix,m;
244 real g,r,r2;
245
246 g=0.0;
247
248 for(i=0; (i<nx); i++) {
249 ix=index[i];
250 r2=0.0;
251 switch (this->type) {
252 case NORMAL:
253 for(m=0; (m<DIM); m++) {
254 r = this->x0[nx0][ix][m]-xc[ix][m];
255 r2 += r*r;
256 }
257 break;
258 case X:
259 case Y:
260 case Z:
261 r = this->x0[nx0][ix][this->type-1]-xc[ix][this->type-1];
262 break;
263 case LATERAL:
264 for(m=0; (m<DIM); m++) {
265 if (m != this->axis) {
266 r = this->x0[nx0][ix][m]-xc[ix][m];
267 r2 += r*r;
268 }
269 }
270 break;
271 default:
272 fatal_error(0,"Error: did not expect option value %d",this->type);
273 }
274 g+=r2;
275 }
276 g/=nx;
277
278 return g;
279 }
280
281 /* this is the mean loop for the correlation type functions
282 * fx and nx are file pointers to things like read_first_x and
283 * read_next_x
284 */
285 void corr_loop(t_corr *this,char *fn,int gnx[],atom_id *index[],
286 t_calc_func *calc1,t_prep_data_func *prep1,
287 int nrestart,real dt,
288 t_first_x *fx,t_next_x *nx)
289 {
290 rvec *x[2];
291 real t;
292 int i,j,status,cur=0,maxframes=0;
293 #define prev (1-cur)
294 matrix box;
295
296 this->natoms=fx(&status,fn,&this->t0,&(x[prev]),box);
297 #ifdef DEBUG
298 fprintf(stderr,"Read %d atoms for first frame\n",this->natoms);
299 #endif
300
301 snew(x[cur],this->natoms);
302
303 init_restart(this,nrestart,dt);
304 memcpy(x[cur],x[prev],this->natoms*sizeof(x[prev][0]));
305 t=this->t0;
306 do {
307 if (this->nframes >= maxframes-1) {
308 if (maxframes==0) {
309 for(i=0; (i<this->ngrp); i++) {
310 this->ndata[i] = NULL;
311 this->data[i] = NULL;
312 }
313 this->time = NULL;
314 }
315 maxframes+=10;
316 for(i=0; (i<this->ngrp); i++) {
317 srenew(this->ndata[i],maxframes);
318 srenew(this->data[i],maxframes);
319 for(j=maxframes-10; j<maxframes; j++) {
320 this->ndata[i][j]=0;
321 this->data[i][j]=0;
322 }
323 }
324 srenew(this->time,maxframes);
325 }
326
327 this->time[this->nframes]=t;
328
329 /* loop over all groups in index file */
330 for(i=0; (i<this->ngrp); i++) {
331 /* nice for putting things in boxes and such */
332 prep1(this,gnx[i],index[i],x[cur],x[prev],box);
333 /* calculate something useful, like mean square displacements */
334 calc_corr(this,i,gnx[i],index[i],x[cur],calc1);
335 }
336 cur=prev;
337
338 this->nframes++;
339 } while (nx(status,&t,this->natoms,x[cur],box));
340 fprintf(stderr,"\n");
341
342 close_trj(status);
343 }
344
345 static void prep_data_mol(t_corr *this,int gnx,atom_id index[],
346 rvec xcur[],rvec xprev[],matrix box)
347 {
348 int i,j,k,m,ind;
349 rvec hbox;
350
351 /* Remove periodicity */
352 for(m=0; (m<DIM); m++)
353 hbox[m]=0.5*box[m][m];
354 for(i=0; (i<gnx); i++) {
355 ind=index[i];
356 for(j=this->mols->index[ind]; (j<this->mols->index[ind+1]); j++) {
357 k=this->mols->a[j];
358 for(m=0; (m<DIM); m++) {
359 while(xcur[k][m]-xprev[k][m] <= hbox[m])
360 xcur[k][m] += box[m][m];
361 while(xcur[k][m]-xprev[k][m] > hbox[m])
362 xcur[k][m] -= box[m][m];
363 }
364 }
365 }
366 }
367
368 static real calc_one_mw(t_corr *this,int ix,int nx0,rvec xc[],real *tm)
369 {
370 real r2,r,mm;
371 int m;
372
373 mm=this->mass[ix];
374 if (mm < 1)
375 return 0;
376 (*tm)+=this->mass[ix];
377 r2=0.0;
378 switch (this->type) {
379 case NORMAL:
380 for(m=0; (m<DIM); m++) {
381 r = this->x0[nx0][ix][m]-xc[ix][m];
382 r2 += this->mass[ix]*r*r;
383 }
384 break;
385 case X:
386 case Y:
387 case Z:
388 r = this->x0[nx0][ix][this->type-1]-xc[ix][this->type-1];
389 r2 = this->mass[ix]*r*r;
390 break;
391 case LATERAL:
392 for(m=0; (m<DIM); m++) {
393 if (m != this->axis) {
394 r = this->x0[nx0][ix][m]-xc[ix][m];
395 r2 += this->mass[ix]*r*r;
396 }
397 }
398 break;
399 default:
400 fatal_error(0,"Options got screwed. Did not expect value %d\n",this->type);
401 } /* end switch */
402 return r2;
403 }
404
405 static real calc1_mw(t_corr *this,int nx,atom_id index[],int nx0,rvec xc[])
406 {
407 int i;
408 real g,tm;
409
410 g=tm=0.0;
411 for(i=0; (i<nx); i++)
412 g+=calc_one_mw(this,index[i],nx0,xc,&tm);
413
414 g/=tm;
415
416 return g;
417 }
418
419 static void prep_data_norm(t_corr *this,int gnx,atom_id index[],
420 rvec xcur[],rvec xprev[],matrix box)
421 {
422 int i,m,ind;
423 rvec hbox;
424
425 /* Remove periodicity */
426 for(m=0; (m<DIM); m++)
427 hbox[m]=0.5*box[m][m];
428 for(i=0; (i<gnx); i++) {
429 ind=index[i];
430 for(m=0; (m<DIM); m++) {
431 while(xcur[ind][m]-xprev[ind][m] <= hbox[m])
432 xcur[ind][m] += box[m][m];
433 while(xcur[ind][m]-xprev[ind][m] > hbox[m])
434 xcur[ind][m] -= box[m][m];
435 }
436 }
437 }
438
439 static real calc1_mol(t_corr *this,int nx,atom_id index[],int nx0,rvec xc[])
440 {
441 int i,ii,j;
442 real g,mm,gtot,tt;
443
444 if (this->lsq == NULL) {
445 this->nnx = nx;
446 snew(this->lsq,this->nrestart);
447 for(i=0; (i<this->nrestart); i++)
448 snew(this->lsq[i],nx);
449 }
450 tt=this->time[in_data(this,nx0)];
451 gtot=0;
452 for(i=0; (i<nx); i++) {
453 ii=index[i];
454 g=mm=0;
455 for(j=this->mols->index[ii]; (j<this->mols->index[ii+1]); j++)
456 g += calc_one_mw(this,this->mols->a[j],nx0,xc,&mm);
457
458 g/=mm;
459 gtot+=g;
460 add_lsq(&(this->lsq[nx0][i]),tt,g);
461 }
462 return gtot/nx;
463 }
464
465 void printdist(t_corr *this,char *fn,char *difn)
466 {
467 #define NDIST 100
468 FILE *out;
469 int *ndist;
470 real *diff,mind=0,maxd=0,dav,a,b;
471 t_lsq lsq1;
472 int i,j;
473
474 out=xvgropen(difn,"Diffusion Coefficients / Molecule","Molecule","D");
475
476 dav=0;
477 snew(diff,this->nnx);
478 for(i=0; (i<this->nnx); i++) {
479 init_lsq(&lsq1);
480 for(j=0; (j<this->nrestart); j++) {
481 lsq1.sx+=this->lsq[j][i].sx;
482 lsq1.sy+=this->lsq[j][i].sy;
483 lsq1.xx+=this->lsq[j][i].xx;
484 lsq1.yx+=this->lsq[j][i].yx;
485 lsq1.np+=this->lsq[j][i].np;
486 }
487 get_lsq_ab(&lsq1,&a,&b);
488 diff[i]=a*FACTOR/this->dim_factor;
489 fprintf(out,"%10d %10g\n",i,diff[i]);
490 if (dav == 0) {
491 mind=maxd=dav=diff[i];
492 }
493 else {
494 mind=min(mind,diff[i]);
495 maxd=max(maxd,diff[i]);
496 dav+=diff[i];
497 }
498 }
499 dav/=this->nnx;
500 fclose(out);
501 xvgr_file(difn,"-graphtype bar");
502
503 ndist=(int *)calloc(NDIST+1,sizeof(*ndist));
504 for(i=0; i<NDIST+1; i++)
505 ndist[i]=0;
506 for(i=0; (i<this->nnx); i++) {
507 int index=(int)(0.5+NDIST*(diff[i]-mind)/(maxd-mind));
508 if ((index >= 0) && (index <= NDIST))
509 ndist[index]++;
510 }
511 out=xvgropen(fn,"Distribution of Diffusion Constants","D ()","Number");
512 for(i=0; (i<=NDIST); i++)
513 fprintf(out,"%10g %10d\n",
514 mind+(i*(maxd-mind))/NDIST,ndist[i]);
515 fclose(out);
516 xvgr_file(fn,NULL);
517
518 }
519
520 void do_corr(int NFILE, t_filenm fnm[],int nrgrp,
521 t_topology *top,bool bMol,bool bMW,
522 int type,real dim_factor,int axis,
523 int nrestart,real dt)
524 {
525 t_corr *msd;
526 t_first_x *fx;
527 t_next_x *nx;
528 int *gnx;
529 atom_id **index;
530 char **grpname;
531
532 fx=read_first_x;
533 nx=read_next_x;
534
535 gnx = (int *)calloc(nrgrp,sizeof(int));
536 index = (atom_id **)calloc(nrgrp,sizeof(atom_id *));
537 grpname = (char **)calloc(nrgrp,sizeof(char *));
538
539 get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),nrgrp,gnx,index,grpname);
540
541 msd = init_corr(nrgrp,type,axis,dim_factor,bMW,bMol,top);
542
543 corr_loop(msd,ftp2fn(efTRX,NFILE,fnm),gnx,index,
544 bMW ? calc1_mw : (bMol ? calc1_mol : calc1_norm),
545 bMol ? prep_data_mol : prep_data_norm,nrestart,dt,
546 fx,nx);
547
548 if (opt2bSet("-d",NFILE,fnm))
549 corr_print(msd,opt2fn("-d",NFILE,fnm),"Diffusion constant",
550 "D (10\\S5\\Ncm\\S2\\Ns\\S-1\\N)",TRUE,TRUE);
551 corr_print(msd,opt2fn("-o",NFILE,fnm),
552 "Mean Square Displacement",
553 "MSD (nm\\S2\\N)",TRUE,FALSE);
554
555 xvgr_file(ftp2fn(efXVG,NFILE,fnm),"-nxy");
556 if (bMol)
557 printdist(msd,opt2fn("-m",NFILE,fnm),opt2fn("-d",NFILE,fnm));
558 }
559
560 int main(int argc,char *argv[])
561 {
562 static char *desc[] = {
563 "g_msd computes the mean square displacement (MSD) of atoms from",
564 "their initial positions. This provides an easy way to compute",
565 "the diffusion constant using the Einstein relation.[PAR]",
566 "If the -d option is given, the diffusion constant will be printed in",
567 "addition to the MSD[PAR]",
568 "Mean Square Displacement calculations and Correlation functions",
569 "can be calculated more accurately, when using multiple starting",
570 "points (see also Gromacs Manual). You can select the number of",
571 "starting points, and the interval (in picoseconds) between starting",
572 "points. More starting points implies more CPU time."
573 };
574 static char *normtype[]= { NULL,"no","x","y","z",NULL };
575 static char *axtitle[] = { NULL,"no","x","y","z",NULL };
576 static int ngroup = 1;
577 static int nrestart = 1;
578 static real dt = 0;
579 static bool bMW = TRUE;
580 t_pargs pa[] = {
581 { "-type", FALSE, etENUM, {normtype},
582 "Compute diffusion coefficient in one direction" },
583 { "-lateral", FALSE, etENUM, {axtitle},
584 "Calculate the lateral diffusion in a plane perpendicular to" },
585 { "-ngroup", FALSE, etINT, {&ngroup},
586 "Number of groups to calculate MSD for" },
587 { "-mw", FALSE, etBOOL, {&bMW},
588 "Mass weighted MSD" },
589 { "-nrestart",FALSE, etINT, {&nrestart},
590 "Number of restarting points in trajectory" },
591 { "-dt", FALSE, etREAL, {&dt},
592 "Time between restarting points in trajectory (only with -nrestart > 1)" }
593 };
594 static char *bugs[] = {
595 "The diffusion constant given in the title of the graph for lateral"
596 " diffusion has to be multiplied by 6/4"
597 };
598
599 t_filenm fnm[] = {
600 { efTRX, NULL, NULL, ffREAD },
601 { efTPS, NULL, NULL, ffREAD },
602 { efNDX, NULL, NULL, ffOPTRD },
603 { efXVG, NULL, "msd", ffWRITE },
604 { efXVG, "-m", "mol", ffOPTWR },
605 { efXVG, "-d", "diff",ffOPTWR }
606 };
607 #define NFILE asize(fnm)
608
609 t_topology top;
610 matrix box;
611 char title[256];
612 rvec *xdum;
613 bool bMol,bTop;
614 int axis,type;
615 real dim_factor;
616
617 CopyRight(stdout,argv[0]);
618
619 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
620 NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
621
622 if (ngroup < 1)
623 fatal_error(0,"Must have at least 1 group (now %d)",ngroup);
624
625 bMol = opt2bSet("-m",NFILE,fnm);
626
627 if (normtype[0][0]!='n') {
628 type = normtype[0][0] - 'x'+1; /* See defines above */
629 dim_factor = 2.0;
630 }
631 else {
632 type = NORMAL;
633 dim_factor = 6.0;
634 }
635 if ((type==NORMAL) && (axtitle[0][0]!='n')) {
636 type=LATERAL;
637 dim_factor = 4.0;
638 axis = (axtitle[0][0] - 'x'); /* See defines above */
639 }
640 else
641 axis = 0;
642 fprintf(stderr,"type = %d, axis = %d, dim_factor = %g\n",
643 type,axis,dim_factor);
644 bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&xdum,NULL,box,bMW);
645 if (bMol && !bTop)
646 fatal_error(0,"Could not read a topology form %s. Try a tpr file instead.",
647 ftp2fn(efTPS,NFILE,fnm));
648
649 do_corr(NFILE,fnm,ngroup,&top,bMol,bMW,type,dim_factor,axis,nrestart,dt);
650
651 thanx(stdout);
652
653 return 0;
654 }
got merge a long time ago