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.
45 #include "gromacs/legacyheaders/macros.h"
46 #include "gromacs/pbcutil/pbc.h"
47 #include "gromacs/utility/futil.h"
48 #include "gromacs/legacyheaders/viewit.h"
49 #include "gromacs/legacyheaders/txtdump.h"
50 #include "gromacs/statistics/statistics.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/random/random.h"
54 #include "gromacs/legacyheaders/names.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/legacyheaders/calcmu.h"
57 #include "gromacs/fileio/enxio.h"
59 #include "gromacs/legacyheaders/copyrite.h"
60 #include "gromacs/fileio/trxio.h"
62 #include "gromacs/bonded/bonded.h"
63 #include "gromacs/commandline/pargs.h"
64 #include "gromacs/fileio/matio.h"
65 #include "gromacs/fileio/xvgr.h"
66 #include "gromacs/linearalgebra/nrjac.h"
67 #include "gromacs/math/vec.h"
68 #include "gromacs/pbcutil/rmpbc.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/smalloc.h"
73 #define e2d(x) ENM2DEBYE*(x)
74 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
75 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
87 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
95 if ((ptr = getenv("GMX_DIPOLE_SPACING")) != NULL)
97 double bw = strtod(ptr, NULL);
102 gb->spacing = 0.01; /* nm */
104 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
111 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
114 snew(gb->elem, gb->nelem);
115 snew(gb->count, gb->nelem);
117 snew(gb->cmap, gb->nx);
118 gb->ny = std::max(2, ndegrees);
119 for (i = 0; (i < gb->nx); i++)
121 snew(gb->cmap[i], gb->ny);
128 static void done_gkrbin(t_gkrbin **gb)
136 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
138 int cy, index = gmx_nint(r/gb->spacing);
141 if (index < gb->nelem)
143 gb->elem[index] += cosa;
151 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
155 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
157 cy = std::min(gb->ny-1, std::max(0, cy));
160 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
162 gb->cmap[index][cy] += 1;
166 static void rvec2sprvec(rvec dipcart, rvec dipsp)
169 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
170 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
171 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
176 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
177 int mindex[], rvec x[], rvec mu[],
178 int ePBC, matrix box, t_atom *atom, int *nAtom)
180 static rvec *xcm[2] = { NULL, NULL};
181 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
182 real qtot, q, cosa, rr, phi;
186 for (n = 0; (n < ncos); n++)
190 snew(xcm[n], ngrp[n]);
192 for (i = 0; (i < ngrp[n]); i++)
194 /* Calculate center of mass of molecule */
200 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
205 clear_rvec(xcm[n][i]);
207 for (j = j0; j < j1; j++)
211 for (k = 0; k < DIM; k++)
213 xcm[n][i][k] += q*x[j][k];
216 svmul(1/qtot, xcm[n][i], xcm[n][i]);
220 set_pbc(&pbc, ePBC, box);
223 for (i = 0; i < ngrp[grp0]; i++)
225 gi = molindex[grp0][i];
226 j0 = (ncos == 2) ? 0 : i+1;
227 for (j = j0; j < ngrp[grp1]; j++)
229 gj = molindex[grp1][j];
230 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
232 /* Compute distance between molecules including PBC */
233 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
239 rvec r_ij, r_kj, r_kl, mm, nn;
243 copy_rvec(xcm[grp0][i], xj);
244 copy_rvec(xcm[grp1][j], xk);
245 rvec_add(xj, mu[gi], xi);
246 rvec_add(xk, mu[gj], xl);
247 phi = dih_angle(xi, xj, xk, xl, &pbc,
248 r_ij, r_kj, r_kl, mm, nn, /* out */
249 &sign, &t1, &t2, &t3);
254 cosa = cos_angle(mu[gi], mu[gj]);
257 if (debug || gmx_isnan(cosa))
259 fprintf(debug ? debug : stderr,
260 "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",
261 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
262 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
266 add2gkr(gb, rr, cosa, phi);
272 static real normalize_cmap(t_gkrbin *gb)
279 for (i = 0; (i < gb->nx); i++)
281 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
282 for (j = 0; (j < gb->ny); j++)
284 gb->cmap[i][j] /= vol;
285 hi = std::max(hi, gb->cmap[i][j]);
290 gmx_fatal(FARGS, "No data in the cmap");
295 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
302 t_rgb rlo = { 1, 1, 1 };
303 t_rgb rhi = { 0, 0, 0 };
305 hi = normalize_cmap(gb);
306 snew(xaxis, gb->nx+1);
307 for (i = 0; (i < gb->nx+1); i++)
309 xaxis[i] = i*gb->spacing;
312 for (j = 0; (j < gb->ny); j++)
316 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
320 yaxis[j] = (180.0*j)/(gb->ny-1.0);
322 /*2.0*j/(gb->ny-1.0)-1.0;*/
324 out = gmx_ffopen(cmap, "w");
326 "Dipole Orientation Distribution", "Fraction", "r (nm)",
327 gb->bPhi ? "Phi" : "Alpha",
328 gb->nx, gb->ny, xaxis, yaxis,
329 gb->cmap, 0, hi, rlo, rhi, nlevels);
335 static void print_gkrbin(const char *fn, t_gkrbin *gb,
336 int ngrp, int nframes, real volume,
337 const output_env_t oenv)
339 /* We compute Gk(r), gOO and hOO according to
340 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
341 * In this implementation the angle between dipoles is stored
342 * rather than their inner product. This allows to take polarizible
343 * models into account. The RDF is calculated as well, almost for free!
346 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
348 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
351 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
352 xvgr_legend(fp, asize(leg), leg, oenv);
354 Gkr = 1; /* Self-dipole inproduct = 1 */
359 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
360 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
364 while (last > 1 && gb->elem[last-1] == 0)
369 /* Divide by dipole squared, by number of frames, by number of origins.
370 * Multiply by 2 because we only take half the matrix of interactions
373 fac = 2.0/((double) ngrp * (double) nframes);
376 for (i = 0; i < last; i++)
378 /* Centre of the coordinate in the spherical layer */
381 /* Volume of the layer */
382 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
385 gOO = gb->count[i]*fac/(rho*vol_s);
387 /* Dipole correlation hOO, normalized by the relative number density, like
388 * in a Radial distribution function.
390 ggg = gb->elem[i]*fac;
391 hOO = 3.0*ggg/(rho*vol_s);
395 cosav = gb->elem[i]/gb->count[i];
401 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
403 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
404 x1, Gkr, cosav, hOO, gOO, ener);
412 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
413 real *t, int nre, t_enxframe *fr)
419 bCont = do_enx(fmu, fr);
422 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
423 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
428 if (Vol != -1) /* we've got Volume in the energy file */
430 *vol = fr->ener[Vol].e;
432 for (i = 0; i < DIM; i++)
434 mu[i] = fr->ener[iMu[i]].e;
442 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
445 int ncharged, m, a0, a1, a;
448 for (m = 0; m < n; m++)
450 a0 = mols->index[index[m]];
451 a1 = mols->index[index[m]+1];
454 for (a = a0; a < a1; a++)
459 /* This check is only for the count print */
460 if (fabs(qtot) > 0.01)
466 /* Remove the net charge at the center of mass */
467 for (a = a0; a < a1; a++)
469 atom[a].q -= qtot*atom[a].m/mtot;
476 printf("There are %d charged molecules in the selection,\n"
477 "will subtract their charge at their center of mass\n", ncharged);
481 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
487 for (k = k0; (k < k1); k++)
490 for (m = 0; (m < DIM); m++)
497 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
499 int i, k, m, n, niter;
500 real q, r2, mass, masstot;
501 rvec com; /* center of mass */
502 rvec r; /* distance of atoms to center of mass */
504 double dd[DIM], **ev, tmp;
508 for (i = 0; (i < DIM); i++)
515 /* Compute center of mass */
518 for (k = k0; (k < k1); k++)
522 for (i = 0; (i < DIM); i++)
524 com[i] += mass*x[k][i];
527 svmul((1.0/masstot), com, com);
529 /* We want traceless quadrupole moments, so let us calculate the complete
530 * quadrupole moment tensor and diagonalize this tensor to get
531 * the individual components on the diagonal.
534 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
536 for (m = 0; (m < DIM); m++)
538 for (n = 0; (n < DIM); n++)
543 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
545 q = (atom[k].q)*100.0;
546 rvec_sub(x[k], com, r);
548 for (m = 0; (m < DIM); m++)
550 for (n = 0; (n < DIM); n++)
552 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
558 for (i = 0; (i < DIM); i++)
560 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
561 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
565 /* We've got the quadrupole tensor, now diagonalize the sucker */
566 jacobi(inten, 3, dd, ev, &niter);
570 for (i = 0; (i < DIM); i++)
572 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
573 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
575 for (i = 0; (i < DIM); i++)
577 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
578 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
581 /* Sort the eigenvalues, for water we know that the order is as follows:
585 * At the moment I have no idea how this will work out for other molecules...
589 if (dd[i+1] > dd[i]) { \
598 for (m = 0; (m < DIM); m++)
600 quad[0] = dd[2]; /* yy */
601 quad[1] = dd[0]; /* zz */
602 quad[2] = dd[1]; /* xx */
607 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
611 for (i = 0; (i < DIM); i++)
621 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
623 real calc_eps(double M_diff, double volume, double epsRF, double temp)
625 double eps, A, teller, noemer;
626 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
627 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
629 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
638 teller = 1 + (A*2*epsRF/(2*epsRF+1));
639 noemer = 1 - (A/(2*epsRF+1));
641 eps = teller / noemer;
646 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
647 int idim, int nslice, rvec slab_dipole[],
653 for (k = k0; (k < k1); k++)
658 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
659 rvec_inc(slab_dipole[k], mu);
662 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
663 rvec slab_dipole[], matrix box, int nframes,
664 const output_env_t oenv)
670 const char *leg_dim[4] = {
671 "\\f{12}m\\f{4}\\sX\\N",
672 "\\f{12}m\\f{4}\\sY\\N",
673 "\\f{12}m\\f{4}\\sZ\\N",
674 "\\f{12}m\\f{4}\\stot\\N"
677 sprintf(buf, "Box-%c (nm)", 'X'+idim);
678 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
680 xvgr_legend(fp, DIM, leg_dim, oenv);
681 for (i = 0; (i < nslice); i++)
683 mutot = norm(slab_dipole[i])/nframes;
684 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
685 ((i+0.5)*box[idim][idim])/nslice,
686 slab_dipole[i][XX]/nframes,
687 slab_dipole[i][YY]/nframes,
688 slab_dipole[i][ZZ]/nframes,
692 do_view(oenv, fn, "-autoscale xy -nxy");
695 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
698 double dc, d, ddc1, ddc2, ddc3;
699 rvec xxx = { 1, 0, 0 };
700 rvec yyy = { 0, 1, 0 };
701 rvec zzz = { 0, 0, 1 };
704 ddc1 = ddc2 = ddc3 = 0;
705 for (i = k = 0; (i < n); i++)
707 ddc1 += fabs(cos_angle(dip[i], xxx));
708 ddc2 += fabs(cos_angle(dip[i], yyy));
709 ddc3 += fabs(cos_angle(dip[i], zzz));
712 for (j = i+1; (j < n); j++, k++)
714 dc = cos_angle(dip[i], dip[j]);
725 static void do_dip(t_topology *top, int ePBC, real volume,
727 const char *out_mtot, const char *out_eps,
728 const char *out_aver, const char *dipdist,
729 const char *cosaver, const char *fndip3d,
730 const char *fnadip, gmx_bool bPairs,
731 const char *corrtype, const char *corf,
732 gmx_bool bGkr, const char *gkrfn,
733 gmx_bool bPhi, int *nlevels, int ndegrees,
735 const char *cmap, real rcmax,
737 gmx_bool bMU, const char *mufn,
738 int *gnx, int *molindex[],
739 real mu_max, real mu_aver,
740 real epsilonRF, real temp,
741 int *gkatom, int skip,
742 gmx_bool bSlab, int nslices,
743 const char *axtitle, const char *slabfn,
744 const output_env_t oenv)
746 const char *leg_mtot[] = {
752 #define NLEGMTOT asize(leg_mtot)
753 const char *leg_eps[] = {
758 #define NLEGEPS asize(leg_eps)
759 const char *leg_aver[] = {
762 "< |M|\\S2\\N > - < |M| >\\S2\\N",
763 "< |M| >\\S2\\N / < |M|\\S2\\N >"
765 #define NLEGAVER asize(leg_aver)
766 const char *leg_cosaver[] = {
767 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
769 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
770 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
771 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
773 #define NLEGCOSAVER asize(leg_cosaver)
774 const char *leg_adip[] = {
779 #define NLEGADIP asize(leg_adip)
781 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
782 FILE *dip3d = NULL, *adip = NULL;
783 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
784 t_gkrbin *gkrbin = NULL;
785 gmx_enxnm_t *enm = NULL;
787 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
788 ener_file_t fmu = NULL;
789 int i, n, m, natom = 0, gnx_tot, teller, tel3;
791 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
793 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
796 gmx_bool bCorr, bTotal, bCont;
797 double M_diff = 0, epsilon, invtel, vol_aver;
798 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
799 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
800 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
803 rvec *slab_dipoles = NULL;
805 t_block *mols = NULL;
806 gmx_rmpbc_t gpbc = NULL;
818 /* initialize to a negative value so we can see that it wasn't picked up */
819 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
822 fmu = open_enx(mufn, "r");
823 do_enxnms(fmu, &nre, &enm);
825 /* Determine the indexes of the energy grps we need */
826 for (i = 0; (i < nre); i++)
828 if (strstr(enm[i].name, "Volume"))
832 else if (strstr(enm[i].name, "Mu-X"))
836 else if (strstr(enm[i].name, "Mu-Y"))
840 else if (strstr(enm[i].name, "Mu-Z"))
845 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
847 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
852 atom = top->atoms.atom;
855 if ((iVol == -1) && bMU)
857 printf("Using Volume from topology: %g nm^3\n", volume);
860 /* Correlation stuff */
861 bCorr = (corrtype[0] != 'n');
862 bTotal = (corrtype[0] == 't');
868 snew(muall[0], nframes*DIM);
873 for (i = 0; (i < gnx[0]); i++)
875 snew(muall[i], nframes*DIM);
880 /* Allocate array which contains for every molecule in a frame the
885 snew(dipole, gnx_tot);
890 for (i = 0; (i < DIM); i++)
892 Qlsq[i] = gmx_stats_init();
894 mulsq = gmx_stats_init();
896 /* Open all the files */
897 outmtot = xvgropen(out_mtot,
898 "Total dipole moment of the simulation box vs. time",
899 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
900 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
901 "Time (ps)", "", oenv);
902 outaver = xvgropen(out_aver, "Total dipole moment",
903 "Time (ps)", "D", oenv);
906 idim = axtitle[0] - 'X';
907 if ((idim < 0) || (idim >= DIM))
909 idim = axtitle[0] - 'x';
911 if ((idim < 0) || (idim >= DIM))
919 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
920 axtitle, nslices, idim);
923 snew(slab_dipoles, nslices);
924 fprintf(stderr, "Doing slab analysis\n");
930 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
931 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
936 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
937 "Average absolute dipole orientation", "Time (ps)", "", oenv);
938 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
944 snew(dipsp, gnx_tot);
946 /* we need a dummy file for gnuplot */
947 dip3d = (FILE *)gmx_ffopen("dummy.dat", "w");
948 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
951 dip3d = (FILE *)gmx_ffopen(fndip3d, "w");
954 gmx::BinaryInformationSettings settings;
955 settings.generatedByHeader(true);
956 settings.linePrefix("# ");
957 gmx::printBinaryInformation(dip3d, output_env_get_program_context(oenv),
960 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
963 /* Write legends to all the files */
964 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
965 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
967 if (bMU && (mu_aver == -1))
969 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
973 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
979 /* Read the first frame from energy or traj file */
984 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
987 timecheck = check_times(t);
992 if ((teller % 10) == 0)
994 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
999 printf("End of %s reached\n", mufn);
1003 while (bCont && (timecheck < 0));
1007 natom = read_first_x(oenv, &status, fn, &t, &x, box);
1010 /* Calculate spacing for dipole bin (simple histogram) */
1011 ndipbin = 1 + static_cast<int>(mu_max/0.01);
1012 snew(dipole_bin, ndipbin);
1014 for (m = 0; (m < DIM); m++)
1016 M[m] = M2[m] = M4[m] = 0.0;
1021 /* Use 0.7 iso 0.5 to account for pressure scaling */
1022 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1023 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1025 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1027 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1029 /* Start while loop over frames */
1034 if (bCorr && (teller >= nframes))
1039 srenew(muall[0], nframes*DIM);
1043 for (i = 0; (i < gnx_tot); i++)
1045 srenew(muall[i], nframes*DIM);
1051 muframelsq = gmx_stats_init();
1054 for (m = 0; (m < DIM); m++)
1061 /* Copy rvec into double precision local variable */
1062 for (m = 0; (m < DIM); m++)
1070 for (m = 0; (m < DIM); m++)
1075 gmx_rmpbc(gpbc, natom, box, x);
1077 /* Begin loop of all molecules in frame */
1078 for (n = 0; (n < ncos); n++)
1080 for (i = 0; (i < gnx[n]); i++)
1084 ind0 = mols->index[molindex[n][i]];
1085 ind1 = mols->index[molindex[n][i]+1];
1087 mol_dip(ind0, ind1, x, atom, dipole[i]);
1088 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1089 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1092 update_slab_dipoles(ind0, ind1, x,
1093 dipole[i], idim, nslices, slab_dipoles, box);
1097 mol_quad(ind0, ind1, x, atom, quad);
1098 for (m = 0; (m < DIM); m++)
1100 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1103 if (bCorr && !bTotal)
1106 muall[i][tel3+XX] = dipole[i][XX];
1107 muall[i][tel3+YY] = dipole[i][YY];
1108 muall[i][tel3+ZZ] = dipole[i][ZZ];
1111 for (m = 0; (m < DIM); m++)
1113 M_av[m] += dipole[i][m]; /* M per frame */
1114 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1116 mu_mol = sqrt(mu_mol);
1118 mu_ave += mu_mol; /* calc. the average mu */
1120 /* Update the dipole distribution */
1121 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1129 rvec2sprvec(dipole[i], dipsp[i]);
1131 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1133 if (dipsp[i][ZZ] < 0.5 * M_PI)
1142 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1144 if (dipsp[i][ZZ] < 0.5 * M_PI)
1153 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1155 if (dipsp[i][ZZ] < 0.5 * M_PI)
1164 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1166 if (dipsp[i][ZZ] < 0.5 * M_PI)
1177 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1182 x[ind0][XX]+dipole[i][XX]/25,
1183 x[ind0][YY]+dipole[i][YY]/25,
1184 x[ind0][ZZ]+dipole[i][ZZ]/25,
1188 } /* End loop of all molecules in frame */
1192 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1193 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1194 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1195 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1196 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1197 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1201 /* Compute square of total dipole */
1202 for (m = 0; (m < DIM); m++)
1204 M_av2[m] = M_av[m]*M_av[m];
1209 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1210 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1211 sqr(dipaxis[YY]-0.5)+
1212 sqr(dipaxis[ZZ]-0.5));
1215 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1216 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1220 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1221 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1227 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1234 muall[0][tel3+XX] = M_av[XX];
1235 muall[0][tel3+YY] = M_av[YY];
1236 muall[0][tel3+ZZ] = M_av[ZZ];
1239 /* Write to file the total dipole moment of the box, and its components
1242 if ((skip == 0) || ((teller % skip) == 0))
1244 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1245 t, M_av[XX], M_av[YY], M_av[ZZ],
1246 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1249 for (m = 0; (m < DIM); m++)
1253 M4[m] += sqr(M_av2[m]);
1255 /* Increment loop counter */
1258 /* Calculate for output the running averages */
1259 invtel = 1.0/teller;
1260 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1261 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1262 M_diff = M2_ave - M_ave2;
1264 /* Compute volume from box in traj, else we use the one from above */
1271 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1273 /* Calculate running average for dipole */
1276 mu_aver = (mu_ave/gnx_tot)*invtel;
1279 if ((skip == 0) || ((teller % skip) == 0))
1281 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1282 * the two. Here M is sum mu_i. Further write the finite system
1283 * Kirkwood G factor and epsilon.
1285 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1286 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1291 gmx_stats_get_average(muframelsq, &aver);
1292 fprintf(adip, "%10g %f \n", t, aver);
1295 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1297 if (!bMU || (mu_aver != -1))
1299 /* Finite system Kirkwood G-factor */
1300 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1301 /* Infinite system Kirkwood G-factor */
1302 if (epsilonRF == 0.0)
1304 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1308 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1309 Gk/(3*epsilon*(2*epsilonRF+1)));
1312 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1317 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1320 gmx_stats_done(muframelsq);
1324 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1328 bCont = read_next_x(oenv, status, &t, x, box);
1330 timecheck = check_times(t);
1332 while (bCont && (timecheck == 0) );
1334 gmx_rmpbc_done(gpbc);
1341 gmx_ffclose(outmtot);
1342 gmx_ffclose(outaver);
1343 gmx_ffclose(outeps);
1357 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1358 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1359 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1360 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1361 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1367 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1368 sfree(slab_dipoles);
1372 printf("Average volume over run is %g\n", vol_aver);
1375 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1376 print_cmap(cmap, gkrbin, nlevels);
1378 /* Autocorrelation function */
1383 printf("Not enough frames for autocorrelation\n");
1387 dt = (t1 - t0)/(teller-1);
1388 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1394 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1395 teller, 1, muall, dt, mode, TRUE);
1399 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1400 teller, gnx_tot, muall, dt,
1401 mode, strcmp(corrtype, "molsep"));
1407 real aver, sigma, error;
1409 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1410 printf("\nDipole moment (Debye)\n");
1411 printf("---------------------\n");
1412 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1413 aver, sigma, error);
1417 for (m = 0; (m < DIM); m++)
1419 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1422 printf("\nQuadrupole moment (Debye-Ang)\n");
1423 printf("-----------------------------\n");
1424 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1425 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1426 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1430 printf("The following averages for the complete trajectory have been calculated:\n\n");
1431 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1432 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1433 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1435 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1436 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1437 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1439 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1440 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1442 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1444 if (!bMU || (mu_aver != -1))
1446 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1447 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1449 printf("Epsilon = %g\n", epsilon);
1453 /* Write to file the dipole moment distibution during the simulation.
1455 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1456 for (i = 0; (i < ndipbin); i++)
1458 fprintf(outdd, "%10g %10f\n",
1459 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1466 done_gkrbin(&gkrbin);
1470 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1479 while (m < mols->nr && index[i] != mols->index[m])
1485 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1487 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1489 if (i >= *n || index[i] != j)
1491 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1495 /* Modify the index in place */
1498 printf("There are %d molecules in the selection\n", nmol);
1502 int gmx_dipoles(int argc, char *argv[])
1504 const char *desc[] = {
1505 "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1506 "system. From this you can compute e.g. the dielectric constant for",
1507 "low-dielectric media.",
1508 "For molecules with a net charge, the net charge is subtracted at",
1509 "center of mass of the molecule.[PAR]",
1510 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1511 "components as well as the norm of the vector.",
1512 "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",
1514 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1516 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1517 "Furthermore, the dipole autocorrelation function will be computed when",
1518 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1520 "The correlation functions can be averaged over all molecules",
1521 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1522 "or it can be computed over the total dipole moment of the simulation box",
1523 "([TT]total[tt]).[PAR]",
1524 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1525 "G-factor, as well as the average cosine of the angle between the dipoles",
1526 "as a function of the distance. The plot also includes gOO and hOO",
1527 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1528 "we also include the energy per scale computed by taking the inner product of",
1529 "the dipoles divided by the distance to the third power.[PAR]",
1532 "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1533 "This will calculate the autocorrelation function of the molecular",
1534 "dipoles using a first order Legendre polynomial of the angle of the",
1535 "dipole vector and itself a time t later. For this calculation 1001",
1536 "frames will be used. Further, the dielectric constant will be calculated",
1537 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1538 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1539 "distribution function a maximum of 5.0 will be used."
1541 real mu_max = 5, mu_aver = -1, rcmax = 0;
1542 real epsilonRF = 0.0, temp = 300;
1543 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1544 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1545 const char *axtitle = "Z";
1546 int nslices = 10; /* nr of slices defined */
1547 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1548 int nlevels = 20, ndegrees = 90;
1551 { "-mu", FALSE, etREAL, {&mu_aver},
1552 "dipole of a single molecule (in Debye)" },
1553 { "-mumax", FALSE, etREAL, {&mu_max},
1554 "max dipole in Debye (for histogram)" },
1555 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1556 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1557 { "-skip", FALSE, etINT, {&skip},
1558 "Skip steps in the output (but not in the computations)" },
1559 { "-temp", FALSE, etREAL, {&temp},
1560 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1561 { "-corr", FALSE, etENUM, {corrtype},
1562 "Correlation function to calculate" },
1563 { "-pairs", FALSE, etBOOL, {&bPairs},
1564 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1565 { "-quad", FALSE, etBOOL, {&bQuad},
1566 "Take quadrupole into account"},
1567 { "-ncos", FALSE, etINT, {&ncos},
1568 "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." },
1569 { "-axis", FALSE, etSTR, {&axtitle},
1570 "Take the normal on the computational box in direction X, Y or Z." },
1571 { "-sl", FALSE, etINT, {&nslices},
1572 "Divide the box into this number of slices." },
1573 { "-gkratom", FALSE, etINT, {&nFA},
1574 "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" },
1575 { "-gkratom2", FALSE, etINT, {&nFB},
1576 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1577 { "-rcmax", FALSE, etREAL, {&rcmax},
1578 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1579 { "-phi", FALSE, etBOOL, {&bPhi},
1580 "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." },
1581 { "-nlevels", FALSE, etINT, {&nlevels},
1582 "Number of colors in the cmap output" },
1583 { "-ndegrees", FALSE, etINT, {&ndegrees},
1584 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1589 char **grpname = NULL;
1590 gmx_bool bGkr, bMU, bSlab;
1592 { efEDR, "-en", NULL, ffOPTRD },
1593 { efTRX, "-f", NULL, ffREAD },
1594 { efTPX, NULL, NULL, ffREAD },
1595 { efNDX, NULL, NULL, ffOPTRD },
1596 { efXVG, "-o", "Mtot", ffWRITE },
1597 { efXVG, "-eps", "epsilon", ffWRITE },
1598 { efXVG, "-a", "aver", ffWRITE },
1599 { efXVG, "-d", "dipdist", ffWRITE },
1600 { efXVG, "-c", "dipcorr", ffOPTWR },
1601 { efXVG, "-g", "gkr", ffOPTWR },
1602 { efXVG, "-adip", "adip", ffOPTWR },
1603 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1604 { efXVG, "-cos", "cosaver", ffOPTWR },
1605 { efXPM, "-cmap", "cmap", ffOPTWR },
1606 { efXVG, "-slab", "slab", ffOPTWR }
1608 #define NFILE asize(fnm)
1617 ppa = add_acf_pargs(&npargs, pa);
1618 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
1619 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1624 printf("Using %g as mu_max and %g as the dipole moment.\n",
1626 if (epsilonRF == 0.0)
1628 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1631 bMU = opt2bSet("-en", NFILE, fnm);
1634 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.");
1636 bGkr = opt2bSet("-g", NFILE, fnm);
1637 if (opt2parg_bSet("-ncos", asize(pa), pa))
1639 if ((ncos != 1) && (ncos != 2))
1641 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1645 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1646 opt2parg_bSet("-axis", asize(pa), pa));
1651 printf("WARNING: Can not determine quadrupoles from energy file\n");
1656 printf("WARNING: Can not determine Gk(r) from energy file\n");
1662 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1663 " not enter a valid dipole for the molecules\n");
1668 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1669 &natoms, NULL, NULL, NULL, top);
1672 snew(grpname, ncos);
1673 snew(grpindex, ncos);
1674 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1675 ncos, gnx, grpindex, grpname);
1676 for (k = 0; (k < ncos); k++)
1678 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1679 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1683 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1684 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1685 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1686 opt2fn_null("-cos", NFILE, fnm),
1687 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1688 bPairs, corrtype[0],
1689 opt2fn("-c", NFILE, fnm),
1690 bGkr, opt2fn("-g", NFILE, fnm),
1691 bPhi, &nlevels, ndegrees,
1693 opt2fn("-cmap", NFILE, fnm), rcmax,
1694 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1695 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1696 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1698 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1699 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1700 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1701 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1702 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");