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, 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.
46 #include "gromacs/commandline/pargs.h"
52 #include "gromacs/fileio/futil.h"
55 #include "gmx_statistics.h"
62 #include "gromacs/fileio/enxio.h"
64 #include "gromacs/fileio/matio.h"
67 #include "gromacs/fileio/trxio.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/programinfo.h"
72 #define e2d(x) ENM2DEBYE*(x)
73 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
74 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
86 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
94 if ((ptr = getenv("GKRWIDTH")) != NULL)
96 double bw = strtod(ptr, NULL);
101 gb->spacing = 0.01; /* nm */
103 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
110 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
113 snew(gb->elem, gb->nelem);
114 snew(gb->count, gb->nelem);
116 snew(gb->cmap, gb->nx);
117 gb->ny = std::max(2, ndegrees);
118 for (i = 0; (i < gb->nx); i++)
120 snew(gb->cmap[i], gb->ny);
127 static void done_gkrbin(t_gkrbin **gb)
135 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
137 int cy, index = gmx_nint(r/gb->spacing);
140 if (index < gb->nelem)
142 gb->elem[index] += cosa;
150 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
154 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
156 cy = std::min(gb->ny-1, std::max(0, cy));
159 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
161 gb->cmap[index][cy] += 1;
165 static void rvec2sprvec(rvec dipcart, rvec dipsp)
168 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
169 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
170 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
175 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
176 int mindex[], rvec x[], rvec mu[],
177 int ePBC, matrix box, t_atom *atom, int *nAtom)
179 static rvec *xcm[2] = { NULL, NULL};
180 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
181 real qtot, q, cosa, rr, phi;
185 for (n = 0; (n < ncos); n++)
189 snew(xcm[n], ngrp[n]);
191 for (i = 0; (i < ngrp[n]); i++)
193 /* Calculate center of mass of molecule */
199 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
204 clear_rvec(xcm[n][i]);
206 for (j = j0; j < j1; j++)
210 for (k = 0; k < DIM; k++)
212 xcm[n][i][k] += q*x[j][k];
215 svmul(1/qtot, xcm[n][i], xcm[n][i]);
219 set_pbc(&pbc, ePBC, box);
222 for (i = 0; i < ngrp[grp0]; i++)
224 gi = molindex[grp0][i];
225 j0 = (ncos == 2) ? 0 : i+1;
226 for (j = j0; j < ngrp[grp1]; j++)
228 gj = molindex[grp1][j];
229 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
231 /* Compute distance between molecules including PBC */
232 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
238 rvec r_ij, r_kj, r_kl, mm, nn;
242 copy_rvec(xcm[grp0][i], xj);
243 copy_rvec(xcm[grp1][j], xk);
244 rvec_add(xj, mu[gi], xi);
245 rvec_add(xk, mu[gj], xl);
246 phi = dih_angle(xi, xj, xk, xl, &pbc,
247 r_ij, r_kj, r_kl, mm, nn, /* out */
248 &sign, &t1, &t2, &t3);
253 cosa = cos_angle(mu[gi], mu[gj]);
256 if (debug || gmx_isnan(cosa))
258 fprintf(debug ? debug : stderr,
259 "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",
260 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
261 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
265 add2gkr(gb, rr, cosa, phi);
271 static real normalize_cmap(t_gkrbin *gb)
278 for (i = 0; (i < gb->nx); i++)
280 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
281 for (j = 0; (j < gb->ny); j++)
283 gb->cmap[i][j] /= vol;
284 hi = std::max(hi, gb->cmap[i][j]);
289 gmx_fatal(FARGS, "No data in the cmap");
294 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
301 t_rgb rlo = { 1, 1, 1 };
302 t_rgb rhi = { 0, 0, 0 };
304 hi = normalize_cmap(gb);
305 snew(xaxis, gb->nx+1);
306 for (i = 0; (i < gb->nx+1); i++)
308 xaxis[i] = i*gb->spacing;
311 for (j = 0; (j < gb->ny); j++)
315 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
319 yaxis[j] = (180.0*j)/(gb->ny-1.0);
321 /*2.0*j/(gb->ny-1.0)-1.0;*/
323 out = ffopen(cmap, "w");
325 "Dipole Orientation Distribution", "Fraction", "r (nm)",
326 gb->bPhi ? "Phi" : "Alpha",
327 gb->nx, gb->ny, xaxis, yaxis,
328 gb->cmap, 0, hi, rlo, rhi, nlevels);
334 static void print_gkrbin(const char *fn, t_gkrbin *gb,
335 int ngrp, int nframes, real volume,
336 const output_env_t oenv)
338 /* We compute Gk(r), gOO and hOO according to
339 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
340 * In this implementation the angle between dipoles is stored
341 * rather than their inner product. This allows to take polarizible
342 * models into account. The RDF is calculated as well, almost for free!
345 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
347 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
350 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
351 xvgr_legend(fp, asize(leg), leg, oenv);
353 Gkr = 1; /* Self-dipole inproduct = 1 */
358 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
359 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
363 while (last > 1 && gb->elem[last-1] == 0)
368 /* Divide by dipole squared, by number of frames, by number of origins.
369 * Multiply by 2 because we only take half the matrix of interactions
372 fac = 2.0/((double) ngrp * (double) nframes);
375 for (i = 0; i < last; i++)
377 /* Centre of the coordinate in the spherical layer */
380 /* Volume of the layer */
381 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
384 gOO = gb->count[i]*fac/(rho*vol_s);
386 /* Dipole correlation hOO, normalized by the relative number density, like
387 * in a Radial distribution function.
389 ggg = gb->elem[i]*fac;
390 hOO = 3.0*ggg/(rho*vol_s);
394 cosav = gb->elem[i]/gb->count[i];
400 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
402 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
403 x1, Gkr, cosav, hOO, gOO, ener);
411 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
412 real *t, int nre, t_enxframe *fr)
418 bCont = do_enx(fmu, fr);
421 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
422 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
427 if (Vol != -1) /* we've got Volume in the energy file */
429 *vol = fr->ener[Vol].e;
431 for (i = 0; i < DIM; i++)
433 mu[i] = fr->ener[iMu[i]].e;
441 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
444 int ncharged, m, a0, a1, a;
447 for (m = 0; m < n; m++)
449 a0 = mols->index[index[m]];
450 a1 = mols->index[index[m]+1];
453 for (a = a0; a < a1; a++)
458 /* This check is only for the count print */
459 if (fabs(qtot) > 0.01)
465 /* Remove the net charge at the center of mass */
466 for (a = a0; a < a1; a++)
468 atom[a].q -= qtot*atom[a].m/mtot;
475 printf("There are %d charged molecules in the selection,\n"
476 "will subtract their charge at their center of mass\n", ncharged);
480 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
486 for (k = k0; (k < k1); k++)
489 for (m = 0; (m < DIM); m++)
496 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
498 int i, k, m, n, niter;
499 real q, r2, mass, masstot;
500 rvec com; /* center of mass */
501 rvec r; /* distance of atoms to center of mass */
503 double dd[DIM], **ev, tmp;
507 for (i = 0; (i < DIM); i++)
514 /* Compute center of mass */
517 for (k = k0; (k < k1); k++)
521 for (i = 0; (i < DIM); i++)
523 com[i] += mass*x[k][i];
526 svmul((1.0/masstot), com, com);
528 /* We want traceless quadrupole moments, so let us calculate the complete
529 * quadrupole moment tensor and diagonalize this tensor to get
530 * the individual components on the diagonal.
533 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
535 for (m = 0; (m < DIM); m++)
537 for (n = 0; (n < DIM); n++)
542 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
544 q = (atom[k].q)*100.0;
545 rvec_sub(x[k], com, r);
547 for (m = 0; (m < DIM); m++)
549 for (n = 0; (n < DIM); n++)
551 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
557 for (i = 0; (i < DIM); i++)
559 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
560 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
564 /* We've got the quadrupole tensor, now diagonalize the sucker */
565 jacobi(inten, 3, dd, ev, &niter);
569 for (i = 0; (i < DIM); i++)
571 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
572 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
574 for (i = 0; (i < DIM); i++)
576 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
577 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
580 /* Sort the eigenvalues, for water we know that the order is as follows:
584 * At the moment I have no idea how this will work out for other molecules...
588 if (dd[i+1] > dd[i]) { \
597 for (m = 0; (m < DIM); m++)
599 quad[0] = dd[2]; /* yy */
600 quad[1] = dd[0]; /* zz */
601 quad[2] = dd[1]; /* xx */
606 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
610 for (i = 0; (i < DIM); i++)
620 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
622 real calc_eps(double M_diff, double volume, double epsRF, double temp)
624 double eps, A, teller, noemer;
625 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
626 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
628 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
637 teller = 1 + (A*2*epsRF/(2*epsRF+1));
638 noemer = 1 - (A/(2*epsRF+1));
640 eps = teller / noemer;
645 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
646 int idim, int nslice, rvec slab_dipole[],
652 for (k = k0; (k < k1); k++)
657 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
658 rvec_inc(slab_dipole[k], mu);
661 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
662 rvec slab_dipole[], matrix box, int nframes,
663 const output_env_t oenv)
669 const char *leg_dim[4] = {
670 "\\f{12}m\\f{4}\\sX\\N",
671 "\\f{12}m\\f{4}\\sY\\N",
672 "\\f{12}m\\f{4}\\sZ\\N",
673 "\\f{12}m\\f{4}\\stot\\N"
676 sprintf(buf, "Box-%c (nm)", 'X'+idim);
677 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
679 xvgr_legend(fp, DIM, leg_dim, oenv);
680 for (i = 0; (i < nslice); i++)
682 mutot = norm(slab_dipole[i])/nframes;
683 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
684 ((i+0.5)*box[idim][idim])/nslice,
685 slab_dipole[i][XX]/nframes,
686 slab_dipole[i][YY]/nframes,
687 slab_dipole[i][ZZ]/nframes,
691 do_view(oenv, fn, "-autoscale xy -nxy");
694 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
697 double dc, d, ddc1, ddc2, ddc3;
698 rvec xxx = { 1, 0, 0 };
699 rvec yyy = { 0, 1, 0 };
700 rvec zzz = { 0, 0, 1 };
703 ddc1 = ddc2 = ddc3 = 0;
704 for (i = k = 0; (i < n); i++)
706 ddc1 += fabs(cos_angle(dip[i], xxx));
707 ddc2 += fabs(cos_angle(dip[i], yyy));
708 ddc3 += fabs(cos_angle(dip[i], zzz));
711 for (j = i+1; (j < n); j++, k++)
713 dc = cos_angle(dip[i], dip[j]);
724 static void do_dip(t_topology *top, int ePBC, real volume,
726 const char *out_mtot, const char *out_eps,
727 const char *out_aver, const char *dipdist,
728 const char *cosaver, const char *fndip3d,
729 const char *fnadip, gmx_bool bPairs,
730 const char *corrtype, const char *corf,
731 gmx_bool bGkr, const char *gkrfn,
732 gmx_bool bPhi, int *nlevels, int ndegrees,
734 const char *cmap, real rcmax,
736 gmx_bool bMU, const char *mufn,
737 int *gnx, int *molindex[],
738 real mu_max, real mu_aver,
739 real epsilonRF, real temp,
740 int *gkatom, int skip,
741 gmx_bool bSlab, int nslices,
742 const char *axtitle, const char *slabfn,
743 const output_env_t oenv)
745 const char *leg_mtot[] = {
751 #define NLEGMTOT asize(leg_mtot)
752 const char *leg_eps[] = {
757 #define NLEGEPS asize(leg_eps)
758 const char *leg_aver[] = {
761 "< |M|\\S2\\N > - < |M| >\\S2\\N",
762 "< |M| >\\S2\\N / < |M|\\S2\\N >"
764 #define NLEGAVER asize(leg_aver)
765 const char *leg_cosaver[] = {
766 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
768 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
769 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
770 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
772 #define NLEGCOSAVER asize(leg_cosaver)
773 const char *leg_adip[] = {
778 #define NLEGADIP asize(leg_adip)
780 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
781 FILE *dip3d = NULL, *adip = NULL;
782 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
783 t_gkrbin *gkrbin = NULL;
784 gmx_enxnm_t *enm = NULL;
786 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
787 ener_file_t fmu = NULL;
788 int i, n, m, natom = 0, gnx_tot, teller, tel3;
790 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
792 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
795 gmx_bool bCorr, bTotal, bCont;
796 double M_diff = 0, epsilon, invtel, vol_aver;
797 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
798 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
799 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
802 rvec *slab_dipoles = NULL;
804 t_block *mols = NULL;
805 gmx_rmpbc_t gpbc = NULL;
817 /* initialize to a negative value so we can see that it wasn't picked up */
818 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
821 fmu = open_enx(mufn, "r");
822 do_enxnms(fmu, &nre, &enm);
824 /* Determine the indexes of the energy grps we need */
825 for (i = 0; (i < nre); i++)
827 if (strstr(enm[i].name, "Volume"))
831 else if (strstr(enm[i].name, "Mu-X"))
835 else if (strstr(enm[i].name, "Mu-Y"))
839 else if (strstr(enm[i].name, "Mu-Z"))
844 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
846 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
851 atom = top->atoms.atom;
854 if ((iVol == -1) && bMU)
856 printf("Using Volume from topology: %g nm^3\n", volume);
859 /* Correlation stuff */
860 bCorr = (corrtype[0] != 'n');
861 bTotal = (corrtype[0] == 't');
867 snew(muall[0], nframes*DIM);
872 for (i = 0; (i < gnx[0]); i++)
874 snew(muall[i], nframes*DIM);
879 /* Allocate array which contains for every molecule in a frame the
884 snew(dipole, gnx_tot);
889 for (i = 0; (i < DIM); i++)
891 Qlsq[i] = gmx_stats_init();
893 mulsq = gmx_stats_init();
895 /* Open all the files */
896 outmtot = xvgropen(out_mtot,
897 "Total dipole moment of the simulation box vs. time",
898 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
899 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
900 "Time (ps)", "", oenv);
901 outaver = xvgropen(out_aver, "Total dipole moment",
902 "Time (ps)", "D", oenv);
905 idim = axtitle[0] - 'X';
906 if ((idim < 0) || (idim >= DIM))
908 idim = axtitle[0] - 'x';
910 if ((idim < 0) || (idim >= DIM))
918 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
919 axtitle, nslices, idim);
922 snew(slab_dipoles, nslices);
923 fprintf(stderr, "Doing slab analysis\n");
929 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
930 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
935 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
936 "Average absolute dipole orientation", "Time (ps)", "", oenv);
937 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
943 snew(dipsp, gnx_tot);
945 /* we need a dummy file for gnuplot */
946 dip3d = (FILE *)ffopen("dummy.dat", "w");
947 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
950 dip3d = (FILE *)ffopen(fndip3d, "w");
953 gmx::BinaryInformationSettings settings;
954 settings.generatedByHeader(true);
955 settings.linePrefix("# ");
956 gmx::printBinaryInformation(dip3d, gmx::ProgramInfo::getInstance(),
959 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
962 /* Write legends to all the files */
963 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
964 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
966 if (bMU && (mu_aver == -1))
968 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
972 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
978 /* Read the first frame from energy or traj file */
983 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
986 timecheck = check_times(t);
991 if ((teller % 10) == 0)
993 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
998 printf("End of %s reached\n", mufn);
1002 while (bCont && (timecheck < 0));
1006 natom = read_first_x(oenv, &status, fn, &t, &x, box);
1009 /* Calculate spacing for dipole bin (simple histogram) */
1010 ndipbin = 1 + static_cast<int>(mu_max/0.01);
1011 snew(dipole_bin, ndipbin);
1013 for (m = 0; (m < DIM); m++)
1015 M[m] = M2[m] = M4[m] = 0.0;
1020 /* Use 0.7 iso 0.5 to account for pressure scaling */
1021 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1022 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1024 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1026 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1028 /* Start while loop over frames */
1033 if (bCorr && (teller >= nframes))
1038 srenew(muall[0], nframes*DIM);
1042 for (i = 0; (i < gnx_tot); i++)
1044 srenew(muall[i], nframes*DIM);
1050 muframelsq = gmx_stats_init();
1053 for (m = 0; (m < DIM); m++)
1060 /* Copy rvec into double precision local variable */
1061 for (m = 0; (m < DIM); m++)
1069 for (m = 0; (m < DIM); m++)
1074 gmx_rmpbc(gpbc, natom, box, x);
1076 /* Begin loop of all molecules in frame */
1077 for (n = 0; (n < ncos); n++)
1079 for (i = 0; (i < gnx[n]); i++)
1083 ind0 = mols->index[molindex[n][i]];
1084 ind1 = mols->index[molindex[n][i]+1];
1086 mol_dip(ind0, ind1, x, atom, dipole[i]);
1087 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1088 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1091 update_slab_dipoles(ind0, ind1, x,
1092 dipole[i], idim, nslices, slab_dipoles, box);
1096 mol_quad(ind0, ind1, x, atom, quad);
1097 for (m = 0; (m < DIM); m++)
1099 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1102 if (bCorr && !bTotal)
1105 muall[i][tel3+XX] = dipole[i][XX];
1106 muall[i][tel3+YY] = dipole[i][YY];
1107 muall[i][tel3+ZZ] = dipole[i][ZZ];
1110 for (m = 0; (m < DIM); m++)
1112 M_av[m] += dipole[i][m]; /* M per frame */
1113 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1115 mu_mol = sqrt(mu_mol);
1117 mu_ave += mu_mol; /* calc. the average mu */
1119 /* Update the dipole distribution */
1120 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1128 rvec2sprvec(dipole[i], dipsp[i]);
1130 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1132 if (dipsp[i][ZZ] < 0.5 * M_PI)
1141 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1143 if (dipsp[i][ZZ] < 0.5 * M_PI)
1152 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1154 if (dipsp[i][ZZ] < 0.5 * M_PI)
1163 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1165 if (dipsp[i][ZZ] < 0.5 * M_PI)
1176 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1181 x[ind0][XX]+dipole[i][XX]/25,
1182 x[ind0][YY]+dipole[i][YY]/25,
1183 x[ind0][ZZ]+dipole[i][ZZ]/25,
1187 } /* End loop of all molecules in frame */
1191 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1192 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1193 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1194 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1195 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1196 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1200 /* Compute square of total dipole */
1201 for (m = 0; (m < DIM); m++)
1203 M_av2[m] = M_av[m]*M_av[m];
1208 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1209 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1210 sqr(dipaxis[YY]-0.5)+
1211 sqr(dipaxis[ZZ]-0.5));
1214 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1215 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1219 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1220 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1226 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1233 muall[0][tel3+XX] = M_av[XX];
1234 muall[0][tel3+YY] = M_av[YY];
1235 muall[0][tel3+ZZ] = M_av[ZZ];
1238 /* Write to file the total dipole moment of the box, and its components
1241 if ((skip == 0) || ((teller % skip) == 0))
1243 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1244 t, M_av[XX], M_av[YY], M_av[ZZ],
1245 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1248 for (m = 0; (m < DIM); m++)
1252 M4[m] += sqr(M_av2[m]);
1254 /* Increment loop counter */
1257 /* Calculate for output the running averages */
1258 invtel = 1.0/teller;
1259 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1260 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1261 M_diff = M2_ave - M_ave2;
1263 /* Compute volume from box in traj, else we use the one from above */
1270 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1272 /* Calculate running average for dipole */
1275 mu_aver = (mu_ave/gnx_tot)*invtel;
1278 if ((skip == 0) || ((teller % skip) == 0))
1280 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1281 * the two. Here M is sum mu_i. Further write the finite system
1282 * Kirkwood G factor and epsilon.
1284 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1285 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1290 gmx_stats_get_average(muframelsq, &aver);
1291 fprintf(adip, "%10g %f \n", t, aver);
1294 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1296 if (!bMU || (mu_aver != -1))
1298 /* Finite system Kirkwood G-factor */
1299 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1300 /* Infinite system Kirkwood G-factor */
1301 if (epsilonRF == 0.0)
1303 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1307 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1308 Gk/(3*epsilon*(2*epsilonRF+1)));
1311 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1316 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1319 gmx_stats_done(muframelsq);
1323 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1327 bCont = read_next_x(oenv, status, &t, x, box);
1329 timecheck = check_times(t);
1331 while (bCont && (timecheck == 0) );
1333 gmx_rmpbc_done(gpbc);
1356 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1357 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1358 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1359 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1360 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1366 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1367 sfree(slab_dipoles);
1371 printf("Average volume over run is %g\n", vol_aver);
1374 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1375 print_cmap(cmap, gkrbin, nlevels);
1377 /* Autocorrelation function */
1382 printf("Not enough frames for autocorrelation\n");
1386 dt = (t1 - t0)/(teller-1);
1387 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1393 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1394 teller, 1, muall, dt, mode, TRUE);
1398 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1399 teller, gnx_tot, muall, dt,
1400 mode, strcmp(corrtype, "molsep"));
1406 real aver, sigma, error;
1408 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1409 printf("\nDipole moment (Debye)\n");
1410 printf("---------------------\n");
1411 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1412 aver, sigma, error);
1416 for (m = 0; (m < DIM); m++)
1418 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1421 printf("\nQuadrupole moment (Debye-Ang)\n");
1422 printf("-----------------------------\n");
1423 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1424 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1425 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1429 printf("The following averages for the complete trajectory have been calculated:\n\n");
1430 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1431 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1432 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1434 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1435 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1436 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1438 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1439 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1441 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1443 if (!bMU || (mu_aver != -1))
1445 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1446 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1448 printf("Epsilon = %g\n", epsilon);
1452 /* Write to file the dipole moment distibution during the simulation.
1454 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1455 for (i = 0; (i < ndipbin); i++)
1457 fprintf(outdd, "%10g %10f\n",
1458 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1465 done_gkrbin(&gkrbin);
1469 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1478 while (m < mols->nr && index[i] != mols->index[m])
1484 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1486 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1488 if (i >= *n || index[i] != j)
1490 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1494 /* Modify the index in place */
1497 printf("There are %d molecules in the selection\n", nmol);
1501 int gmx_dipoles(int argc, char *argv[])
1503 const char *desc[] = {
1504 "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1505 "system. From this you can compute e.g. the dielectric constant for",
1506 "low-dielectric media.",
1507 "For molecules with a net charge, the net charge is subtracted at",
1508 "center of mass of the molecule.[PAR]",
1509 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1510 "components as well as the norm of the vector.",
1511 "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",
1513 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1515 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1516 "Furthermore, the dipole autocorrelation function will be computed when",
1517 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1519 "The correlation functions can be averaged over all molecules",
1520 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1521 "or it can be computed over the total dipole moment of the simulation box",
1522 "([TT]total[tt]).[PAR]",
1523 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1524 "G-factor, as well as the average cosine of the angle between the dipoles",
1525 "as a function of the distance. The plot also includes gOO and hOO",
1526 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1527 "we also include the energy per scale computed by taking the inner product of",
1528 "the dipoles divided by the distance to the third power.[PAR]",
1531 "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1532 "This will calculate the autocorrelation function of the molecular",
1533 "dipoles using a first order Legendre polynomial of the angle of the",
1534 "dipole vector and itself a time t later. For this calculation 1001",
1535 "frames will be used. Further, the dielectric constant will be calculated",
1536 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1537 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1538 "distribution function a maximum of 5.0 will be used."
1540 real mu_max = 5, mu_aver = -1, rcmax = 0;
1541 real epsilonRF = 0.0, temp = 300;
1542 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1543 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1544 const char *axtitle = "Z";
1545 int nslices = 10; /* nr of slices defined */
1546 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1547 int nlevels = 20, ndegrees = 90;
1550 { "-mu", FALSE, etREAL, {&mu_aver},
1551 "dipole of a single molecule (in Debye)" },
1552 { "-mumax", FALSE, etREAL, {&mu_max},
1553 "max dipole in Debye (for histogram)" },
1554 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1555 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1556 { "-skip", FALSE, etINT, {&skip},
1557 "Skip steps in the output (but not in the computations)" },
1558 { "-temp", FALSE, etREAL, {&temp},
1559 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1560 { "-corr", FALSE, etENUM, {corrtype},
1561 "Correlation function to calculate" },
1562 { "-pairs", FALSE, etBOOL, {&bPairs},
1563 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1564 { "-quad", FALSE, etBOOL, {&bQuad},
1565 "Take quadrupole into account"},
1566 { "-ncos", FALSE, etINT, {&ncos},
1567 "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." },
1568 { "-axis", FALSE, etSTR, {&axtitle},
1569 "Take the normal on the computational box in direction X, Y or Z." },
1570 { "-sl", FALSE, etINT, {&nslices},
1571 "Divide the box into this number of slices." },
1572 { "-gkratom", FALSE, etINT, {&nFA},
1573 "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" },
1574 { "-gkratom2", FALSE, etINT, {&nFB},
1575 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1576 { "-rcmax", FALSE, etREAL, {&rcmax},
1577 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1578 { "-phi", FALSE, etBOOL, {&bPhi},
1579 "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." },
1580 { "-nlevels", FALSE, etINT, {&nlevels},
1581 "Number of colors in the cmap output" },
1582 { "-ndegrees", FALSE, etINT, {&ndegrees},
1583 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1588 char **grpname = NULL;
1589 gmx_bool bGkr, bMU, bSlab;
1591 { efEDR, "-en", NULL, ffOPTRD },
1592 { efTRX, "-f", NULL, ffREAD },
1593 { efTPX, NULL, NULL, ffREAD },
1594 { efNDX, NULL, NULL, ffOPTRD },
1595 { efXVG, "-o", "Mtot", ffWRITE },
1596 { efXVG, "-eps", "epsilon", ffWRITE },
1597 { efXVG, "-a", "aver", ffWRITE },
1598 { efXVG, "-d", "dipdist", ffWRITE },
1599 { efXVG, "-c", "dipcorr", ffOPTWR },
1600 { efXVG, "-g", "gkr", ffOPTWR },
1601 { efXVG, "-adip", "adip", ffOPTWR },
1602 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1603 { efXVG, "-cos", "cosaver", ffOPTWR },
1604 { efXPM, "-cmap", "cmap", ffOPTWR },
1605 { efXVG, "-slab", "slab", ffOPTWR }
1607 #define NFILE asize(fnm)
1616 ppa = add_acf_pargs(&npargs, pa);
1617 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1618 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1623 printf("Using %g as mu_max and %g as the dipole moment.\n",
1625 if (epsilonRF == 0.0)
1627 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1630 bMU = opt2bSet("-en", NFILE, fnm);
1633 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.");
1635 bGkr = opt2bSet("-g", NFILE, fnm);
1636 if (opt2parg_bSet("-ncos", asize(pa), pa))
1638 if ((ncos != 1) && (ncos != 2))
1640 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1644 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1645 opt2parg_bSet("-axis", asize(pa), pa));
1650 printf("WARNING: Can not determine quadrupoles from energy file\n");
1655 printf("WARNING: Can not determine Gk(r) from energy file\n");
1661 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1662 " not enter a valid dipole for the molecules\n");
1667 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1668 &natoms, NULL, NULL, NULL, top);
1671 snew(grpname, ncos);
1672 snew(grpindex, ncos);
1673 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1674 ncos, gnx, grpindex, grpname);
1675 for (k = 0; (k < ncos); k++)
1677 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1678 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1682 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1683 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1684 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1685 opt2fn_null("-cos", NFILE, fnm),
1686 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1687 bPairs, corrtype[0],
1688 opt2fn("-c", NFILE, fnm),
1689 bGkr, opt2fn("-g", NFILE, fnm),
1690 bPhi, &nlevels, ndegrees,
1692 opt2fn("-cmap", NFILE, fnm), rcmax,
1693 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1694 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1695 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1697 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1698 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1699 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1700 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1701 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");