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.
46 #include "gromacs/pbcutil/pbc.h"
48 #include "gromacs/utility/futil.h"
51 #include "gromacs/statistics/statistics.h"
53 #include "gromacs/topology/index.h"
54 #include "gromacs/random/random.h"
56 #include "gromacs/math/units.h"
58 #include "gromacs/fileio/enxio.h"
61 #include "gromacs/fileio/trxio.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/programcontext.h"
72 #include "gromacs/utility/smalloc.h"
74 #define e2d(x) ENM2DEBYE*(x)
75 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
76 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
88 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
96 if ((ptr = getenv("GMX_DIPOLE_SPACING")) != NULL)
98 double bw = strtod(ptr, NULL);
103 gb->spacing = 0.01; /* nm */
105 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
112 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
115 snew(gb->elem, gb->nelem);
116 snew(gb->count, gb->nelem);
118 snew(gb->cmap, gb->nx);
119 gb->ny = std::max(2, ndegrees);
120 for (i = 0; (i < gb->nx); i++)
122 snew(gb->cmap[i], gb->ny);
129 static void done_gkrbin(t_gkrbin **gb)
137 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
139 int cy, index = gmx_nint(r/gb->spacing);
142 if (index < gb->nelem)
144 gb->elem[index] += cosa;
152 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
156 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
158 cy = std::min(gb->ny-1, std::max(0, cy));
161 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
163 gb->cmap[index][cy] += 1;
167 static void rvec2sprvec(rvec dipcart, rvec dipsp)
170 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
171 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
172 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
177 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
178 int mindex[], rvec x[], rvec mu[],
179 int ePBC, matrix box, t_atom *atom, int *nAtom)
181 static rvec *xcm[2] = { NULL, NULL};
182 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
183 real qtot, q, cosa, rr, phi;
187 for (n = 0; (n < ncos); n++)
191 snew(xcm[n], ngrp[n]);
193 for (i = 0; (i < ngrp[n]); i++)
195 /* Calculate center of mass of molecule */
201 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
206 clear_rvec(xcm[n][i]);
208 for (j = j0; j < j1; j++)
212 for (k = 0; k < DIM; k++)
214 xcm[n][i][k] += q*x[j][k];
217 svmul(1/qtot, xcm[n][i], xcm[n][i]);
221 set_pbc(&pbc, ePBC, box);
224 for (i = 0; i < ngrp[grp0]; i++)
226 gi = molindex[grp0][i];
227 j0 = (ncos == 2) ? 0 : i+1;
228 for (j = j0; j < ngrp[grp1]; j++)
230 gj = molindex[grp1][j];
231 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
233 /* Compute distance between molecules including PBC */
234 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
240 rvec r_ij, r_kj, r_kl, mm, nn;
244 copy_rvec(xcm[grp0][i], xj);
245 copy_rvec(xcm[grp1][j], xk);
246 rvec_add(xj, mu[gi], xi);
247 rvec_add(xk, mu[gj], xl);
248 phi = dih_angle(xi, xj, xk, xl, &pbc,
249 r_ij, r_kj, r_kl, mm, nn, /* out */
250 &sign, &t1, &t2, &t3);
255 cosa = cos_angle(mu[gi], mu[gj]);
258 if (debug || gmx_isnan(cosa))
260 fprintf(debug ? debug : stderr,
261 "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",
262 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
263 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
267 add2gkr(gb, rr, cosa, phi);
273 static real normalize_cmap(t_gkrbin *gb)
280 for (i = 0; (i < gb->nx); i++)
282 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
283 for (j = 0; (j < gb->ny); j++)
285 gb->cmap[i][j] /= vol;
286 hi = std::max(hi, gb->cmap[i][j]);
291 gmx_fatal(FARGS, "No data in the cmap");
296 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
303 t_rgb rlo = { 1, 1, 1 };
304 t_rgb rhi = { 0, 0, 0 };
306 hi = normalize_cmap(gb);
307 snew(xaxis, gb->nx+1);
308 for (i = 0; (i < gb->nx+1); i++)
310 xaxis[i] = i*gb->spacing;
313 for (j = 0; (j < gb->ny); j++)
317 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
321 yaxis[j] = (180.0*j)/(gb->ny-1.0);
323 /*2.0*j/(gb->ny-1.0)-1.0;*/
325 out = gmx_ffopen(cmap, "w");
327 "Dipole Orientation Distribution", "Fraction", "r (nm)",
328 gb->bPhi ? "Phi" : "Alpha",
329 gb->nx, gb->ny, xaxis, yaxis,
330 gb->cmap, 0, hi, rlo, rhi, nlevels);
336 static void print_gkrbin(const char *fn, t_gkrbin *gb,
337 int ngrp, int nframes, real volume,
338 const output_env_t oenv)
340 /* We compute Gk(r), gOO and hOO according to
341 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
342 * In this implementation the angle between dipoles is stored
343 * rather than their inner product. This allows to take polarizible
344 * models into account. The RDF is calculated as well, almost for free!
347 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
349 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
352 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
353 xvgr_legend(fp, asize(leg), leg, oenv);
355 Gkr = 1; /* Self-dipole inproduct = 1 */
360 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
361 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
365 while (last > 1 && gb->elem[last-1] == 0)
370 /* Divide by dipole squared, by number of frames, by number of origins.
371 * Multiply by 2 because we only take half the matrix of interactions
374 fac = 2.0/((double) ngrp * (double) nframes);
377 for (i = 0; i < last; i++)
379 /* Centre of the coordinate in the spherical layer */
382 /* Volume of the layer */
383 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
386 gOO = gb->count[i]*fac/(rho*vol_s);
388 /* Dipole correlation hOO, normalized by the relative number density, like
389 * in a Radial distribution function.
391 ggg = gb->elem[i]*fac;
392 hOO = 3.0*ggg/(rho*vol_s);
396 cosav = gb->elem[i]/gb->count[i];
402 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
404 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
405 x1, Gkr, cosav, hOO, gOO, ener);
413 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
414 real *t, int nre, t_enxframe *fr)
420 bCont = do_enx(fmu, fr);
423 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
424 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
429 if (Vol != -1) /* we've got Volume in the energy file */
431 *vol = fr->ener[Vol].e;
433 for (i = 0; i < DIM; i++)
435 mu[i] = fr->ener[iMu[i]].e;
443 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
446 int ncharged, m, a0, a1, a;
449 for (m = 0; m < n; m++)
451 a0 = mols->index[index[m]];
452 a1 = mols->index[index[m]+1];
455 for (a = a0; a < a1; a++)
460 /* This check is only for the count print */
461 if (fabs(qtot) > 0.01)
467 /* Remove the net charge at the center of mass */
468 for (a = a0; a < a1; a++)
470 atom[a].q -= qtot*atom[a].m/mtot;
477 printf("There are %d charged molecules in the selection,\n"
478 "will subtract their charge at their center of mass\n", ncharged);
482 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
488 for (k = k0; (k < k1); k++)
491 for (m = 0; (m < DIM); m++)
498 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
500 int i, k, m, n, niter;
501 real q, r2, mass, masstot;
502 rvec com; /* center of mass */
503 rvec r; /* distance of atoms to center of mass */
505 double dd[DIM], **ev, tmp;
509 for (i = 0; (i < DIM); i++)
516 /* Compute center of mass */
519 for (k = k0; (k < k1); k++)
523 for (i = 0; (i < DIM); i++)
525 com[i] += mass*x[k][i];
528 svmul((1.0/masstot), com, com);
530 /* We want traceless quadrupole moments, so let us calculate the complete
531 * quadrupole moment tensor and diagonalize this tensor to get
532 * the individual components on the diagonal.
535 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
537 for (m = 0; (m < DIM); m++)
539 for (n = 0; (n < DIM); n++)
544 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
546 q = (atom[k].q)*100.0;
547 rvec_sub(x[k], com, r);
549 for (m = 0; (m < DIM); m++)
551 for (n = 0; (n < DIM); n++)
553 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
559 for (i = 0; (i < DIM); i++)
561 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
562 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
566 /* We've got the quadrupole tensor, now diagonalize the sucker */
567 jacobi(inten, 3, dd, ev, &niter);
571 for (i = 0; (i < DIM); i++)
573 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
574 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
576 for (i = 0; (i < DIM); i++)
578 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
579 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
582 /* Sort the eigenvalues, for water we know that the order is as follows:
586 * At the moment I have no idea how this will work out for other molecules...
590 if (dd[i+1] > dd[i]) { \
599 for (m = 0; (m < DIM); m++)
601 quad[0] = dd[2]; /* yy */
602 quad[1] = dd[0]; /* zz */
603 quad[2] = dd[1]; /* xx */
608 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
612 for (i = 0; (i < DIM); i++)
622 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
624 real calc_eps(double M_diff, double volume, double epsRF, double temp)
626 double eps, A, teller, noemer;
627 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
628 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
630 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
639 teller = 1 + (A*2*epsRF/(2*epsRF+1));
640 noemer = 1 - (A/(2*epsRF+1));
642 eps = teller / noemer;
647 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
648 int idim, int nslice, rvec slab_dipole[],
654 for (k = k0; (k < k1); k++)
659 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
660 rvec_inc(slab_dipole[k], mu);
663 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
664 rvec slab_dipole[], matrix box, int nframes,
665 const output_env_t oenv)
671 const char *leg_dim[4] = {
672 "\\f{12}m\\f{4}\\sX\\N",
673 "\\f{12}m\\f{4}\\sY\\N",
674 "\\f{12}m\\f{4}\\sZ\\N",
675 "\\f{12}m\\f{4}\\stot\\N"
678 sprintf(buf, "Box-%c (nm)", 'X'+idim);
679 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
681 xvgr_legend(fp, DIM, leg_dim, oenv);
682 for (i = 0; (i < nslice); i++)
684 mutot = norm(slab_dipole[i])/nframes;
685 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
686 ((i+0.5)*box[idim][idim])/nslice,
687 slab_dipole[i][XX]/nframes,
688 slab_dipole[i][YY]/nframes,
689 slab_dipole[i][ZZ]/nframes,
693 do_view(oenv, fn, "-autoscale xy -nxy");
696 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
699 double dc, d, ddc1, ddc2, ddc3;
700 rvec xxx = { 1, 0, 0 };
701 rvec yyy = { 0, 1, 0 };
702 rvec zzz = { 0, 0, 1 };
705 ddc1 = ddc2 = ddc3 = 0;
706 for (i = k = 0; (i < n); i++)
708 ddc1 += fabs(cos_angle(dip[i], xxx));
709 ddc2 += fabs(cos_angle(dip[i], yyy));
710 ddc3 += fabs(cos_angle(dip[i], zzz));
713 for (j = i+1; (j < n); j++, k++)
715 dc = cos_angle(dip[i], dip[j]);
726 static void do_dip(t_topology *top, int ePBC, real volume,
728 const char *out_mtot, const char *out_eps,
729 const char *out_aver, const char *dipdist,
730 const char *cosaver, const char *fndip3d,
731 const char *fnadip, gmx_bool bPairs,
732 const char *corrtype, const char *corf,
733 gmx_bool bGkr, const char *gkrfn,
734 gmx_bool bPhi, int *nlevels, int ndegrees,
736 const char *cmap, real rcmax,
738 gmx_bool bMU, const char *mufn,
739 int *gnx, int *molindex[],
740 real mu_max, real mu_aver,
741 real epsilonRF, real temp,
742 int *gkatom, int skip,
743 gmx_bool bSlab, int nslices,
744 const char *axtitle, const char *slabfn,
745 const output_env_t oenv)
747 const char *leg_mtot[] = {
753 #define NLEGMTOT asize(leg_mtot)
754 const char *leg_eps[] = {
759 #define NLEGEPS asize(leg_eps)
760 const char *leg_aver[] = {
763 "< |M|\\S2\\N > - < |M| >\\S2\\N",
764 "< |M| >\\S2\\N / < |M|\\S2\\N >"
766 #define NLEGAVER asize(leg_aver)
767 const char *leg_cosaver[] = {
768 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
770 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
771 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
772 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
774 #define NLEGCOSAVER asize(leg_cosaver)
775 const char *leg_adip[] = {
780 #define NLEGADIP asize(leg_adip)
782 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
783 FILE *dip3d = NULL, *adip = NULL;
784 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
785 t_gkrbin *gkrbin = NULL;
786 gmx_enxnm_t *enm = NULL;
788 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
789 ener_file_t fmu = NULL;
790 int i, n, m, natom = 0, gnx_tot, teller, tel3;
792 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
794 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
797 gmx_bool bCorr, bTotal, bCont;
798 double M_diff = 0, epsilon, invtel, vol_aver;
799 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
800 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
801 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
804 rvec *slab_dipoles = NULL;
806 t_block *mols = NULL;
807 gmx_rmpbc_t gpbc = NULL;
819 /* initialize to a negative value so we can see that it wasn't picked up */
820 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
823 fmu = open_enx(mufn, "r");
824 do_enxnms(fmu, &nre, &enm);
826 /* Determine the indexes of the energy grps we need */
827 for (i = 0; (i < nre); i++)
829 if (strstr(enm[i].name, "Volume"))
833 else if (strstr(enm[i].name, "Mu-X"))
837 else if (strstr(enm[i].name, "Mu-Y"))
841 else if (strstr(enm[i].name, "Mu-Z"))
846 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
848 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
853 atom = top->atoms.atom;
856 if ((iVol == -1) && bMU)
858 printf("Using Volume from topology: %g nm^3\n", volume);
861 /* Correlation stuff */
862 bCorr = (corrtype[0] != 'n');
863 bTotal = (corrtype[0] == 't');
869 snew(muall[0], nframes*DIM);
874 for (i = 0; (i < gnx[0]); i++)
876 snew(muall[i], nframes*DIM);
881 /* Allocate array which contains for every molecule in a frame the
886 snew(dipole, gnx_tot);
891 for (i = 0; (i < DIM); i++)
893 Qlsq[i] = gmx_stats_init();
895 mulsq = gmx_stats_init();
897 /* Open all the files */
898 outmtot = xvgropen(out_mtot,
899 "Total dipole moment of the simulation box vs. time",
900 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
901 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
902 "Time (ps)", "", oenv);
903 outaver = xvgropen(out_aver, "Total dipole moment",
904 "Time (ps)", "D", oenv);
907 idim = axtitle[0] - 'X';
908 if ((idim < 0) || (idim >= DIM))
910 idim = axtitle[0] - 'x';
912 if ((idim < 0) || (idim >= DIM))
920 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
921 axtitle, nslices, idim);
924 snew(slab_dipoles, nslices);
925 fprintf(stderr, "Doing slab analysis\n");
931 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
932 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
937 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
938 "Average absolute dipole orientation", "Time (ps)", "", oenv);
939 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
945 snew(dipsp, gnx_tot);
947 /* we need a dummy file for gnuplot */
948 dip3d = (FILE *)gmx_ffopen("dummy.dat", "w");
949 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
952 dip3d = (FILE *)gmx_ffopen(fndip3d, "w");
955 gmx::BinaryInformationSettings settings;
956 settings.generatedByHeader(true);
957 settings.linePrefix("# ");
958 gmx::printBinaryInformation(dip3d, gmx::getProgramContext(),
961 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
964 /* Write legends to all the files */
965 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
966 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
968 if (bMU && (mu_aver == -1))
970 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
974 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
980 /* Read the first frame from energy or traj file */
985 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
988 timecheck = check_times(t);
993 if ((teller % 10) == 0)
995 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
1000 printf("End of %s reached\n", mufn);
1004 while (bCont && (timecheck < 0));
1008 natom = read_first_x(oenv, &status, fn, &t, &x, box);
1011 /* Calculate spacing for dipole bin (simple histogram) */
1012 ndipbin = 1 + static_cast<int>(mu_max/0.01);
1013 snew(dipole_bin, ndipbin);
1015 for (m = 0; (m < DIM); m++)
1017 M[m] = M2[m] = M4[m] = 0.0;
1022 /* Use 0.7 iso 0.5 to account for pressure scaling */
1023 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1024 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1026 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1028 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1030 /* Start while loop over frames */
1035 if (bCorr && (teller >= nframes))
1040 srenew(muall[0], nframes*DIM);
1044 for (i = 0; (i < gnx_tot); i++)
1046 srenew(muall[i], nframes*DIM);
1052 muframelsq = gmx_stats_init();
1055 for (m = 0; (m < DIM); m++)
1062 /* Copy rvec into double precision local variable */
1063 for (m = 0; (m < DIM); m++)
1071 for (m = 0; (m < DIM); m++)
1076 gmx_rmpbc(gpbc, natom, box, x);
1078 /* Begin loop of all molecules in frame */
1079 for (n = 0; (n < ncos); n++)
1081 for (i = 0; (i < gnx[n]); i++)
1085 ind0 = mols->index[molindex[n][i]];
1086 ind1 = mols->index[molindex[n][i]+1];
1088 mol_dip(ind0, ind1, x, atom, dipole[i]);
1089 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1090 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1093 update_slab_dipoles(ind0, ind1, x,
1094 dipole[i], idim, nslices, slab_dipoles, box);
1098 mol_quad(ind0, ind1, x, atom, quad);
1099 for (m = 0; (m < DIM); m++)
1101 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1104 if (bCorr && !bTotal)
1107 muall[i][tel3+XX] = dipole[i][XX];
1108 muall[i][tel3+YY] = dipole[i][YY];
1109 muall[i][tel3+ZZ] = dipole[i][ZZ];
1112 for (m = 0; (m < DIM); m++)
1114 M_av[m] += dipole[i][m]; /* M per frame */
1115 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1117 mu_mol = sqrt(mu_mol);
1119 mu_ave += mu_mol; /* calc. the average mu */
1121 /* Update the dipole distribution */
1122 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1130 rvec2sprvec(dipole[i], dipsp[i]);
1132 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1134 if (dipsp[i][ZZ] < 0.5 * M_PI)
1143 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1145 if (dipsp[i][ZZ] < 0.5 * M_PI)
1154 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1156 if (dipsp[i][ZZ] < 0.5 * M_PI)
1165 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1167 if (dipsp[i][ZZ] < 0.5 * M_PI)
1178 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1183 x[ind0][XX]+dipole[i][XX]/25,
1184 x[ind0][YY]+dipole[i][YY]/25,
1185 x[ind0][ZZ]+dipole[i][ZZ]/25,
1189 } /* End loop of all molecules in frame */
1193 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1194 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1195 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1196 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1197 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1198 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1202 /* Compute square of total dipole */
1203 for (m = 0; (m < DIM); m++)
1205 M_av2[m] = M_av[m]*M_av[m];
1210 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1211 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1212 sqr(dipaxis[YY]-0.5)+
1213 sqr(dipaxis[ZZ]-0.5));
1216 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1217 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1221 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1222 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1228 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1235 muall[0][tel3+XX] = M_av[XX];
1236 muall[0][tel3+YY] = M_av[YY];
1237 muall[0][tel3+ZZ] = M_av[ZZ];
1240 /* Write to file the total dipole moment of the box, and its components
1243 if ((skip == 0) || ((teller % skip) == 0))
1245 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1246 t, M_av[XX], M_av[YY], M_av[ZZ],
1247 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1250 for (m = 0; (m < DIM); m++)
1254 M4[m] += sqr(M_av2[m]);
1256 /* Increment loop counter */
1259 /* Calculate for output the running averages */
1260 invtel = 1.0/teller;
1261 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1262 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1263 M_diff = M2_ave - M_ave2;
1265 /* Compute volume from box in traj, else we use the one from above */
1272 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1274 /* Calculate running average for dipole */
1277 mu_aver = (mu_ave/gnx_tot)*invtel;
1280 if ((skip == 0) || ((teller % skip) == 0))
1282 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1283 * the two. Here M is sum mu_i. Further write the finite system
1284 * Kirkwood G factor and epsilon.
1286 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1287 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1292 gmx_stats_get_average(muframelsq, &aver);
1293 fprintf(adip, "%10g %f \n", t, aver);
1296 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1298 if (!bMU || (mu_aver != -1))
1300 /* Finite system Kirkwood G-factor */
1301 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1302 /* Infinite system Kirkwood G-factor */
1303 if (epsilonRF == 0.0)
1305 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1309 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1310 Gk/(3*epsilon*(2*epsilonRF+1)));
1313 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1318 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1321 gmx_stats_done(muframelsq);
1325 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1329 bCont = read_next_x(oenv, status, &t, x, box);
1331 timecheck = check_times(t);
1333 while (bCont && (timecheck == 0) );
1335 gmx_rmpbc_done(gpbc);
1342 gmx_ffclose(outmtot);
1343 gmx_ffclose(outaver);
1344 gmx_ffclose(outeps);
1358 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1359 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1360 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1361 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1362 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1368 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1369 sfree(slab_dipoles);
1373 printf("Average volume over run is %g\n", vol_aver);
1376 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1377 print_cmap(cmap, gkrbin, nlevels);
1379 /* Autocorrelation function */
1384 printf("Not enough frames for autocorrelation\n");
1388 dt = (t1 - t0)/(teller-1);
1389 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1395 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1396 teller, 1, muall, dt, mode, TRUE);
1400 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1401 teller, gnx_tot, muall, dt,
1402 mode, strcmp(corrtype, "molsep"));
1408 real aver, sigma, error;
1410 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1411 printf("\nDipole moment (Debye)\n");
1412 printf("---------------------\n");
1413 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1414 aver, sigma, error);
1418 for (m = 0; (m < DIM); m++)
1420 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1423 printf("\nQuadrupole moment (Debye-Ang)\n");
1424 printf("-----------------------------\n");
1425 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1426 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1427 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1431 printf("The following averages for the complete trajectory have been calculated:\n\n");
1432 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1433 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1434 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1436 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1437 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1438 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1440 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1441 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1443 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1445 if (!bMU || (mu_aver != -1))
1447 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1448 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1450 printf("Epsilon = %g\n", epsilon);
1454 /* Write to file the dipole moment distibution during the simulation.
1456 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1457 for (i = 0; (i < ndipbin); i++)
1459 fprintf(outdd, "%10g %10f\n",
1460 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1467 done_gkrbin(&gkrbin);
1471 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1480 while (m < mols->nr && index[i] != mols->index[m])
1486 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1488 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1490 if (i >= *n || index[i] != j)
1492 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1496 /* Modify the index in place */
1499 printf("There are %d molecules in the selection\n", nmol);
1503 int gmx_dipoles(int argc, char *argv[])
1505 const char *desc[] = {
1506 "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1507 "system. From this you can compute e.g. the dielectric constant for",
1508 "low-dielectric media.",
1509 "For molecules with a net charge, the net charge is subtracted at",
1510 "center of mass of the molecule.[PAR]",
1511 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1512 "components as well as the norm of the vector.",
1513 "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",
1515 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1517 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1518 "Furthermore, the dipole autocorrelation function will be computed when",
1519 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1521 "The correlation functions can be averaged over all molecules",
1522 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1523 "or it can be computed over the total dipole moment of the simulation box",
1524 "([TT]total[tt]).[PAR]",
1525 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1526 "G-factor, as well as the average cosine of the angle between the dipoles",
1527 "as a function of the distance. The plot also includes gOO and hOO",
1528 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1529 "we also include the energy per scale computed by taking the inner product of",
1530 "the dipoles divided by the distance to the third power.[PAR]",
1533 "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1534 "This will calculate the autocorrelation function of the molecular",
1535 "dipoles using a first order Legendre polynomial of the angle of the",
1536 "dipole vector and itself a time t later. For this calculation 1001",
1537 "frames will be used. Further, the dielectric constant will be calculated",
1538 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1539 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1540 "distribution function a maximum of 5.0 will be used."
1542 real mu_max = 5, mu_aver = -1, rcmax = 0;
1543 real epsilonRF = 0.0, temp = 300;
1544 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1545 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1546 const char *axtitle = "Z";
1547 int nslices = 10; /* nr of slices defined */
1548 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1549 int nlevels = 20, ndegrees = 90;
1552 { "-mu", FALSE, etREAL, {&mu_aver},
1553 "dipole of a single molecule (in Debye)" },
1554 { "-mumax", FALSE, etREAL, {&mu_max},
1555 "max dipole in Debye (for histogram)" },
1556 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1557 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1558 { "-skip", FALSE, etINT, {&skip},
1559 "Skip steps in the output (but not in the computations)" },
1560 { "-temp", FALSE, etREAL, {&temp},
1561 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1562 { "-corr", FALSE, etENUM, {corrtype},
1563 "Correlation function to calculate" },
1564 { "-pairs", FALSE, etBOOL, {&bPairs},
1565 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1566 { "-quad", FALSE, etBOOL, {&bQuad},
1567 "Take quadrupole into account"},
1568 { "-ncos", FALSE, etINT, {&ncos},
1569 "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." },
1570 { "-axis", FALSE, etSTR, {&axtitle},
1571 "Take the normal on the computational box in direction X, Y or Z." },
1572 { "-sl", FALSE, etINT, {&nslices},
1573 "Divide the box into this number of slices." },
1574 { "-gkratom", FALSE, etINT, {&nFA},
1575 "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" },
1576 { "-gkratom2", FALSE, etINT, {&nFB},
1577 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1578 { "-rcmax", FALSE, etREAL, {&rcmax},
1579 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1580 { "-phi", FALSE, etBOOL, {&bPhi},
1581 "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." },
1582 { "-nlevels", FALSE, etINT, {&nlevels},
1583 "Number of colors in the cmap output" },
1584 { "-ndegrees", FALSE, etINT, {&ndegrees},
1585 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1590 char **grpname = NULL;
1591 gmx_bool bGkr, bMU, bSlab;
1593 { efEDR, "-en", NULL, ffOPTRD },
1594 { efTRX, "-f", NULL, ffREAD },
1595 { efTPX, NULL, NULL, ffREAD },
1596 { efNDX, NULL, NULL, ffOPTRD },
1597 { efXVG, "-o", "Mtot", ffWRITE },
1598 { efXVG, "-eps", "epsilon", ffWRITE },
1599 { efXVG, "-a", "aver", ffWRITE },
1600 { efXVG, "-d", "dipdist", ffWRITE },
1601 { efXVG, "-c", "dipcorr", ffOPTWR },
1602 { efXVG, "-g", "gkr", ffOPTWR },
1603 { efXVG, "-adip", "adip", ffOPTWR },
1604 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1605 { efXVG, "-cos", "cosaver", ffOPTWR },
1606 { efXPM, "-cmap", "cmap", ffOPTWR },
1607 { efXVG, "-slab", "slab", ffOPTWR }
1609 #define NFILE asize(fnm)
1618 ppa = add_acf_pargs(&npargs, pa);
1619 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1620 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1625 printf("Using %g as mu_max and %g as the dipole moment.\n",
1627 if (epsilonRF == 0.0)
1629 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1632 bMU = opt2bSet("-en", NFILE, fnm);
1635 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.");
1637 bGkr = opt2bSet("-g", NFILE, fnm);
1638 if (opt2parg_bSet("-ncos", asize(pa), pa))
1640 if ((ncos != 1) && (ncos != 2))
1642 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1646 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1647 opt2parg_bSet("-axis", asize(pa), pa));
1652 printf("WARNING: Can not determine quadrupoles from energy file\n");
1657 printf("WARNING: Can not determine Gk(r) from energy file\n");
1663 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1664 " not enter a valid dipole for the molecules\n");
1669 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1670 &natoms, NULL, NULL, NULL, top);
1673 snew(grpname, ncos);
1674 snew(grpindex, ncos);
1675 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1676 ncos, gnx, grpindex, grpname);
1677 for (k = 0; (k < ncos); k++)
1679 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1680 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1684 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1685 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1686 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1687 opt2fn_null("-cos", NFILE, fnm),
1688 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1689 bPairs, corrtype[0],
1690 opt2fn("-c", NFILE, fnm),
1691 bGkr, opt2fn("-g", NFILE, fnm),
1692 bPhi, &nlevels, ndegrees,
1694 opt2fn("-cmap", NFILE, fnm), rcmax,
1695 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1696 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1697 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1699 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1700 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1701 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1702 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1703 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");