Bug Summary

File:gromacs/gmxana/gmx_current.c
Location:line 284, column 5
Description:Value stored to 'corint' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
8 *
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
13 *
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
23 *
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
31 *
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
34 */
35#ifdef HAVE_CONFIG_H1
36#include <config.h>
37#endif
38
39#include <stdlib.h>
40
41#include "gromacs/commandline/pargs.h"
42#include "typedefs.h"
43#include "gromacs/utility/smalloc.h"
44#include "gromacs/math/vec.h"
45#include "gromacs/fileio/tpxio.h"
46#include "gromacs/fileio/trxio.h"
47#include "gromacs/fileio/xvgr.h"
48#include "rmpbc.h"
49#include "pbc.h"
50#include "physics.h"
51#include "index.h"
52#include "gromacs/statistics/statistics.h"
53#include "gmx_ana.h"
54#include "macros.h"
55
56#include "gromacs/utility/fatalerror.h"
57
58#define SQR(x)(pow(x, 2.0)) (pow(x, 2.0))
59#define EPSI0((5.72765E-4)*(1.60217733e-19)*(1.60217733e-19)*(6.0221367e23
)/((1e3)*(1e-9))) (EPSILON0(5.72765E-4)*E_CHARGE(1.60217733e-19)*E_CHARGE(1.60217733e-19)*AVOGADRO(6.0221367e23)/(KILO(1e3)*NANO(1e-9))) /* EPSILON0 in SI units */
60
61static void index_atom2mol(int *n, int *index, t_block *mols)
62{
63 int nat, i, nmol, mol, j;
64
65 nat = *n;
66 i = 0;
67 nmol = 0;
68 mol = 0;
69 while (i < nat)
70 {
71 while (index[i] > mols->index[mol])
72 {
73 mol++;
74 if (mol >= mols->nr)
75 {
76 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 76, "Atom index out of range: %d", index[i]+1);
77 }
78 }
79 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
80 {
81 if (i >= nat || index[i] != j)
82 {
83 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 83, "The index group does not consist of whole molecules");
84 }
85 i++;
86 }
87 index[nmol++] = mol;
88 }
89
90 fprintf(stderrstderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
91
92 *n = nmol;
93}
94
95static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
96{
97
98 real mtot;
99 real qtot;
100 real qall;
101 int i, j, k, l;
102 int ai, ci;
103 gmx_bool bNEU;
104 ai = 0;
105 ci = 0;
106 qall = 0.0;
107
108
109
110 for (i = 0; i < top.mols.nr; i++)
111 {
112 k = top.mols.index[i];
113 l = top.mols.index[i+1];
114 mtot = 0.0;
115 qtot = 0.0;
116
117 for (j = k; j < l; j++)
118 {
119 mtot += top.atoms.atom[j].m;
120 qtot += top.atoms.atom[j].q;
121 }
122
123 for (j = k; j < l; j++)
124 {
125 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
126 mass2[j] = top.atoms.atom[j].m/mtot;
127 qmol[j] = qtot;
128 }
129
130
131 qall += qtot;
132
133 if (qtot < 0.0)
134 {
135 ai++;
136 }
137 if (qtot > 0.0)
138 {
139 ci++;
140 }
141 }
142
143 if (abs(qall) > 0.01)
144 {
145 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
146 bNEU = FALSE0;
147 }
148 else
149 {
150 bNEU = TRUE1;
151 }
152
153 return bNEU;
154
155}
156
157static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
158{
159
160 rvec hbox;
161 int d, i, m;
162
163 for (d = 0; d < DIM3; d++)
164 {
165 hbox[d] = 0.5*box[d][d];
166 }
167 for (i = 0; i < natoms; i++)
168 {
169 for (m = DIM3-1; m >= 0; m--)
170 {
171 while (x[i][m]-xp[i][m] <= -hbox[m])
172 {
173 for (d = 0; d <= m; d++)
174 {
175 x[i][d] += box[m][d];
176 }
177 }
178 while (x[i][m]-xp[i][m] > hbox[m])
179 {
180 for (d = 0; d <= m; d++)
181 {
182 x[i][d] -= box[m][d];
183 }
184 }
185 }
186 }
187}
188
189static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
190 rvec fr[], rvec mj, real mass2[], real qmol[])
191{
192
193 int i, j, k, l;
194 rvec tmp;
195 rvec center;
196 rvec mt1, mt2;
197 t_pbc pbc;
198
199
200 calc_box_center(ecenterRECT, box, center);
201
202 if (!bNoJump)
203 {
204 set_pbc(&pbc, ePBC, box);
205 }
206
207 clear_rvec(tmp);
208
209
210 for (i = 0; i < isize; i++)
211 {
212 clear_rvec(mt1);
213 clear_rvec(mt2);
214 k = top.mols.index[index0[i]];
215 l = top.mols.index[index0[i+1]];
216 for (j = k; j < l; j++)
217 {
218 svmul(mass2[j], fr[j], tmp);
219 rvec_inc(mt1, tmp);
220 }
221
222 if (bNoJump)
223 {
224 svmul(qmol[k], mt1, mt1);
225 }
226 else
227 {
228 pbc_dx(&pbc, mt1, center, mt2);
229 svmul(qmol[k], mt2, mt1);
230 }
231
232 rvec_inc(mj, mt1);
233
234 }
235
236}
237
238
239static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
240{
241
242 /* bCOR determines if the correlation is computed via
243 * static properties (FALSE) or the correlation integral (TRUE)
244 */
245
246 real epsilon = 0.0;
247
248
249 if (bCOR)
250 {
251 epsilon = md2-2.0*cor+mj2;
252 }
253 else
254 {
255 epsilon = md2+mj2+2.0*cor;
256 }
257
258 if (eps_rf == 0.0)
259 {
260 epsilon = 1.0+prefactor*epsilon;
261
262 }
263 else
264 {
265 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
266 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
267 }
268
269
270 return epsilon;
271
272}
273
274
275static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
276{
277
278 int i;
279 real corint;
280 real deltat = 0.0;
281 real rfr;
282 real sigma = 0.0;
283 real sigma_ret = 0.0;
284 corint = 0.0;
Value stored to 'corint' is never read
285
286 if (nfr > 1)
287 {
288 i = 0;
289
290 while (i < nfr)
291 {
292
293 rfr = (real) (nfr/nshift-i/nshift);
294 cacf[i] /= rfr;
295
296 if (time[vfr[i]] <= time[vfr[ei]])
297 {
298 sigma_ret = sigma;
299 }
300
301 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
302
303 if ((i+1) < nfr)
304 {
305 deltat = (time[vfr[i+1]]-time[vfr[i]]);
306 }
307 corint = 2.0*deltat*cacf[i]*prefactor;
308 if (i == 0 || (i+1) == nfr)
309 {
310 corint *= 0.5;
311 }
312 sigma += corint;
313
314 i++;
315 }
316
317 }
318 else
319 {
320 printf("Too less points.\n");
321 }
322
323 return sigma_ret;
324
325}
326
327static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
328{
329
330 int i;
331 real rtmp;
332 real rfr;
333
334
335 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
336
337
338
339 for (i = 0; i < nfr; i++)
340 {
341
342 if (refr[i] != 0.0)
343 {
344 dsp2[i] *= prefactor/refr[i];
345 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
346 }
347
348
349 }
350
351}
352
353
354static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
355 FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
356 int ePBC, t_topology top, t_trxframe fr, real temp,
357 real trust, real bfit, real efit, real bvit, real evit,
358 t_trxstatus *status, int isize, int nmols, int nshift,
359 atom_id *index0, int indexm[], real mass2[],
360 real qmol[], real eps_rf, const output_env_t oenv)
361{
362 int i, j, k, l, f;
363 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
364 int vshfr;
365 real *xshfr = NULL((void*)0);
366 int *vfr = NULL((void*)0);
367 real refr = 0.0;
368 real rfr = 0.0;
369 real *cacf = NULL((void*)0);
370 real *time = NULL((void*)0);
371 real *djc = NULL((void*)0);
372 real corint = 0.0;
373 real prefactorav = 0.0;
374 real prefactor = 0.0;
375 real volume;
376 real volume_av = 0.0;
377 real dk_s, dk_d, dk_f;
378 real dm_s, dm_d;
379 real mj = 0.0;
380 real mj2 = 0.0;
381 real mjd = 0.0;
382 real mjdav = 0.0;
383 real md2 = 0.0;
384 real mdav2 = 0.0;
385 real sgk;
386 rvec mja_tmp;
387 rvec mjd_tmp;
388 rvec mdvec;
389 rvec *mu = NULL((void*)0);
390 rvec *xp = NULL((void*)0);
391 rvec *v0 = NULL((void*)0);
392 rvec *mjdsp = NULL((void*)0);
393 real *dsp2 = NULL((void*)0);
394 real t0;
395 real rtmp;
396 real qtmp;
397
398
399
400 rvec tmp;
401 rvec *mtrans = NULL((void*)0);
402
403 /*
404 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
405 */
406
407 int bi, ei, ie, ii;
408 real rest = 0.0;
409 real sigma = 0.0;
410 real malt = 0.0;
411 real err = 0.0;
412 real *xfit;
413 real *yfit;
414 gmx_rmpbc_t gpbc = NULL((void*)0);
415
416 /*
417 * indices for EH
418 */
419
420 ei = 0;
421 bi = 0;
422
423 /*
424 * indices for GK
425 */
426
427 ii = 0;
428 ie = 0;
429 t0 = 0;
430 sgk = 0.0;
431
432
433 /* This is the main loop over frames */
434
435
436 nfr = 0;
437
438
439 nvfr = 0;
440 vshfr = 0;
441 nalloc = 0;
442 valloc = 0;
443
444 clear_rvec(mja_tmp);
445 clear_rvec(mjd_tmp);
446 clear_rvec(mdvec);
447 clear_rvec(tmp);
448 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
449
450 do
451 {
452
453 refr = (real)(nfr+1);
454
455 if (nfr >= nalloc)
456 {
457 nalloc += 100;
458 srenew(time, nalloc)(time) = save_realloc("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 458, (time), (nalloc), sizeof(*(time)));
459 srenew(mu, nalloc)(mu) = save_realloc("mu", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 459, (mu), (nalloc), sizeof(*(mu)));
460 srenew(mjdsp, nalloc)(mjdsp) = save_realloc("mjdsp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 460, (mjdsp), (nalloc), sizeof(*(mjdsp)));
461 srenew(dsp2, nalloc)(dsp2) = save_realloc("dsp2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 461, (dsp2), (nalloc), sizeof(*(dsp2)));
462 srenew(mtrans, nalloc)(mtrans) = save_realloc("mtrans", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 462, (mtrans), (nalloc), sizeof(*(mtrans)));
463 srenew(xshfr, nalloc)(xshfr) = save_realloc("xshfr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 463, (xshfr), (nalloc), sizeof(*(xshfr)));
464
465 for (i = nfr; i < nalloc; i++)
466 {
467 clear_rvec(mjdsp[i]);
468 clear_rvec(mu[i]);
469 clear_rvec(mtrans[i]);
470 dsp2[i] = 0.0;
471 xshfr[i] = 0.0;
472 }
473 }
474
475
476
477 if (nfr == 0)
478 {
479 t0 = fr.time;
480
481 }
482
483 time[nfr] = fr.time-t0;
484
485 if (time[nfr] <= bfit)
486 {
487 bi = nfr;
488 }
489 if (time[nfr] <= efit)
490 {
491 ei = nfr;
492 }
493
494 if (bNoJump)
495 {
496
497 if (xp)
498 {
499 remove_jump(fr.box, fr.natoms, xp, fr.x);
500 }
501 else
502 {
503 snew(xp, fr.natoms)(xp) = save_calloc("xp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 503, (fr.natoms), sizeof(*(xp)));
504 }
505
506 for (i = 0; i < fr.natoms; i++)
507 {
508 copy_rvec(fr.x[i], xp[i]);
509 }
510
511 }
512
513 gmx_rmpbc_trxfr(gpbc, &fr);
514
515 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
516
517 for (i = 0; i < isize; i++)
518 {
519 j = index0[i];
520 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
521 rvec_inc(mu[nfr], fr.x[j]);
522 }
523
524 /*if(mod(nfr,nshift)==0){*/
525 if (nfr%nshift == 0)
526 {
527 for (j = nfr; j >= 0; j--)
528 {
529 rvec_sub(mtrans[nfr], mtrans[j], tmp);
530 dsp2[nfr-j] += norm2(tmp);
531 xshfr[nfr-j] += 1.0;
532 }
533 }
534
535 if (fr.bV)
536 {
537 if (nvfr >= valloc)
538 {
539 valloc += 100;
540 srenew(vfr, valloc)(vfr) = save_realloc("vfr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 540, (vfr), (valloc), sizeof(*(vfr)));
541 if (bINT)
542 {
543 srenew(djc, valloc)(djc) = save_realloc("djc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 543, (djc), (valloc), sizeof(*(djc)));
544 }
545 srenew(v0, valloc)(v0) = save_realloc("v0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 545, (v0), (valloc), sizeof(*(v0)));
546 if (bACF)
547 {
548 srenew(cacf, valloc)(cacf) = save_realloc("cacf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 548, (cacf), (valloc), sizeof(*(cacf)));
549 }
550 }
551 if (time[nfr] <= bvit)
552 {
553 ii = nvfr;
554 }
555 if (time[nfr] <= evit)
556 {
557 ie = nvfr;
558 }
559 vfr[nvfr] = nfr;
560 clear_rvec(v0[nvfr]);
561 if (bACF)
562 {
563 cacf[nvfr] = 0.0;
564 }
565 if (bINT)
566 {
567 djc[nvfr] = 0.0;
568 }
569 for (i = 0; i < isize; i++)
570 {
571 j = index0[i];
572 svmul(mass2[j], fr.v[j], fr.v[j]);
573 svmul(qmol[j], fr.v[j], fr.v[j]);
574 rvec_inc(v0[nvfr], fr.v[j]);
575 }
576
577 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX0], v0[nfr][YY1], v0[nfr][ZZ2]);
578 if (bACF || bINT)
579 {
580 /*if(mod(nvfr,nshift)==0){*/
581 if (nvfr%nshift == 0)
582 {
583 for (j = nvfr; j >= 0; j--)
584 {
585 if (bACF)
586 {
587 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
588 }
589 if (bINT)
590 {
591 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
592 }
593 }
594 vshfr++;
595 }
596 }
597 nvfr++;
598 }
599
600 volume = det(fr.box);
601 volume_av += volume;
602
603 rvec_inc(mja_tmp, mtrans[nfr]);
604 mjd += iprod(mu[nfr], mtrans[nfr]);
605 rvec_inc(mdvec, mu[nfr]);
606
607 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
608 md2 += iprod(mu[nfr], mu[nfr]);
609
610 fprintf(fmj, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mtrans[nfr][XX0], mtrans[nfr][YY1], mtrans[nfr][ZZ2], mj2/refr, norm(mja_tmp)/refr);
611 fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
612 time[nfr], mu[nfr][XX0], mu[nfr][YY1], mu[nfr][ZZ2], md2/refr, norm(mdvec)/refr);
613
614 nfr++;
615
616 }
617 while (read_next_frame(oenv, status, &fr));
618
619 gmx_rmpbc_done(gpbc);
620
621 volume_av /= refr;
622
623 prefactor = 1.0;
624 prefactor /= 3.0*EPSILON0(5.72765E-4)*volume_av*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp;
625
626
627 prefactorav = E_CHARGE(1.60217733e-19)*E_CHARGE(1.60217733e-19);
628 prefactorav /= volume_av*BOLTZMANN(1.380658e-23)*temp*NANO(1e-9)*6.0;
629
630 fprintf(stderrstderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
631
632 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
633
634 /*
635 * Now we can average and calculate the correlation functions
636 */
637
638
639 mj2 /= refr;
640 mjd /= refr;
641 md2 /= refr;
642
643 svmul(1.0/refr, mdvec, mdvec);
644 svmul(1.0/refr, mja_tmp, mja_tmp);
645
646 mdav2 = norm2(mdvec);
647 mj = norm2(mja_tmp);
648 mjdav = iprod(mdvec, mja_tmp);
649
650
651 printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mja_tmp[XX0], mja_tmp[YY1], mja_tmp[ZZ2], mj2);
652 printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mdvec[XX0], mdvec[YY1], mdvec[ZZ2], md2);
653
654 if (v0 != NULL((void*)0))
655 {
656 trust *= (real) nvfr;
657 itrust = (int) trust;
658
659 if (bINT)
660 {
661
662 printf("\nCalculating M_D - current correlation integral ... \n");
663 corint = calc_cacf(mcor, prefactorav/EPSI0((5.72765E-4)*(1.60217733e-19)*(1.60217733e-19)*(6.0221367e23
)/((1e3)*(1e-9))), djc, time, nvfr, vfr, ie, nshift);
664
665 }
666
667 if (bACF)
668 {
669
670 printf("\nCalculating current autocorrelation ... \n");
671 sgk = calc_cacf(outf, prefactorav/PICO(1e-12), cacf, time, nvfr, vfr, ie, nshift);
672
673 if (ie > ii)
674 {
675
676 snew(xfit, ie-ii+1)(xfit) = save_calloc("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 676, (ie-ii+1), sizeof(*(xfit)));
677 snew(yfit, ie-ii+1)(yfit) = save_calloc("yfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 677, (ie-ii+1), sizeof(*(yfit)));
678
679 for (i = ii; i <= ie; i++)
680 {
681
682 xfit[i-ii] = log(time[vfr[i]]);
683 rtmp = fabs(cacf[i]);
684 yfit[i-ii] = log(rtmp);
685
686 }
687
688 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
689
690 malt = exp(malt);
691
692 sigma += 1.0;
693
694 malt *= prefactorav*2.0e12/sigma;
695
696 sfree(xfit)save_free("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 696, (xfit));
697 sfree(yfit)save_free("yfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 697, (yfit));
698
699 }
700 }
701 }
702
703
704 /* Calculation of the dielectric constant */
705
706 fprintf(stderrstderr, "\n********************************************\n");
707 dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE0);
708 fprintf(stderrstderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
709 fprintf(stderrstderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
710 fprintf(stderrstderr, "********************************************\n");
711
712
713 dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE0);
714 fprintf(stderrstderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
715 fprintf(stderrstderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
716 fprintf(stderrstderr, "\n********************************************\n");
717 if (bINT)
718 {
719 dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE1);
720 fprintf(stderrstderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
721 fprintf(stderrstderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
722 }
723
724 fprintf(stderrstderr, "\n***************************************************");
725 fprintf(stderrstderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
726 fprintf(stderrstderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
727
728
729
730 if (bACF)
731 {
732 fprintf(stderrstderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
733 fprintf(stderrstderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
734 }
735
736 if (ei > bi)
737 {
738 fprintf(stderrstderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
739 fprintf(stderrstderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
740
741 snew(xfit, ei-bi+1)(xfit) = save_calloc("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 741, (ei-bi+1), sizeof(*(xfit)));
742 snew(yfit, ei-bi+1)(yfit) = save_calloc("yfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 742, (ei-bi+1), sizeof(*(yfit)));
743
744 for (i = bi; i <= ei; i++)
745 {
746 xfit[i-bi] = time[i];
747 yfit[i-bi] = dsp2[i];
748 }
749
750 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
751
752 sigma *= 1e12;
753 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE1);
754
755 fprintf(stderrstderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
756 fprintf(stderrstderr, "sigma=%.4f\n", sigma);
757 fprintf(stderrstderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
758 fprintf(stderrstderr, "Dielectric constant using EH: %.4f\n", dk_d);
759
760 sfree(xfit)save_free("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 760, (xfit));
761 sfree(yfit)save_free("yfit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 761, (yfit));
762 }
763 else
764 {
765 fprintf(stderrstderr, "Too less points for a fit.\n");
766 }
767
768
769 if (v0 != NULL((void*)0))
770 {
771 sfree(v0)save_free("v0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 771, (v0));
772 }
773 if (bACF)
774 {
775 sfree(cacf)save_free("cacf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 775, (cacf));
776 }
777 if (bINT)
778 {
779 sfree(djc)save_free("djc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 779, (djc));
780 }
781
782 sfree(time)save_free("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 782, (time));
783
784
785 sfree(mjdsp)save_free("mjdsp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 785, (mjdsp));
786 sfree(mu)save_free("mu", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 786, (mu));
787}
788
789
790
791int gmx_current(int argc, char *argv[])
792{
793
794 static int nshift = 1000;
795 static real temp = 300.0;
796 static real trust = 0.25;
797 static real eps_rf = 0.0;
798 static gmx_bool bNoJump = TRUE1;
799 static real bfit = 100.0;
800 static real bvit = 0.5;
801 static real efit = 400.0;
802 static real evit = 5.0;
803 t_pargs pa[] = {
804 { "-sh", FALSE0, etINT, {&nshift},
805 "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
806 { "-nojump", FALSE0, etBOOL, {&bNoJump},
807 "Removes jumps of atoms across the box."},
808 { "-eps", FALSE0, etREAL, {&eps_rf},
809 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
810 { "-bfit", FALSE0, etREAL, {&bfit},
811 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
812 { "-efit", FALSE0, etREAL, {&efit},
813 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
814 { "-bvit", FALSE0, etREAL, {&bvit},
815 "Begin of the fit of the current autocorrelation function to a*t^b."},
816 { "-evit", FALSE0, etREAL, {&evit},
817 "End of the fit of the current autocorrelation function to a*t^b."},
818 { "-tr", FALSE0, etREAL, {&trust},
819 "Fraction of the trajectory taken into account for the integral."},
820 { "-temp", FALSE0, etREAL, {&temp},
821 "Temperature for calculating epsilon."}
822 };
823
824 output_env_t oenv;
825 t_topology top;
826 char title[STRLEN4096];
827 char **grpname = NULL((void*)0);
828 const char *indexfn;
829 t_trxframe fr;
830 real *mass2 = NULL((void*)0);
831 rvec *xtop, *vtop;
832 matrix box;
833 atom_id *index0;
834 int *indexm = NULL((void*)0);
835 int isize;
836 t_trxstatus *status;
837 int flags = 0;
838 gmx_bool bTop;
839 gmx_bool bNEU;
840 gmx_bool bACF;
841 gmx_bool bINT;
842 int ePBC = -1;
843 int natoms;
844 int nmols;
845 int i, j, k = 0, l;
846 int step;
847 real t;
848 real lambda;
849 real *qmol;
850 FILE *outf = NULL((void*)0);
851 FILE *outi = NULL((void*)0);
852 FILE *tfile = NULL((void*)0);
853 FILE *mcor = NULL((void*)0);
854 FILE *fmj = NULL((void*)0);
855 FILE *fmd = NULL((void*)0);
856 FILE *fmjdsp = NULL((void*)0);
857 FILE *fcur = NULL((void*)0);
858 t_filenm fnm[] = {
859 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, /* this is for the topology */
860 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
861 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, /* and this for the trajectory */
862 { efXVG, "-o", "current.xvg", ffWRITE1<<2 },
863 { efXVG, "-caf", "caf.xvg", ffOPTWR(1<<2| 1<<3) },
864 { efXVG, "-dsp", "dsp.xvg", ffWRITE1<<2 },
865 { efXVG, "-md", "md.xvg", ffWRITE1<<2 },
866 { efXVG, "-mj", "mj.xvg", ffWRITE1<<2},
867 { efXVG, "-mc", "mc.xvg", ffOPTWR(1<<2| 1<<3) }
868 };
869
870#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
871
872
873 const char *desc[] = {
874 "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
875 "of the rotational and translational dipole moment of the system, and the resulting static",
876 "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
877 "Furthermore, the routine is capable of extracting the static conductivity from the current ",
878 "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
879 "can be used to obtain the static conductivity."
880 "[PAR]",
881 "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
882 "correlation of the rotational and translational part of the dipole moment in the corresponding",
883 "file. However, this option is only available for trajectories containing velocities.",
884 "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
885 "autocorrelation functions. Since averaging proceeds by shifting the starting point",
886 "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
887 "starting points. Towards the end, statistical inaccuracy grows and integrating the",
888 "correlation function only yields reliable values until a certain point, depending on",
889 "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
890 "for calculating the static dielectric constant.",
891 "[PAR]",
892 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
893 "[PAR]",
894 "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
895 "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
896 "tin-foil boundary conditions).",
897 "[PAR]",
898 "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
899 "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
900 "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
901 "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
902 "option has to be used with care, since only very short time spans fulfill the approximation that the density",
903 "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
904 "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
905 "the translational part of the dielectric constant."
906 };
907
908
909 /* At first the arguments will be parsed and the system information processed */
910 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5),
911 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
912 {
913 return 0;
914 }
915
916 bACF = opt2bSet("-caf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
917 bINT = opt2bSet("-mc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
918
919 bTop = read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE1);
920
921
922
923 sfree(xtop)save_free("xtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 923, (xtop));
924 sfree(vtop)save_free("vtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 924, (vtop));
925 indexfn = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
926 snew(grpname, 1)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 926, (1), sizeof(*(grpname)));
927
928 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
929
930 flags = flags | TRX_READ_X(1<<0) | TRX_READ_V(1<<2);
931
932 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &fr, flags);
933
934 snew(mass2, top.atoms.nr)(mass2) = save_calloc("mass2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 934, (top.atoms.nr), sizeof(*(mass2)));
935 snew(qmol, top.atoms.nr)(qmol) = save_calloc("qmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 935, (top.atoms.nr), sizeof(*(qmol)));
936
937 bNEU = precalc(top, mass2, qmol);
938
939
940 snew(indexm, isize)(indexm) = save_calloc("indexm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_current.c"
, 940, (isize), sizeof(*(indexm)));
941
942 for (i = 0; i < isize; i++)
943 {
944 indexm[i] = index0[i];
945 }
946
947 nmols = isize;
948
949
950 index_atom2mol(&nmols, indexm, &top.mols);
951
952 if (fr.bV)
953 {
954 if (bACF)
955 {
956 outf = xvgropen(opt2fn("-caf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
957 "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
958 "ACF (e nm/ps)\\S2", oenv);
959 fprintf(outf, "# time\t acf\t average \t std.dev\n");
960 }
961 fcur = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
962 "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
963 fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
964 if (bINT)
965 {
966 mcor = xvgropen(opt2fn("-mc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
967 "M\\sD\\N - current autocorrelation function",
968 output_env_get_xvgr_tlabel(oenv),
969 "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
970 fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
971 }
972 }
973
974 fmj = xvgropen(opt2fn("-mj", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
975 "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
976 "< M\\sJ\\N > (enm)", oenv);
977 fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
978 fmd = xvgropen(opt2fn("-md", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
979 "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
980 "< M\\sD\\N > (enm)", oenv);
981 fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
982
983 fmjdsp = xvgropen(opt2fn("-dsp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
984 "MSD of the squared translational dipole moment M",
985 output_env_get_xvgr_tlabel(oenv),
986 "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
987 oenv);
988
989
990 /* System information is read and prepared, dielectric() processes the frames
991 * and calculates the requested quantities */
992
993 dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
994 temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
995 index0, indexm, mass2, qmol, eps_rf, oenv);
996
997 gmx_ffclose(fmj);
998 gmx_ffclose(fmd);
999 gmx_ffclose(fmjdsp);
1000 if (bACF)
1001 {
1002 gmx_ffclose(outf);
1003 }
1004 if (bINT)
1005 {
1006 gmx_ffclose(mcor);
1007 }
1008
1009 return 0;
1010}