2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/pbcutil/pbc.h"
46 #include "gromacs/utility/futil.h"
49 #include "gromacs/statistics/statistics.h"
51 #include "gromacs/topology/index.h"
52 #include "gromacs/random/random.h"
54 #include "gromacs/math/units.h"
56 #include "gromacs/fileio/enxio.h"
59 #include "gromacs/fileio/trxio.h"
61 #include "gromacs/commandline/pargs.h"
62 #include "gromacs/fileio/matio.h"
63 #include "gromacs/fileio/xvgr.h"
64 #include "gromacs/linearalgebra/nrjac.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/pbcutil/rmpbc.h"
67 #include "gromacs/utility/exceptions.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/smalloc.h"
71 #define e2d(x) ENM2DEBYE*(x)
72 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
73 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
85 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
93 if ((ptr = getenv("GMX_DIPOLE_SPACING")) != NULL)
95 double bw = strtod(ptr, NULL);
100 gb->spacing = 0.01; /* nm */
102 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
109 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
112 snew(gb->elem, gb->nelem);
113 snew(gb->count, gb->nelem);
115 snew(gb->cmap, gb->nx);
116 gb->ny = std::max(2, ndegrees);
117 for (i = 0; (i < gb->nx); i++)
119 snew(gb->cmap[i], gb->ny);
126 static void done_gkrbin(t_gkrbin **gb)
134 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
136 int cy, index = gmx_nint(r/gb->spacing);
139 if (index < gb->nelem)
141 gb->elem[index] += cosa;
149 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
153 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
155 cy = std::min(gb->ny-1, std::max(0, cy));
158 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
160 gb->cmap[index][cy] += 1;
164 static void rvec2sprvec(rvec dipcart, rvec dipsp)
167 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
168 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
169 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
174 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
175 int mindex[], rvec x[], rvec mu[],
176 int ePBC, matrix box, t_atom *atom, int *nAtom)
178 static rvec *xcm[2] = { NULL, NULL};
179 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
180 real qtot, q, cosa, rr, phi;
184 for (n = 0; (n < ncos); n++)
188 snew(xcm[n], ngrp[n]);
190 for (i = 0; (i < ngrp[n]); i++)
192 /* Calculate center of mass of molecule */
198 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
203 clear_rvec(xcm[n][i]);
205 for (j = j0; j < j1; j++)
209 for (k = 0; k < DIM; k++)
211 xcm[n][i][k] += q*x[j][k];
214 svmul(1/qtot, xcm[n][i], xcm[n][i]);
218 set_pbc(&pbc, ePBC, box);
221 for (i = 0; i < ngrp[grp0]; i++)
223 gi = molindex[grp0][i];
224 j0 = (ncos == 2) ? 0 : i+1;
225 for (j = j0; j < ngrp[grp1]; j++)
227 gj = molindex[grp1][j];
228 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
230 /* Compute distance between molecules including PBC */
231 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
237 rvec r_ij, r_kj, r_kl, mm, nn;
241 copy_rvec(xcm[grp0][i], xj);
242 copy_rvec(xcm[grp1][j], xk);
243 rvec_add(xj, mu[gi], xi);
244 rvec_add(xk, mu[gj], xl);
245 phi = dih_angle(xi, xj, xk, xl, &pbc,
246 r_ij, r_kj, r_kl, mm, nn, /* out */
247 &sign, &t1, &t2, &t3);
252 cosa = cos_angle(mu[gi], mu[gj]);
255 if (debug || gmx_isnan(cosa))
257 fprintf(debug ? debug : stderr,
258 "mu[%d] = %5.2f %5.2f %5.2f |mi| = %5.2f, mu[%d] = %5.2f %5.2f %5.2f |mj| = %5.2f rr = %5.2f cosa = %5.2f\n",
259 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
260 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
264 add2gkr(gb, rr, cosa, phi);
270 static real normalize_cmap(t_gkrbin *gb)
277 for (i = 0; (i < gb->nx); i++)
279 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
280 for (j = 0; (j < gb->ny); j++)
282 gb->cmap[i][j] /= vol;
283 hi = std::max(hi, gb->cmap[i][j]);
288 gmx_fatal(FARGS, "No data in the cmap");
293 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
300 t_rgb rlo = { 1, 1, 1 };
301 t_rgb rhi = { 0, 0, 0 };
303 hi = normalize_cmap(gb);
304 snew(xaxis, gb->nx+1);
305 for (i = 0; (i < gb->nx+1); i++)
307 xaxis[i] = i*gb->spacing;
310 for (j = 0; (j < gb->ny); j++)
314 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
318 yaxis[j] = (180.0*j)/(gb->ny-1.0);
320 /*2.0*j/(gb->ny-1.0)-1.0;*/
322 out = gmx_ffopen(cmap, "w");
324 "Dipole Orientation Distribution", "Fraction", "r (nm)",
325 gb->bPhi ? "Phi" : "Alpha",
326 gb->nx, gb->ny, xaxis, yaxis,
327 gb->cmap, 0, hi, rlo, rhi, nlevels);
333 static void print_gkrbin(const char *fn, t_gkrbin *gb,
334 int ngrp, int nframes, real volume,
335 const output_env_t oenv)
337 /* We compute Gk(r), gOO and hOO according to
338 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
339 * In this implementation the angle between dipoles is stored
340 * rather than their inner product. This allows to take polarizible
341 * models into account. The RDF is calculated as well, almost for free!
344 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
346 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
349 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
350 xvgr_legend(fp, asize(leg), leg, oenv);
352 Gkr = 1; /* Self-dipole inproduct = 1 */
357 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
358 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
362 while (last > 1 && gb->elem[last-1] == 0)
367 /* Divide by dipole squared, by number of frames, by number of origins.
368 * Multiply by 2 because we only take half the matrix of interactions
371 fac = 2.0/((double) ngrp * (double) nframes);
374 for (i = 0; i < last; i++)
376 /* Centre of the coordinate in the spherical layer */
379 /* Volume of the layer */
380 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
383 gOO = gb->count[i]*fac/(rho*vol_s);
385 /* Dipole correlation hOO, normalized by the relative number density, like
386 * in a Radial distribution function.
388 ggg = gb->elem[i]*fac;
389 hOO = 3.0*ggg/(rho*vol_s);
393 cosav = gb->elem[i]/gb->count[i];
399 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
401 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
402 x1, Gkr, cosav, hOO, gOO, ener);
410 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
411 real *t, int nre, t_enxframe *fr)
417 bCont = do_enx(fmu, fr);
420 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
421 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
426 if (Vol != -1) /* we've got Volume in the energy file */
428 *vol = fr->ener[Vol].e;
430 for (i = 0; i < DIM; i++)
432 mu[i] = fr->ener[iMu[i]].e;
440 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
443 int ncharged, m, a0, a1, a;
446 for (m = 0; m < n; m++)
448 a0 = mols->index[index[m]];
449 a1 = mols->index[index[m]+1];
452 for (a = a0; a < a1; a++)
457 /* This check is only for the count print */
458 if (fabs(qtot) > 0.01)
464 /* Remove the net charge at the center of mass */
465 for (a = a0; a < a1; a++)
467 atom[a].q -= qtot*atom[a].m/mtot;
474 printf("There are %d charged molecules in the selection,\n"
475 "will subtract their charge at their center of mass\n", ncharged);
479 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
485 for (k = k0; (k < k1); k++)
488 for (m = 0; (m < DIM); m++)
495 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
497 int i, k, m, n, niter;
498 real q, r2, mass, masstot;
499 rvec com; /* center of mass */
500 rvec r; /* distance of atoms to center of mass */
502 double dd[DIM], **ev, tmp;
506 for (i = 0; (i < DIM); i++)
513 /* Compute center of mass */
516 for (k = k0; (k < k1); k++)
520 for (i = 0; (i < DIM); i++)
522 com[i] += mass*x[k][i];
525 svmul((1.0/masstot), com, com);
527 /* We want traceless quadrupole moments, so let us calculate the complete
528 * quadrupole moment tensor and diagonalize this tensor to get
529 * the individual components on the diagonal.
532 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
534 for (m = 0; (m < DIM); m++)
536 for (n = 0; (n < DIM); n++)
541 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
543 q = (atom[k].q)*100.0;
544 rvec_sub(x[k], com, r);
546 for (m = 0; (m < DIM); m++)
548 for (n = 0; (n < DIM); n++)
550 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
556 for (i = 0; (i < DIM); i++)
558 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
559 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
563 /* We've got the quadrupole tensor, now diagonalize the sucker */
564 jacobi(inten, 3, dd, ev, &niter);
568 for (i = 0; (i < DIM); i++)
570 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
571 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
573 for (i = 0; (i < DIM); i++)
575 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
576 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
579 /* Sort the eigenvalues, for water we know that the order is as follows:
583 * At the moment I have no idea how this will work out for other molecules...
587 if (dd[i+1] > dd[i]) { \
596 for (m = 0; (m < DIM); m++)
598 quad[0] = dd[2]; /* yy */
599 quad[1] = dd[0]; /* zz */
600 quad[2] = dd[1]; /* xx */
605 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
609 for (i = 0; (i < DIM); i++)
619 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
621 real calc_eps(double M_diff, double volume, double epsRF, double temp)
623 double eps, A, teller, noemer;
624 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
625 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
627 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
636 teller = 1 + (A*2*epsRF/(2*epsRF+1));
637 noemer = 1 - (A/(2*epsRF+1));
639 eps = teller / noemer;
644 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
645 int idim, int nslice, rvec slab_dipole[],
651 for (k = k0; (k < k1); k++)
656 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
657 rvec_inc(slab_dipole[k], mu);
660 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
661 rvec slab_dipole[], matrix box, int nframes,
662 const output_env_t oenv)
668 const char *leg_dim[4] = {
669 "\\f{12}m\\f{4}\\sX\\N",
670 "\\f{12}m\\f{4}\\sY\\N",
671 "\\f{12}m\\f{4}\\sZ\\N",
672 "\\f{12}m\\f{4}\\stot\\N"
675 sprintf(buf, "Box-%c (nm)", 'X'+idim);
676 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
678 xvgr_legend(fp, DIM, leg_dim, oenv);
679 for (i = 0; (i < nslice); i++)
681 mutot = norm(slab_dipole[i])/nframes;
682 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
683 ((i+0.5)*box[idim][idim])/nslice,
684 slab_dipole[i][XX]/nframes,
685 slab_dipole[i][YY]/nframes,
686 slab_dipole[i][ZZ]/nframes,
690 do_view(oenv, fn, "-autoscale xy -nxy");
693 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
696 double dc, d, ddc1, ddc2, ddc3;
697 rvec xxx = { 1, 0, 0 };
698 rvec yyy = { 0, 1, 0 };
699 rvec zzz = { 0, 0, 1 };
702 ddc1 = ddc2 = ddc3 = 0;
703 for (i = k = 0; (i < n); i++)
705 ddc1 += fabs(cos_angle(dip[i], xxx));
706 ddc2 += fabs(cos_angle(dip[i], yyy));
707 ddc3 += fabs(cos_angle(dip[i], zzz));
710 for (j = i+1; (j < n); j++, k++)
712 dc = cos_angle(dip[i], dip[j]);
723 static void do_dip(t_topology *top, int ePBC, real volume,
725 const char *out_mtot, const char *out_eps,
726 const char *out_aver, const char *dipdist,
727 const char *cosaver, const char *fndip3d,
728 const char *fnadip, gmx_bool bPairs,
729 const char *corrtype, const char *corf,
730 gmx_bool bGkr, const char *gkrfn,
731 gmx_bool bPhi, int *nlevels, int ndegrees,
733 const char *cmap, real rcmax,
735 gmx_bool bMU, const char *mufn,
736 int *gnx, int *molindex[],
737 real mu_max, real mu_aver,
738 real epsilonRF, real temp,
739 int *gkatom, int skip,
740 gmx_bool bSlab, int nslices,
741 const char *axtitle, const char *slabfn,
742 const output_env_t oenv)
744 const char *leg_mtot[] = {
750 #define NLEGMTOT asize(leg_mtot)
751 const char *leg_eps[] = {
756 #define NLEGEPS asize(leg_eps)
757 const char *leg_aver[] = {
760 "< |M|\\S2\\N > - < |M| >\\S2\\N",
761 "< |M| >\\S2\\N / < |M|\\S2\\N >"
763 #define NLEGAVER asize(leg_aver)
764 const char *leg_cosaver[] = {
765 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
767 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
768 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
769 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
771 #define NLEGCOSAVER asize(leg_cosaver)
772 const char *leg_adip[] = {
777 #define NLEGADIP asize(leg_adip)
779 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
780 FILE *dip3d = NULL, *adip = NULL;
781 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
782 t_gkrbin *gkrbin = NULL;
783 gmx_enxnm_t *enm = NULL;
785 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
786 ener_file_t fmu = NULL;
787 int i, n, m, natom = 0, gnx_tot, teller, tel3;
789 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
791 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
794 gmx_bool bCorr, bTotal, bCont;
795 double M_diff = 0, epsilon, invtel, vol_aver;
796 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
797 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
798 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
801 rvec *slab_dipoles = NULL;
803 t_block *mols = NULL;
804 gmx_rmpbc_t gpbc = NULL;
816 /* initialize to a negative value so we can see that it wasn't picked up */
817 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
820 fmu = open_enx(mufn, "r");
821 do_enxnms(fmu, &nre, &enm);
823 /* Determine the indexes of the energy grps we need */
824 for (i = 0; (i < nre); i++)
826 if (strstr(enm[i].name, "Volume"))
830 else if (strstr(enm[i].name, "Mu-X"))
834 else if (strstr(enm[i].name, "Mu-Y"))
838 else if (strstr(enm[i].name, "Mu-Z"))
843 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
845 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
850 atom = top->atoms.atom;
853 if ((iVol == -1) && bMU)
855 printf("Using Volume from topology: %g nm^3\n", volume);
858 /* Correlation stuff */
859 bCorr = (corrtype[0] != 'n');
860 bTotal = (corrtype[0] == 't');
866 snew(muall[0], nframes*DIM);
871 for (i = 0; (i < gnx[0]); i++)
873 snew(muall[i], nframes*DIM);
878 /* Allocate array which contains for every molecule in a frame the
883 snew(dipole, gnx_tot);
888 for (i = 0; (i < DIM); i++)
890 Qlsq[i] = gmx_stats_init();
892 mulsq = gmx_stats_init();
894 /* Open all the files */
895 outmtot = xvgropen(out_mtot,
896 "Total dipole moment of the simulation box vs. time",
897 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
898 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
899 "Time (ps)", "", oenv);
900 outaver = xvgropen(out_aver, "Total dipole moment",
901 "Time (ps)", "D", oenv);
904 idim = axtitle[0] - 'X';
905 if ((idim < 0) || (idim >= DIM))
907 idim = axtitle[0] - 'x';
909 if ((idim < 0) || (idim >= DIM))
917 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
918 axtitle, nslices, idim);
921 snew(slab_dipoles, nslices);
922 fprintf(stderr, "Doing slab analysis\n");
928 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
929 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
934 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
935 "Average absolute dipole orientation", "Time (ps)", "", oenv);
936 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
942 snew(dipsp, gnx_tot);
944 /* we need a dummy file for gnuplot */
945 dip3d = (FILE *)gmx_ffopen("dummy.dat", "w");
946 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
949 dip3d = (FILE *)gmx_ffopen(fndip3d, "w");
952 gmx::BinaryInformationSettings settings;
953 settings.generatedByHeader(true);
954 settings.linePrefix("# ");
955 gmx::printBinaryInformation(dip3d, output_env_get_program_context(oenv),
958 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
961 /* Write legends to all the files */
962 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
963 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
965 if (bMU && (mu_aver == -1))
967 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
971 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
977 /* Read the first frame from energy or traj file */
982 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
985 timecheck = check_times(t);
990 if ((teller % 10) == 0)
992 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
997 printf("End of %s reached\n", mufn);
1001 while (bCont && (timecheck < 0));
1005 natom = read_first_x(oenv, &status, fn, &t, &x, box);
1008 /* Calculate spacing for dipole bin (simple histogram) */
1009 ndipbin = 1 + static_cast<int>(mu_max/0.01);
1010 snew(dipole_bin, ndipbin);
1012 for (m = 0; (m < DIM); m++)
1014 M[m] = M2[m] = M4[m] = 0.0;
1019 /* Use 0.7 iso 0.5 to account for pressure scaling */
1020 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1021 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1023 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1025 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1027 /* Start while loop over frames */
1032 if (bCorr && (teller >= nframes))
1037 srenew(muall[0], nframes*DIM);
1041 for (i = 0; (i < gnx_tot); i++)
1043 srenew(muall[i], nframes*DIM);
1049 muframelsq = gmx_stats_init();
1052 for (m = 0; (m < DIM); m++)
1059 /* Copy rvec into double precision local variable */
1060 for (m = 0; (m < DIM); m++)
1068 for (m = 0; (m < DIM); m++)
1073 gmx_rmpbc(gpbc, natom, box, x);
1075 /* Begin loop of all molecules in frame */
1076 for (n = 0; (n < ncos); n++)
1078 for (i = 0; (i < gnx[n]); i++)
1082 ind0 = mols->index[molindex[n][i]];
1083 ind1 = mols->index[molindex[n][i]+1];
1085 mol_dip(ind0, ind1, x, atom, dipole[i]);
1086 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1087 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1090 update_slab_dipoles(ind0, ind1, x,
1091 dipole[i], idim, nslices, slab_dipoles, box);
1095 mol_quad(ind0, ind1, x, atom, quad);
1096 for (m = 0; (m < DIM); m++)
1098 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1101 if (bCorr && !bTotal)
1104 muall[i][tel3+XX] = dipole[i][XX];
1105 muall[i][tel3+YY] = dipole[i][YY];
1106 muall[i][tel3+ZZ] = dipole[i][ZZ];
1109 for (m = 0; (m < DIM); m++)
1111 M_av[m] += dipole[i][m]; /* M per frame */
1112 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1114 mu_mol = sqrt(mu_mol);
1116 mu_ave += mu_mol; /* calc. the average mu */
1118 /* Update the dipole distribution */
1119 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1127 rvec2sprvec(dipole[i], dipsp[i]);
1129 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1131 if (dipsp[i][ZZ] < 0.5 * M_PI)
1140 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1142 if (dipsp[i][ZZ] < 0.5 * M_PI)
1151 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1153 if (dipsp[i][ZZ] < 0.5 * M_PI)
1162 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1164 if (dipsp[i][ZZ] < 0.5 * M_PI)
1175 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1180 x[ind0][XX]+dipole[i][XX]/25,
1181 x[ind0][YY]+dipole[i][YY]/25,
1182 x[ind0][ZZ]+dipole[i][ZZ]/25,
1186 } /* End loop of all molecules in frame */
1190 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1191 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1192 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1193 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1194 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1195 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1199 /* Compute square of total dipole */
1200 for (m = 0; (m < DIM); m++)
1202 M_av2[m] = M_av[m]*M_av[m];
1207 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1208 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1209 sqr(dipaxis[YY]-0.5)+
1210 sqr(dipaxis[ZZ]-0.5));
1213 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1214 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1218 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1219 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1225 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1232 muall[0][tel3+XX] = M_av[XX];
1233 muall[0][tel3+YY] = M_av[YY];
1234 muall[0][tel3+ZZ] = M_av[ZZ];
1237 /* Write to file the total dipole moment of the box, and its components
1240 if ((skip == 0) || ((teller % skip) == 0))
1242 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1243 t, M_av[XX], M_av[YY], M_av[ZZ],
1244 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1247 for (m = 0; (m < DIM); m++)
1251 M4[m] += sqr(M_av2[m]);
1253 /* Increment loop counter */
1256 /* Calculate for output the running averages */
1257 invtel = 1.0/teller;
1258 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1259 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1260 M_diff = M2_ave - M_ave2;
1262 /* Compute volume from box in traj, else we use the one from above */
1269 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1271 /* Calculate running average for dipole */
1274 mu_aver = (mu_ave/gnx_tot)*invtel;
1277 if ((skip == 0) || ((teller % skip) == 0))
1279 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1280 * the two. Here M is sum mu_i. Further write the finite system
1281 * Kirkwood G factor and epsilon.
1283 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1284 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1289 gmx_stats_get_average(muframelsq, &aver);
1290 fprintf(adip, "%10g %f \n", t, aver);
1293 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1295 if (!bMU || (mu_aver != -1))
1297 /* Finite system Kirkwood G-factor */
1298 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1299 /* Infinite system Kirkwood G-factor */
1300 if (epsilonRF == 0.0)
1302 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1306 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1307 Gk/(3*epsilon*(2*epsilonRF+1)));
1310 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1315 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1318 gmx_stats_done(muframelsq);
1322 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1326 bCont = read_next_x(oenv, status, &t, x, box);
1328 timecheck = check_times(t);
1330 while (bCont && (timecheck == 0) );
1332 gmx_rmpbc_done(gpbc);
1339 gmx_ffclose(outmtot);
1340 gmx_ffclose(outaver);
1341 gmx_ffclose(outeps);
1355 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1356 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1357 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1358 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1359 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1365 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1366 sfree(slab_dipoles);
1370 printf("Average volume over run is %g\n", vol_aver);
1373 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1374 print_cmap(cmap, gkrbin, nlevels);
1376 /* Autocorrelation function */
1381 printf("Not enough frames for autocorrelation\n");
1385 dt = (t1 - t0)/(teller-1);
1386 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1392 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1393 teller, 1, muall, dt, mode, TRUE);
1397 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1398 teller, gnx_tot, muall, dt,
1399 mode, strcmp(corrtype, "molsep"));
1405 real aver, sigma, error;
1407 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1408 printf("\nDipole moment (Debye)\n");
1409 printf("---------------------\n");
1410 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1411 aver, sigma, error);
1415 for (m = 0; (m < DIM); m++)
1417 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1420 printf("\nQuadrupole moment (Debye-Ang)\n");
1421 printf("-----------------------------\n");
1422 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1423 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1424 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1428 printf("The following averages for the complete trajectory have been calculated:\n\n");
1429 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1430 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1431 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1433 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1434 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1435 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1437 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1438 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1440 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1442 if (!bMU || (mu_aver != -1))
1444 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1445 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1447 printf("Epsilon = %g\n", epsilon);
1451 /* Write to file the dipole moment distibution during the simulation.
1453 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1454 for (i = 0; (i < ndipbin); i++)
1456 fprintf(outdd, "%10g %10f\n",
1457 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1464 done_gkrbin(&gkrbin);
1468 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1477 while (m < mols->nr && index[i] != mols->index[m])
1483 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1485 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1487 if (i >= *n || index[i] != j)
1489 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1493 /* Modify the index in place */
1496 printf("There are %d molecules in the selection\n", nmol);
1500 int gmx_dipoles(int argc, char *argv[])
1502 const char *desc[] = {
1503 "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1504 "system. From this you can compute e.g. the dielectric constant for",
1505 "low-dielectric media.",
1506 "For molecules with a net charge, the net charge is subtracted at",
1507 "center of mass of the molecule.[PAR]",
1508 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1509 "components as well as the norm of the vector.",
1510 "The file [TT]aver.xvg[tt] contains [CHEVRON][MAG][GRK]mu[grk][mag]^2[chevron] and [MAG][CHEVRON][GRK]mu[grk][chevron][mag]^2 during the",
1512 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1514 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1515 "Furthermore, the dipole autocorrelation function will be computed when",
1516 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1518 "The correlation functions can be averaged over all molecules",
1519 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1520 "or it can be computed over the total dipole moment of the simulation box",
1521 "([TT]total[tt]).[PAR]",
1522 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1523 "G-factor, as well as the average cosine of the angle between the dipoles",
1524 "as a function of the distance. The plot also includes gOO and hOO",
1525 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1526 "we also include the energy per scale computed by taking the inner product of",
1527 "the dipoles divided by the distance to the third power.[PAR]",
1530 "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1531 "This will calculate the autocorrelation function of the molecular",
1532 "dipoles using a first order Legendre polynomial of the angle of the",
1533 "dipole vector and itself a time t later. For this calculation 1001",
1534 "frames will be used. Further, the dielectric constant will be calculated",
1535 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1536 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1537 "distribution function a maximum of 5.0 will be used."
1539 real mu_max = 5, mu_aver = -1, rcmax = 0;
1540 real epsilonRF = 0.0, temp = 300;
1541 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1542 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1543 const char *axtitle = "Z";
1544 int nslices = 10; /* nr of slices defined */
1545 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1546 int nlevels = 20, ndegrees = 90;
1549 { "-mu", FALSE, etREAL, {&mu_aver},
1550 "dipole of a single molecule (in Debye)" },
1551 { "-mumax", FALSE, etREAL, {&mu_max},
1552 "max dipole in Debye (for histogram)" },
1553 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1554 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1555 { "-skip", FALSE, etINT, {&skip},
1556 "Skip steps in the output (but not in the computations)" },
1557 { "-temp", FALSE, etREAL, {&temp},
1558 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1559 { "-corr", FALSE, etENUM, {corrtype},
1560 "Correlation function to calculate" },
1561 { "-pairs", FALSE, etBOOL, {&bPairs},
1562 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1563 { "-quad", FALSE, etBOOL, {&bQuad},
1564 "Take quadrupole into account"},
1565 { "-ncos", FALSE, etINT, {&ncos},
1566 "Must be 1 or 2. Determines whether the [CHEVRON][COS][GRK]theta[grk][cos][chevron] is computed between all molecules in one group, or between molecules in two different groups. This turns on the [TT]-g[tt] flag." },
1567 { "-axis", FALSE, etSTR, {&axtitle},
1568 "Take the normal on the computational box in direction X, Y or Z." },
1569 { "-sl", FALSE, etINT, {&nslices},
1570 "Divide the box into this number of slices." },
1571 { "-gkratom", FALSE, etINT, {&nFA},
1572 "Use the n-th atom of a molecule (starting from 1) to calculate the distance between molecules rather than the center of charge (when 0) in the calculation of distance dependent Kirkwood factors" },
1573 { "-gkratom2", FALSE, etINT, {&nFB},
1574 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1575 { "-rcmax", FALSE, etREAL, {&rcmax},
1576 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1577 { "-phi", FALSE, etBOOL, {&bPhi},
1578 "Plot the 'torsion angle' defined as the rotation of the two dipole vectors around the distance vector between the two molecules in the [TT].xpm[tt] file from the [TT]-cmap[tt] option. By default the cosine of the angle between the dipoles is plotted." },
1579 { "-nlevels", FALSE, etINT, {&nlevels},
1580 "Number of colors in the cmap output" },
1581 { "-ndegrees", FALSE, etINT, {&ndegrees},
1582 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1587 char **grpname = NULL;
1588 gmx_bool bGkr, bMU, bSlab;
1590 { efEDR, "-en", NULL, ffOPTRD },
1591 { efTRX, "-f", NULL, ffREAD },
1592 { efTPX, NULL, NULL, ffREAD },
1593 { efNDX, NULL, NULL, ffOPTRD },
1594 { efXVG, "-o", "Mtot", ffWRITE },
1595 { efXVG, "-eps", "epsilon", ffWRITE },
1596 { efXVG, "-a", "aver", ffWRITE },
1597 { efXVG, "-d", "dipdist", ffWRITE },
1598 { efXVG, "-c", "dipcorr", ffOPTWR },
1599 { efXVG, "-g", "gkr", ffOPTWR },
1600 { efXVG, "-adip", "adip", ffOPTWR },
1601 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1602 { efXVG, "-cos", "cosaver", ffOPTWR },
1603 { efXPM, "-cmap", "cmap", ffOPTWR },
1604 { efXVG, "-slab", "slab", ffOPTWR }
1606 #define NFILE asize(fnm)
1615 ppa = add_acf_pargs(&npargs, pa);
1616 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1617 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1622 printf("Using %g as mu_max and %g as the dipole moment.\n",
1624 if (epsilonRF == 0.0)
1626 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1629 bMU = opt2bSet("-en", NFILE, fnm);
1632 gmx_fatal(FARGS, "Due to new ways of treating molecules in GROMACS the total dipole in the energy file may be incorrect, because molecules can be split over periodic boundary conditions before computing the dipole. Please use your trajectory file.");
1634 bGkr = opt2bSet("-g", NFILE, fnm);
1635 if (opt2parg_bSet("-ncos", asize(pa), pa))
1637 if ((ncos != 1) && (ncos != 2))
1639 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1643 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1644 opt2parg_bSet("-axis", asize(pa), pa));
1649 printf("WARNING: Can not determine quadrupoles from energy file\n");
1654 printf("WARNING: Can not determine Gk(r) from energy file\n");
1660 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1661 " not enter a valid dipole for the molecules\n");
1666 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1667 &natoms, NULL, NULL, NULL, top);
1670 snew(grpname, ncos);
1671 snew(grpindex, ncos);
1672 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1673 ncos, gnx, grpindex, grpname);
1674 for (k = 0; (k < ncos); k++)
1676 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1677 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1681 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1682 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1683 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1684 opt2fn_null("-cos", NFILE, fnm),
1685 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1686 bPairs, corrtype[0],
1687 opt2fn("-c", NFILE, fnm),
1688 bGkr, opt2fn("-g", NFILE, fnm),
1689 bPhi, &nlevels, ndegrees,
1691 opt2fn("-cmap", NFILE, fnm), rcmax,
1692 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1693 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1694 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1696 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1697 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1698 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1699 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1700 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");