4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
55 #define FACTOR 1000.0 /* Convert nm^2/ps to 10e-5 cm^2/s */
56 /* NORMAL = total diffusion coefficient (default). X,Y,Z is diffusion
57 coefficient in X,Y,Z direction. LATERAL is diffusion coefficient in
58 plane perpendicular to axis
60 enum { NOT_USED
, NORMAL
, X
, Y
, Z
, LATERAL
};
63 real t0
,delta_t
,dim_factor
;
64 real
**data
,*time
,*data_x
,*data_y
,*data_z
,*data_xy
,*mass
;
68 int type
,axis
,natoms
,nrestart
,nnx
,nframes
,nlast
,ngrp
;
73 typedef real
t_calc_func(t_corr
*,int,atom_id
[],int,rvec
[]);
74 typedef void t_prep_data_func(t_corr
*this,int gnx
,atom_id index
[],
75 rvec xcur
[],rvec xprev
[],matrix box
);
77 static real
thistime(t_corr
*this)
79 return this->time
[this->nframes
];
82 static bool in_data(t_corr
*this,int nx00
)
84 return this->nframes
-this->n_offs
[nx00
];
87 t_corr
*init_corr(int nrgrp
,int type
,int axis
,real dim_factor
,
88 bool bMass
,bool bMol
,real dt
,t_topology
*top
)
104 this->dim_factor
= dim_factor
;
106 snew(this->ndata
,nrgrp
);
107 snew(this->data
,nrgrp
);
108 for(i
=0; (i
<nrgrp
); i
++) {
109 this->ndata
[i
] = NULL
;
110 this->data
[i
] = NULL
;
116 snew(this->mass
,atoms
->nr
);
117 for(i
=0; (i
<atoms
->nr
); i
++) {
118 this->mass
[i
]=atoms
->atom
[i
].m
;
122 this->mols
= &(top
->blocks
[ebMOLS
]);
127 static void done_corr(t_corr
*this)
132 for(i
=0; (i
<this->nrestart
); i
++)
137 static void init_restart(t_corr
*this)
143 srenew(this->x0
,this->nrestart
);
144 snew(this->x0
[this->nrestart
-1],this->natoms
);
145 srenew(this->n_offs
,this->nrestart
);
148 static void corr_print(t_corr
*this,char *fn
,char *title
,char *yaxis
,
149 real beginfit
,real endfit
,
150 real
*DD
,real
*SigmaD
,char *grpname
[])
155 out
=xvgropen(fn
,title
,xvgr_tlabel(),yaxis
);
157 fprintf(out
,"# Diffusion constants fitted from time %g to %g (%s)\n",
158 beginfit
,endfit
,time_unit());
159 for(i
=0; i
<this->ngrp
; i
++)
160 fprintf(out
,"# D[%10s] = %.3f (+/- %.3f) (1e-5 cm^2/s)\n",
161 grpname
[i
],DD
[i
],SigmaD
[i
]);
163 for(i
=0; i
<this->nframes
; i
++) {
164 fprintf(out
,"%10g",convert_time(this->time
[i
]));
165 for(j
=0; j
<this->ngrp
; j
++)
166 fprintf(out
," %10g",this->data
[j
][i
]);
172 /* called from corr_loop, to do the main calculations */
173 static void calc_corr(t_corr
*this,int nr
,int nx
,atom_id index
[],rvec xc
[],
179 /* Check for new starting point */
180 if (this->nlast
< this->nrestart
) {
181 if ((thistime(this) >= (this->nlast
*this->delta_t
)) && (nr
==0)) {
182 memcpy(this->x0
[this->nlast
],xc
,this->natoms
*sizeof(xc
[0]));
183 this->n_offs
[this->nlast
]=this->nframes
;
188 /* nx0 appears to be the number of new starting points,
189 * so for all starting points, call calc1.
191 for(nx0
=0; (nx0
<this->nlast
); nx0
++) {
192 g
= calc1(this,nx
,index
,nx0
,xc
);
194 printf("g[%d]=%g\n",nx0
,g
);
196 this->data
[nr
][in_data(this,nx0
)]+=g
;
197 this->ndata
[nr
][in_data(this,nx0
)]++;
201 static real
calc1_norm(t_corr
*this,int nx
,atom_id index
[],int nx0
,rvec xc
[])
208 for(i
=0; (i
<nx
); i
++) {
211 switch (this->type
) {
213 for(m
=0; (m
<DIM
); m
++) {
214 r
= this->x0
[nx0
][ix
][m
]-xc
[ix
][m
];
221 r
= this->x0
[nx0
][ix
][this->type
-X
]-xc
[ix
][this->type
-X
];
224 for(m
=0; (m
<DIM
); m
++) {
225 if (m
!= this->axis
) {
226 r
= this->x0
[nx0
][ix
][m
]-xc
[ix
][m
];
232 fatal_error(0,"Error: did not expect option value %d",this->type
);
241 /* this is the mean loop for the correlation type functions
242 * fx and nx are file pointers to things like read_first_x and
245 void corr_loop(t_corr
*this,char *fn
,int gnx
[],atom_id
*index
[],
246 t_calc_func
*calc1
,t_prep_data_func
*prep1
,real dt
)
250 int i
,j
,status
,cur
=0,maxframes
=0;
254 this->natoms
=read_first_x(&status
,fn
,&this->t0
,&(x
[prev
]),box
);
256 fprintf(stderr
,"Read %d atoms for first frame\n",this->natoms
);
259 snew(x
[cur
],this->natoms
);
261 /* init_restart(this,nrestart,dt); */
262 memcpy(x
[cur
],x
[prev
],this->natoms
*sizeof(x
[prev
][0]));
265 if (bRmod(t
,this->t0
,dt
))
267 if (this->nframes
>= maxframes
-1) {
269 for(i
=0; (i
<this->ngrp
); i
++) {
270 this->ndata
[i
] = NULL
;
271 this->data
[i
] = NULL
;
276 for(i
=0; (i
<this->ngrp
); i
++) {
277 srenew(this->ndata
[i
],maxframes
);
278 srenew(this->data
[i
],maxframes
);
279 for(j
=maxframes
-10; j
<maxframes
; j
++) {
284 srenew(this->time
,maxframes
);
287 this->time
[this->nframes
] = t
- this->t0
;
289 /* loop over all groups in index file */
290 for(i
=0; (i
<this->ngrp
); i
++) {
291 /* nice for putting things in boxes and such */
292 prep1(this,gnx
[i
],index
[i
],x
[cur
],x
[prev
],box
);
293 /* calculate something useful, like mean square displacements */
294 calc_corr(this,i
,gnx
[i
],index
[i
],x
[cur
],calc1
);
299 } while (read_next_x(status
,&t
,this->natoms
,x
[cur
],box
));
300 fprintf(stderr
,"\nUsed %d restart points spaced %g %s over %g %s\n\n",
302 convert_time(dt
), time_unit(),
303 convert_time(this->time
[this->nframes
-1]), time_unit() );
308 static void prep_data_mol(t_corr
*this,int gnx
,atom_id index
[],
309 rvec xcur
[],rvec xprev
[],matrix box
)
314 /* Remove periodicity */
315 for(m
=0; (m
<DIM
); m
++)
316 hbox
[m
]=0.5*box
[m
][m
];
317 for(i
=0; i
<gnx
; i
++) {
319 for(j
=this->mols
->index
[ind
]; j
<this->mols
->index
[ind
+1]; j
++) {
321 for(m
=DIM
-1; m
>=0; m
--) {
322 while (xcur
[k
][m
]-xprev
[k
][m
] <= -hbox
[m
])
324 xcur
[k
][d
] += box
[m
][d
];
325 while (xcur
[k
][m
]-xprev
[k
][m
] > hbox
[m
])
327 xcur
[k
][d
] -= box
[m
][d
];
333 static real
calc_one_mw(t_corr
*this,int ix
,int nx0
,rvec xc
[],real
*tm
)
341 (*tm
)+=this->mass
[ix
];
343 switch (this->type
) {
345 for(m
=0; (m
<DIM
); m
++) {
346 r
= this->x0
[nx0
][ix
][m
]-xc
[ix
][m
];
347 r2
+= this->mass
[ix
]*r
*r
;
353 r
= this->x0
[nx0
][ix
][this->type
-X
]-xc
[ix
][this->type
-X
];
354 r2
= this->mass
[ix
]*r
*r
;
357 for(m
=0; (m
<DIM
); m
++) {
358 if (m
!= this->axis
) {
359 r
= this->x0
[nx0
][ix
][m
]-xc
[ix
][m
];
360 r2
+= this->mass
[ix
]*r
*r
;
365 fatal_error(0,"Options got screwed. Did not expect value %d\n",this->type
);
370 static real
calc1_mw(t_corr
*this,int nx
,atom_id index
[],int nx0
,rvec xc
[])
376 for(i
=0; (i
<nx
); i
++)
377 g
+=calc_one_mw(this,index
[i
],nx0
,xc
,&tm
);
384 static void prep_data_norm(t_corr
*this,int gnx
,atom_id index
[],
385 rvec xcur
[],rvec xprev
[],matrix box
)
390 /* Remove periodicity */
391 for(m
=0; (m
<DIM
); m
++)
392 hbox
[m
]=0.5*box
[m
][m
];
393 for(i
=0; (i
<gnx
); i
++) {
395 for(m
=DIM
-1; m
>=0; m
--) {
396 while(xcur
[ind
][m
]-xprev
[ind
][m
] <= hbox
[m
])
397 rvec_inc(xcur
[ind
],box
[m
]);
398 while(xcur
[ind
][m
]-xprev
[ind
][m
] > hbox
[m
])
399 rvec_dec(xcur
[ind
],box
[m
]);
404 static real
calc1_mol(t_corr
*this,int nx
,atom_id index
[],int nx0
,rvec xc
[])
409 if (this->lsq
== NULL
) {
411 snew(this->lsq
,this->nrestart
);
412 for(i
=0; (i
<this->nrestart
); i
++)
413 snew(this->lsq
[i
],nx
);
415 tt
=this->time
[in_data(this,nx0
)];
417 for(i
=0; (i
<nx
); i
++) {
420 for(j
=this->mols
->index
[ii
]; (j
<this->mols
->index
[ii
+1]); j
++)
421 g
+= calc_one_mw(this,this->mols
->a
[j
],nx0
,xc
,&mm
);
425 add_lsq(&(this->lsq
[nx0
][i
]),tt
,g
);
430 void printmol(t_corr
*this,char *fn
)
436 real a
,b
,*D
,Dav
,D2av
,VarD
;
438 out
=xvgropen(fn
,"Diffusion Coefficients / Molecule","Molecule","D");
442 for(i
=0; (i
<this->nnx
); i
++) {
444 for(j
=0; (j
<this->nrestart
); j
++) {
445 lsq1
.sx
+=this->lsq
[j
][i
].sx
;
446 lsq1
.sy
+=this->lsq
[j
][i
].sy
;
447 lsq1
.xx
+=this->lsq
[j
][i
].xx
;
448 lsq1
.yx
+=this->lsq
[j
][i
].yx
;
449 lsq1
.np
+=this->lsq
[j
][i
].np
;
451 get_lsq_ab(&lsq1
,&a
,&b
);
452 D
[i
] = a
*FACTOR
/this->dim_factor
;
455 fprintf(out
,"%10d %10g\n",i
,D
[i
]);
458 do_view(fn
,"-graphtype bar");
460 /* Compute variance, stddev and error */
463 VarD
= D2av
- sqr(Dav
);
464 printf("<D> = %.3f Std. Dev. = %.3f Error = %.3f\n",
465 Dav
,sqrt(VarD
),sqrt(VarD
/this->nnx
));
470 void do_corr(char *trx_file
, char *ndx_file
, char *msd_file
, char *mol_file
,
471 int nrgrp
, t_topology
*top
,bool bMW
,
472 int type
,real dim_factor
,int axis
,
473 real dt
,real beginfit
,real endfit
)
480 real
*DD
,*SigmaD
,a
,a2
,b
;
486 if (mol_file
&& !ndx_file
) {
487 gnx
[0] = top
->blocks
[ebMOLS
].nr
;
488 snew(index
[0],gnx
[0]);
489 grpname
[0] = "Molecules";
490 for(i
=0; (i
<gnx
[0]); i
++)
494 get_index(&top
->atoms
,ndx_file
,nrgrp
,gnx
,index
,grpname
);
496 msd
= init_corr(nrgrp
,type
,axis
,dim_factor
,bMW
,(mol_file
!=NULL
),dt
,top
);
498 corr_loop(msd
,trx_file
,gnx
,index
,
499 (mol_file
!=NULL
) ? calc1_mol
: (bMW
? calc1_mw
: calc1_norm
),
500 (mol_file
!=NULL
) ? prep_data_mol
: prep_data_norm
,dt
);
502 /* Correct for the number of points */
503 for(j
=0; (j
<msd
->ngrp
); j
++)
504 for(i
=0; (i
<msd
->nframes
); i
++)
505 msd
->data
[j
][i
] /= msd
->ndata
[j
][i
];
508 printmol(msd
,mol_file
);
512 for(i0
=0; i0
<msd
->nframes
&& msd
->time
[i0
]<beginfit
; i0
++)
516 endfit
= msd
->time
[i1
-1];
518 for(i1
=i0
; i1
<msd
->nframes
&& msd
->time
[i1
]<=endfit
; i1
++)
522 fprintf(stdout
,"Not enough points for fitting (%d).\n"
523 "Can not determine the diffusion constant.\n",N
);
526 snew(SigmaD
,msd
->ngrp
);
527 for(j
=0; j
<msd
->ngrp
; j
++) {
529 lsq_y_ax_b(N
/2,&(msd
->time
[i0
]),&(msd
->data
[j
][i0
]),&a
,&b
);
530 lsq_y_ax_b(N
/2,&(msd
->time
[i0
+N
/2]),&(msd
->data
[j
][i0
+N
/2]),&a2
,&b
);
531 SigmaD
[j
] = fabs(a
-a2
);
534 lsq_y_ax_b(N
,&(msd
->time
[i0
]),&(msd
->data
[j
][i0
]),&(DD
[j
]),&b
);
535 DD
[j
] *= FACTOR
/msd
->dim_factor
;
536 SigmaD
[j
] *= FACTOR
/msd
->dim_factor
;
537 fprintf(stdout
,"D[%10s] %.3f (+/- %.3f) 1e-5 cm^2/s\n",
538 grpname
[j
],DD
[j
],SigmaD
[j
]);
541 /* Print mean square displacement */
542 corr_print(msd
,msd_file
,
543 "Mean Square Displacement",
545 beginfit
,endfit
,DD
,SigmaD
,grpname
);
548 int gmx_msd(int argc
,char *argv
[])
550 static char *desc
[] = {
551 "g_msd computes the mean square displacement (MSD) of atoms from",
552 "their initial positions. This provides an easy way to compute",
553 "the diffusion constant using the Einstein relation.",
554 "The diffusion constant is calculated by least squares fitting a",
555 "straight line through the MSD from [TT]-beginfit[tt] to",
556 "[TT]-endfit[tt]. An error estimate given, which is the difference",
557 "of the diffusion coefficients obtained from fits over the two halfs",
558 "of the fit interval.[PAR]",
559 "Option [TT]-mol[tt] plots the MSD for molecules, this implies",
560 "[TT]-mw[tt], i.e. for each inidividual molecule an diffusion constant",
561 "is computed. When using an index file, it should contain molecule",
562 "numbers instead of atom numbers.",
563 "Using this option one also gets an accurate error estimate",
564 "based on the statistics between individual molecules. Since one usually",
565 "is interested in self-diffusion at infinite dilution this is probably",
566 "the most useful number.[PAR]",
568 static char *normtype
[]= { NULL
,"no","x","y","z",NULL
};
569 static char *axtitle
[] = { NULL
,"no","x","y","z",NULL
};
570 static int ngroup
= 1;
572 static real beginfit
= 0;
573 static real endfit
= -1;
574 static bool bMW
= TRUE
;
576 { "-type", FALSE
, etENUM
, {normtype
},
577 "Compute diffusion coefficient in one direction" },
578 { "-lateral", FALSE
, etENUM
, {axtitle
},
579 "Calculate the lateral diffusion in a plane perpendicular to" },
580 { "-ngroup", FALSE
, etINT
, {&ngroup
},
581 "Number of groups to calculate MSD for" },
582 { "-mw", FALSE
, etBOOL
, {&bMW
},
583 "Mass weighted MSD" },
584 { "-trestart",FALSE
, etTIME
, {&dt
},
585 "Time between restarting points in trajectory (%t)" },
586 { "-beginfit",FALSE
, etTIME
, {&beginfit
},
587 "Start time for fitting the MSD (%t)" },
588 { "-endfit",FALSE
, etTIME
, {&endfit
},
589 "End time for fitting the MSD (%t), -1 is till end" }
593 { efTRX
, NULL
, NULL
, ffREAD
},
594 { efTPS
, NULL
, NULL
, ffREAD
},
595 { efNDX
, NULL
, NULL
, ffOPTRD
},
596 { efXVG
, NULL
, "msd", ffWRITE
},
597 { efXVG
, "-mol", "diff_mol",ffOPTWR
}
599 #define NFILE asize(fnm)
604 char *trx_file
, *tps_file
, *ndx_file
, *msd_file
, *mol_file
;
610 CopyRight(stderr
,argv
[0]);
612 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_BE_NICE
,
613 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
);
614 trx_file
= ftp2fn_null(efTRX
,NFILE
,fnm
);
615 tps_file
= ftp2fn_null(efTPS
,NFILE
,fnm
);
616 ndx_file
= ftp2fn_null(efNDX
,NFILE
,fnm
);
617 msd_file
= ftp2fn_null(efXVG
,NFILE
,fnm
);
618 mol_file
= opt2fn_null("-mol",NFILE
,fnm
);
621 fatal_error(0,"Must have at least 1 group (now %d)",ngroup
);
625 fprintf(stderr
,"Calculating diffusion coefficients for molecules.\n");
628 if (normtype
[0][0]!='n') {
629 type
= normtype
[0][0] - 'x' + X
; /* See defines above */
636 if ((type
==NORMAL
) && (axtitle
[0][0]!='n')) {
639 axis
= (axtitle
[0][0] - 'x'); /* See defines above */
643 bTop
=read_tps_conf(tps_file
,title
,&top
,&xdum
,NULL
,box
,bMW
);
644 if (mol_file
&& !bTop
)
645 fatal_error(0,"Could not read a topology from %s. Try a tpr file instead.",
648 do_corr(trx_file
,ndx_file
,msd_file
,mol_file
,ngroup
,
649 &top
,bMW
,type
,dim_factor
,axis
,dt
,beginfit
,endfit
);
651 view_all(NFILE
, fnm
);