?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Add replacements for pbc enumerations
[gromacs.git] / src / gromacs / gmxana / gmx_current.cpp
blobfd6c082ae809e82fc72ccc0388457275b79c8993
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2017,2018, 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 #include "gmxpre.h"
36
37 #include <cmath>
38 #include <cstdlib>
39
40 #include "gromacs/commandline/pargs.h"
41 #include "gromacs/fileio/confio.h"
42 #include "gromacs/fileio/trxio.h"
43 #include "gromacs/fileio/xvgr.h"
44 #include "gromacs/gmxana/gmx_ana.h"
45 #include "gromacs/math/units.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/pbcutil/pbc.h"
48 #include "gromacs/pbcutil/rmpbc.h"
49 #include "gromacs/statistics/statistics.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/topology/topology.h"
52 #include "gromacs/trajectory/trajectoryframe.h"
53 #include "gromacs/utility/arraysize.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/gmxassert.h"
57 #include "gromacs/utility/smalloc.h"
58
59 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
60
61 static 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(FARGS, "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(FARGS, "The index group does not consist of whole molecules");
84 }
85 i++;
86 }
87 index[nmol++] = mol;
88 }
89
90 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
91
92 *n = nmol;
93 }
94
95 static 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 (std::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 = FALSE;
147 }
148 else
149 {
150 bNEU = TRUE;
151 }
152
153 return bNEU;
154
155 }
156
157 static 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 < DIM; d++)
164 {
165 hbox[d] = 0.5*box[d][d];
166 }
167 for (i = 0; i < natoms; i++)
168 {
169 for (m = DIM-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
189 static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, const 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
239 static 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
275 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, const int vfr[], int ei, int nshift)
276 {
277
278 int i;
279 real deltat = 0.0;
280 real rfr;
281 real sigma = 0.0;
282 real sigma_ret = 0.0;
283
284 if (nfr > 1)
285 {
286 i = 0;
287
288 while (i < nfr)
289 {
290 real corint;
291
292 rfr = static_cast<real>(nfr+i)/nshift;
293 cacf[i] /= rfr;
294
295 if (time[vfr[i]] <= time[vfr[ei]])
296 {
297 sigma_ret = sigma;
298 }
299
300 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
301
302 if ((i+1) < nfr)
303 {
304 deltat = (time[vfr[i+1]]-time[vfr[i]]);
305 }
306 corint = 2.0*deltat*cacf[i]*prefactor;
307 if (i == 0 || (i+1) == nfr)
308 {
309 corint *= 0.5;
310 }
311 sigma += corint;
312
313 i++;
314 }
315
316 }
317 else
318 {
319 printf("Too less points.\n");
320 }
321
322 return sigma_ret;
323
324 }
325
326 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, const real refr[])
327 {
328
329 int i;
330
331 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
332
333 for (i = 0; i < nfr; i++)
334 {
335
336 if (refr[i] != 0.0)
337 {
338 dsp2[i] *= prefactor/refr[i];
339 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
340 }
341
342
343 }
344
345 }
346
347
348 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
349 FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
350 int ePBC, t_topology top, t_trxframe fr, real temp,
351 real bfit, real efit, real bvit, real evit,
352 t_trxstatus *status, int isize, int nmols, int nshift,
353 const int *index0, int indexm[], real mass2[],
354 real qmol[], real eps_rf, const gmx_output_env_t *oenv)
355 {
356 int i, j;
357 int valloc, nalloc, nfr, nvfr;
358 int vshfr;
359 real *xshfr = nullptr;
360 int *vfr = nullptr;
361 real refr = 0.0;
362 real *cacf = nullptr;
363 real *time = nullptr;
364 real *djc = nullptr;
365 real corint = 0.0;
366 real prefactorav = 0.0;
367 real prefactor = 0.0;
368 real volume;
369 real volume_av = 0.0;
370 real dk_s, dk_d, dk_f;
371 real mj = 0.0;
372 real mj2 = 0.0;
373 real mjd = 0.0;
374 real mjdav = 0.0;
375 real md2 = 0.0;
376 real mdav2 = 0.0;
377 real sgk;
378 rvec mja_tmp;
379 rvec mjd_tmp;
380 rvec mdvec;
381 rvec *mu = nullptr;
382 rvec *xp = nullptr;
383 rvec *v0 = nullptr;
384 rvec *mjdsp = nullptr;
385 real *dsp2 = nullptr;
386 real t0;
387 real rtmp;
388
389 rvec tmp;
390 rvec *mtrans = nullptr;
391
392 /*
393 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
394 */
395
396 int bi, ei, ie, ii;
397 real rest = 0.0;
398 real sigma = 0.0;
399 real malt = 0.0;
400 real err = 0.0;
401 real *xfit;
402 real *yfit;
403 gmx_rmpbc_t gpbc = nullptr;
404
405 /*
406 * indices for EH
407 */
408
409 ei = 0;
410 bi = 0;
411
412 /*
413 * indices for GK
414 */
415
416 ii = 0;
417 ie = 0;
418 t0 = 0;
419 sgk = 0.0;
420
421
422 /* This is the main loop over frames */
423
424
425 nfr = 0;
426
427
428 nvfr = 0;
429 vshfr = 0;
430 nalloc = 0;
431 valloc = 0;
432
433 clear_rvec(mja_tmp);
434 clear_rvec(mjd_tmp);
435 clear_rvec(mdvec);
436 clear_rvec(tmp);
437 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
438
439 do
440 {
441
442 refr = (nfr+1);
443
444 if (nfr >= nalloc)
445 {
446 nalloc += 100;
447 srenew(time, nalloc);
448 srenew(mu, nalloc);
449 srenew(mjdsp, nalloc);
450 srenew(dsp2, nalloc);
451 srenew(mtrans, nalloc);
452 srenew(xshfr, nalloc);
453
454 for (i = nfr; i < nalloc; i++)
455 {
456 clear_rvec(mjdsp[i]);
457 clear_rvec(mu[i]);
458 clear_rvec(mtrans[i]);
459 dsp2[i] = 0.0;
460 xshfr[i] = 0.0;
461 }
462 }
463 GMX_RELEASE_ASSERT(time != nullptr, "Memory not allocated correctly - time array is NULL");
464
465 if (nfr == 0)
466 {
467 t0 = fr.time;
468
469 }
470
471 time[nfr] = fr.time-t0;
472
473 if (time[nfr] <= bfit)
474 {
475 bi = nfr;
476 }
477 if (time[nfr] <= efit)
478 {
479 ei = nfr;
480 }
481
482 if (bNoJump)
483 {
484
485 if (xp)
486 {
487 remove_jump(fr.box, fr.natoms, xp, fr.x);
488 }
489 else
490 {
491 snew(xp, fr.natoms);
492 }
493
494 for (i = 0; i < fr.natoms; i++)
495 {
496 copy_rvec(fr.x[i], xp[i]);
497 }
498
499 }
500
501 gmx_rmpbc_trxfr(gpbc, &fr);
502
503 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
504
505 for (i = 0; i < isize; i++)
506 {
507 j = index0[i];
508 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
509 rvec_inc(mu[nfr], fr.x[j]);
510 }
511
512 /*if(mod(nfr,nshift)==0){*/
513 if (nfr%nshift == 0)
514 {
515 for (j = nfr; j >= 0; j--)
516 {
517 rvec_sub(mtrans[nfr], mtrans[j], tmp);
518 dsp2[nfr-j] += norm2(tmp);
519 xshfr[nfr-j] += 1.0;
520 }
521 }
522
523 if (fr.bV)
524 {
525 if (nvfr >= valloc)
526 {
527 valloc += 100;
528 srenew(vfr, valloc);
529 if (bINT)
530 {
531 srenew(djc, valloc);
532 }
533 srenew(v0, valloc);
534 if (bACF)
535 {
536 srenew(cacf, valloc);
537 }
538 }
539 if (time[nfr] <= bvit)
540 {
541 ii = nvfr;
542 }
543 if (time[nfr] <= evit)
544 {
545 ie = nvfr;
546 }
547 vfr[nvfr] = nfr;
548 clear_rvec(v0[nvfr]);
549 if (bACF)
550 {
551 cacf[nvfr] = 0.0;
552 }
553 if (bINT)
554 {
555 djc[nvfr] = 0.0;
556 }
557 for (i = 0; i < isize; i++)
558 {
559 j = index0[i];
560 svmul(mass2[j], fr.v[j], fr.v[j]);
561 svmul(qmol[j], fr.v[j], fr.v[j]);
562 rvec_inc(v0[nvfr], fr.v[j]);
563 }
564
565 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
566 if (bACF || bINT)
567 {
568 /*if(mod(nvfr,nshift)==0){*/
569 if (nvfr%nshift == 0)
570 {
571 for (j = nvfr; j >= 0; j--)
572 {
573 if (bACF)
574 {
575 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
576 }
577 if (bINT)
578 {
579 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
580 }
581 }
582 vshfr++;
583 }
584 }
585 nvfr++;
586 }
587
588 volume = det(fr.box);
589 volume_av += volume;
590
591 rvec_inc(mja_tmp, mtrans[nfr]);
592 mjd += iprod(mu[nfr], mtrans[nfr]);
593 rvec_inc(mdvec, mu[nfr]);
594
595 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
596 md2 += iprod(mu[nfr], mu[nfr]);
597
598 fprintf(fmj, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mtrans[nfr][XX], mtrans[nfr][YY], mtrans[nfr][ZZ], mj2/refr, norm(mja_tmp)/refr);
599 fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
600 time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
601
602 nfr++;
603
604 }
605 while (read_next_frame(oenv, status, &fr));
606
607 gmx_rmpbc_done(gpbc);
608
609 volume_av /= refr;
610
611 prefactor = 1.0;
612 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
613
614
615 prefactorav = E_CHARGE*E_CHARGE;
616 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
617
618 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
619
620 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
621
622 /*
623 * Now we can average and calculate the correlation functions
624 */
625
626
627 mj2 /= refr;
628 mjd /= refr;
629 md2 /= refr;
630
631 svmul(1.0/refr, mdvec, mdvec);
632 svmul(1.0/refr, mja_tmp, mja_tmp);
633
634 mdav2 = norm2(mdvec);
635 mj = norm2(mja_tmp);
636 mjdav = iprod(mdvec, mja_tmp);
637
638
639 printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mja_tmp[XX], mja_tmp[YY], mja_tmp[ZZ], mj2);
640 printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n", nfr, mdvec[XX], mdvec[YY], mdvec[ZZ], md2);
641
642 if (v0 != nullptr)
643 {
644 if (bINT)
645 {
646
647 printf("\nCalculating M_D - current correlation integral ... \n");
648 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
649
650 }
651
652 if (bACF)
653 {
654
655 printf("\nCalculating current autocorrelation ... \n");
656 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
657
658 if (ie > ii)
659 {
660
661 snew(xfit, ie-ii+1);
662 snew(yfit, ie-ii+1);
663
664 for (i = ii; i <= ie; i++)
665 {
666
667 xfit[i-ii] = std::log(time[vfr[i]]);
668 rtmp = std::abs(cacf[i]);
669 yfit[i-ii] = std::log(rtmp);
670
671 }
672
673 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
674
675 malt = std::exp(malt);
676
677 sigma += 1.0;
678
679 malt *= prefactorav*2.0e12/sigma;
680
681 sfree(xfit);
682 sfree(yfit);
683
684 }
685 }
686 }
687
688
689 /* Calculation of the dielectric constant */
690
691 fprintf(stderr, "\n********************************************\n");
692 dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
693 fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
694 fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
695 fprintf(stderr, "********************************************\n");
696
697
698 dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
699 fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
700 fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
701 fprintf(stderr, "\n********************************************\n");
702 if (bINT)
703 {
704 dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
705 fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
706 fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
707 }
708
709 fprintf(stderr, "\n***************************************************");
710 fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
711 fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
712
713
714
715 if (bACF && (ii < nvfr))
716 {
717 fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
718 fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*std::pow(time[vfr[ii]], sigma), sgk);
719 }
720
721 if (ei > bi)
722 {
723 fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
724 fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
725
726 snew(xfit, ei-bi+1);
727 snew(yfit, ei-bi+1);
728
729 for (i = bi; i <= ei; i++)
730 {
731 xfit[i-bi] = time[i];
732 yfit[i-bi] = dsp2[i];
733 }
734
735 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
736
737 sigma *= 1e12;
738 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
739
740 fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
741 fprintf(stderr, "sigma=%.4f\n", sigma);
742 fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
743 fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
744
745 sfree(xfit);
746 sfree(yfit);
747 }
748 else
749 {
750 fprintf(stderr, "Too few points for a fit.\n");
751 }
752
753
754 if (v0 != nullptr)
755 {
756 sfree(v0);
757 }
758 if (bACF)
759 {
760 sfree(cacf);
761 }
762 if (bINT)
763 {
764 sfree(djc);
765 }
766
767 sfree(time);
768
769
770 sfree(mjdsp);
771 sfree(mu);
772 }
773
774
775
776 int gmx_current(int argc, char *argv[])
777 {
778
779 static int nshift = 1000;
780 static real temp = 300.0;
781 static real eps_rf = 0.0;
782 static gmx_bool bNoJump = TRUE;
783 static real bfit = 100.0;
784 static real bvit = 0.5;
785 static real efit = 400.0;
786 static real evit = 5.0;
787 t_pargs pa[] = {
788 { "-sh", FALSE, etINT, {&nshift},
789 "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
790 { "-nojump", FALSE, etBOOL, {&bNoJump},
791 "Removes jumps of atoms across the box."},
792 { "-eps", FALSE, etREAL, {&eps_rf},
793 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
794 { "-bfit", FALSE, etREAL, {&bfit},
795 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
796 { "-efit", FALSE, etREAL, {&efit},
797 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
798 { "-bvit", FALSE, etREAL, {&bvit},
799 "Begin of the fit of the current autocorrelation function to a*t^b."},
800 { "-evit", FALSE, etREAL, {&evit},
801 "End of the fit of the current autocorrelation function to a*t^b."},
802 { "-temp", FALSE, etREAL, {&temp},
803 "Temperature for calculating epsilon."}
804 };
805
806 gmx_output_env_t *oenv;
807 t_topology top;
808 char **grpname = nullptr;
809 const char *indexfn;
810 t_trxframe fr;
811 real *mass2 = nullptr;
812 matrix box;
813 int *index0;
814 int *indexm = nullptr;
815 int isize;
816 t_trxstatus *status;
817 int flags = 0;
818 gmx_bool bACF;
819 gmx_bool bINT;
820 int ePBC = -1;
821 int nmols;
822 int i;
823 real *qmol;
824 FILE *outf = nullptr;
825 FILE *mcor = nullptr;
826 FILE *fmj = nullptr;
827 FILE *fmd = nullptr;
828 FILE *fmjdsp = nullptr;
829 FILE *fcur = nullptr;
830 t_filenm fnm[] = {
831 { efTPS, nullptr, nullptr, ffREAD }, /* this is for the topology */
832 { efNDX, nullptr, nullptr, ffOPTRD },
833 { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
834 { efXVG, "-o", "current", ffWRITE },
835 { efXVG, "-caf", "caf", ffOPTWR },
836 { efXVG, "-dsp", "dsp", ffWRITE },
837 { efXVG, "-md", "md", ffWRITE },
838 { efXVG, "-mj", "mj", ffWRITE },
839 { efXVG, "-mc", "mc", ffOPTWR }
840 };
841
842 #define NFILE asize(fnm)
843
844
845 const char *desc[] = {
846 "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
847 "of the rotational and translational dipole moment of the system, and the resulting static",
848 "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
849 "Furthermore, the routine is capable of extracting the static conductivity from the current ",
850 "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
851 "can be used to obtain the static conductivity.",
852 "[PAR]",
853 "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
854 "correlation of the rotational and translational part of the dipole moment in the corresponding",
855 "file. However, this option is only available for trajectories containing velocities.",
856 "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
857 "autocorrelation functions. Since averaging proceeds by shifting the starting point",
858 "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
859 "starting points. Towards the end, statistical inaccuracy grows and integrating the",
860 "correlation function only yields reliable values until a certain point, depending on",
861 "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
862 "for calculating the static dielectric constant.",
863 "[PAR]",
864 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
865 "[PAR]",
866 "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
867 "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]\\=0 corresponds to",
868 "tin-foil boundary conditions).",
869 "[PAR]",
870 "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
871 "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
872 "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
873 "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
874 "option has to be used with care, since only very short time spans fulfill the approximation that the density",
875 "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
876 "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
877 "the translational part of the dielectric constant."
878 };
879
880
881 /* At first the arguments will be parsed and the system information processed */
882 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
883 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
884 {
885 return 0;
886 }
887
888 bACF = opt2bSet("-caf", NFILE, fnm);
889 bINT = opt2bSet("-mc", NFILE, fnm);
890
891 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, nullptr, nullptr, box, TRUE);
892
893 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
894 snew(grpname, 1);
895
896 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
897
898 flags = flags | TRX_READ_X | TRX_READ_V;
899
900 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
901
902 snew(mass2, top.atoms.nr);
903 snew(qmol, top.atoms.nr);
904
905 precalc(top, mass2, qmol);
906
907
908 snew(indexm, isize);
909
910 for (i = 0; i < isize; i++)
911 {
912 indexm[i] = index0[i];
913 }
914
915 nmols = isize;
916
917
918 index_atom2mol(&nmols, indexm, &top.mols);
919
920 if (fr.bV)
921 {
922 if (bACF)
923 {
924 outf = xvgropen(opt2fn("-caf", NFILE, fnm),
925 "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
926 "ACF (e nm/ps)\\S2", oenv);
927 fprintf(outf, "# time\t acf\t average \t std.dev\n");
928 }
929 fcur = xvgropen(opt2fn("-o", NFILE, fnm),
930 "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
931 fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
932 if (bINT)
933 {
934 mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
935 "M\\sD\\N - current autocorrelation function",
936 output_env_get_xvgr_tlabel(oenv),
937 "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
938 fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
939 }
940 }
941
942 fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
943 "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
944 "< M\\sJ\\N > (enm)", oenv);
945 fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
946 fmd = xvgropen(opt2fn("-md", NFILE, fnm),
947 "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
948 "< M\\sD\\N > (enm)", oenv);
949 fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
950
951 fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
952 "MSD of the squared translational dipole moment M",
953 output_env_get_xvgr_tlabel(oenv),
954 "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
955 oenv);
956
957
958 /* System information is read and prepared, dielectric() processes the frames
959 * and calculates the requested quantities */
960
961 dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
962 temp, bfit, efit, bvit, evit, status, isize, nmols, nshift,
963 index0, indexm, mass2, qmol, eps_rf, oenv);
964
965 xvgrclose(fmj);
966 xvgrclose(fmd);
967 xvgrclose(fmjdsp);
968 if (fr.bV)
969 {
970 if (bACF)
971 {
972 xvgrclose(outf);
973 }
974 xvgrclose(fcur);
975 if (bINT)
976 {
977 xvgrclose(mcor);
978 }
979 }
980
981 return 0;
982 }
The ultimate molecular dynamics simulation package