1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
53 #include "gmx_statistics.h"
66 #define e2d(x) ENM2DEBYE*(x)
67 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
68 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
80 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
88 if ((ptr = getenv("GKRWIDTH")) != NULL)
90 double bw = strtod(ptr, NULL);
95 gb->spacing = 0.01; /* nm */
97 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
104 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
107 snew(gb->elem, gb->nelem);
108 snew(gb->count, gb->nelem);
110 snew(gb->cmap, gb->nx);
111 gb->ny = std::max(2, ndegrees);
112 for (i = 0; (i < gb->nx); i++)
114 snew(gb->cmap[i], gb->ny);
121 static void done_gkrbin(t_gkrbin **gb)
129 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
131 int cy, index = gmx_nint(r/gb->spacing);
134 if (index < gb->nelem)
136 gb->elem[index] += cosa;
144 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
148 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
150 cy = std::min(gb->ny-1, std::max(0, cy));
153 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
155 gb->cmap[index][cy] += 1;
159 static void rvec2sprvec(rvec dipcart, rvec dipsp)
162 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
163 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
164 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
169 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
170 int mindex[], rvec x[], rvec mu[],
171 int ePBC, matrix box, t_atom *atom, int *nAtom)
173 static rvec *xcm[2] = { NULL, NULL};
174 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
175 real qtot, q, cosa, rr, phi;
179 for (n = 0; (n < ncos); n++)
183 snew(xcm[n], ngrp[n]);
185 for (i = 0; (i < ngrp[n]); i++)
187 /* Calculate center of mass of molecule */
193 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
198 clear_rvec(xcm[n][i]);
200 for (j = j0; j < j1; j++)
204 for (k = 0; k < DIM; k++)
206 xcm[n][i][k] += q*x[j][k];
209 svmul(1/qtot, xcm[n][i], xcm[n][i]);
213 set_pbc(&pbc, ePBC, box);
216 for (i = 0; i < ngrp[grp0]; i++)
218 gi = molindex[grp0][i];
219 j0 = (ncos == 2) ? 0 : i+1;
220 for (j = j0; j < ngrp[grp1]; j++)
222 gj = molindex[grp1][j];
223 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
225 /* Compute distance between molecules including PBC */
226 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
232 rvec r_ij, r_kj, r_kl, mm, nn;
236 copy_rvec(xcm[grp0][i], xj);
237 copy_rvec(xcm[grp1][j], xk);
238 rvec_add(xj, mu[gi], xi);
239 rvec_add(xk, mu[gj], xl);
240 phi = dih_angle(xi, xj, xk, xl, &pbc,
241 r_ij, r_kj, r_kl, mm, nn, /* out */
242 &sign, &t1, &t2, &t3);
247 cosa = cos_angle(mu[gi], mu[gj]);
250 if (debug || gmx_isnan(cosa))
252 fprintf(debug ? debug : stderr,
253 "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",
254 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
255 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
259 add2gkr(gb, rr, cosa, phi);
265 static real normalize_cmap(t_gkrbin *gb)
272 for (i = 0; (i < gb->nx); i++)
274 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
275 for (j = 0; (j < gb->ny); j++)
277 gb->cmap[i][j] /= vol;
278 hi = std::max(hi, gb->cmap[i][j]);
283 gmx_fatal(FARGS, "No data in the cmap");
288 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
295 t_rgb rlo = { 1, 1, 1 };
296 t_rgb rhi = { 0, 0, 0 };
298 hi = normalize_cmap(gb);
299 snew(xaxis, gb->nx+1);
300 for (i = 0; (i < gb->nx+1); i++)
302 xaxis[i] = i*gb->spacing;
305 for (j = 0; (j < gb->ny); j++)
309 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
313 yaxis[j] = (180.0*j)/(gb->ny-1.0);
315 /*2.0*j/(gb->ny-1.0)-1.0;*/
317 out = ffopen(cmap, "w");
319 "Dipole Orientation Distribution", "Fraction", "r (nm)",
320 gb->bPhi ? "Phi" : "Alpha",
321 gb->nx, gb->ny, xaxis, yaxis,
322 gb->cmap, 0, hi, rlo, rhi, nlevels);
328 static void print_gkrbin(const char *fn, t_gkrbin *gb,
329 int ngrp, int nframes, real volume,
330 const output_env_t oenv)
332 /* We compute Gk(r), gOO and hOO according to
333 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
334 * In this implementation the angle between dipoles is stored
335 * rather than their inner product. This allows to take polarizible
336 * models into account. The RDF is calculated as well, almost for free!
339 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
341 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
344 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
345 xvgr_legend(fp, asize(leg), leg, oenv);
347 Gkr = 1; /* Self-dipole inproduct = 1 */
352 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
353 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
357 while (last > 1 && gb->elem[last-1] == 0)
362 /* Divide by dipole squared, by number of frames, by number of origins.
363 * Multiply by 2 because we only take half the matrix of interactions
366 fac = 2.0/((double) ngrp * (double) nframes);
369 for (i = 0; i < last; i++)
371 /* Centre of the coordinate in the spherical layer */
374 /* Volume of the layer */
375 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
378 gOO = gb->count[i]*fac/(rho*vol_s);
380 /* Dipole correlation hOO, normalized by the relative number density, like
381 * in a Radial distribution function.
383 ggg = gb->elem[i]*fac;
384 hOO = 3.0*ggg/(rho*vol_s);
388 cosav = gb->elem[i]/gb->count[i];
394 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
396 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
397 x1, Gkr, cosav, hOO, gOO, ener);
405 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
406 real *t, int nre, t_enxframe *fr)
412 bCont = do_enx(fmu, fr);
415 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
416 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
421 if (Vol != -1) /* we've got Volume in the energy file */
423 *vol = fr->ener[Vol].e;
425 for (i = 0; i < DIM; i++)
427 mu[i] = fr->ener[iMu[i]].e;
435 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
438 int ncharged, m, a0, a1, a;
441 for (m = 0; m < n; m++)
443 a0 = mols->index[index[m]];
444 a1 = mols->index[index[m]+1];
447 for (a = a0; a < a1; a++)
452 /* This check is only for the count print */
453 if (fabs(qtot) > 0.01)
459 /* Remove the net charge at the center of mass */
460 for (a = a0; a < a1; a++)
462 atom[a].q -= qtot*atom[a].m/mtot;
469 printf("There are %d charged molecules in the selection,\n"
470 "will subtract their charge at their center of mass\n", ncharged);
474 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
480 for (k = k0; (k < k1); k++)
483 for (m = 0; (m < DIM); m++)
490 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
492 int i, k, m, n, niter;
493 real q, r2, mass, masstot;
494 rvec com; /* center of mass */
495 rvec r; /* distance of atoms to center of mass */
497 double dd[DIM], **ev, tmp;
501 for (i = 0; (i < DIM); i++)
508 /* Compute center of mass */
511 for (k = k0; (k < k1); k++)
515 for (i = 0; (i < DIM); i++)
517 com[i] += mass*x[k][i];
520 svmul((1.0/masstot), com, com);
522 /* We want traceless quadrupole moments, so let us calculate the complete
523 * quadrupole moment tensor and diagonalize this tensor to get
524 * the individual components on the diagonal.
527 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
529 for (m = 0; (m < DIM); m++)
531 for (n = 0; (n < DIM); n++)
536 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
538 q = (atom[k].q)*100.0;
539 rvec_sub(x[k], com, r);
541 for (m = 0; (m < DIM); m++)
543 for (n = 0; (n < DIM); n++)
545 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
551 for (i = 0; (i < DIM); i++)
553 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
554 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
558 /* We've got the quadrupole tensor, now diagonalize the sucker */
559 jacobi(inten, 3, dd, ev, &niter);
563 for (i = 0; (i < DIM); i++)
565 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
566 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
568 for (i = 0; (i < DIM); i++)
570 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
571 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
574 /* Sort the eigenvalues, for water we know that the order is as follows:
578 * At the moment I have no idea how this will work out for other molecules...
582 if (dd[i+1] > dd[i]) { \
591 for (m = 0; (m < DIM); m++)
593 quad[0] = dd[2]; /* yy */
594 quad[1] = dd[0]; /* zz */
595 quad[2] = dd[1]; /* xx */
600 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
604 for (i = 0; (i < DIM); i++)
614 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
616 real calc_eps(double M_diff, double volume, double epsRF, double temp)
618 double eps, A, teller, noemer;
619 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
620 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
622 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
631 teller = 1 + (A*2*epsRF/(2*epsRF+1));
632 noemer = 1 - (A/(2*epsRF+1));
634 eps = teller / noemer;
639 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
640 int idim, int nslice, rvec slab_dipole[],
646 for (k = k0; (k < k1); k++)
651 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
652 rvec_inc(slab_dipole[k], mu);
655 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
656 rvec slab_dipole[], matrix box, int nframes,
657 const output_env_t oenv)
663 const char *leg_dim[4] = {
664 "\\f{12}m\\f{4}\\sX\\N",
665 "\\f{12}m\\f{4}\\sY\\N",
666 "\\f{12}m\\f{4}\\sZ\\N",
667 "\\f{12}m\\f{4}\\stot\\N"
670 sprintf(buf, "Box-%c (nm)", 'X'+idim);
671 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
673 xvgr_legend(fp, DIM, leg_dim, oenv);
674 for (i = 0; (i < nslice); i++)
676 mutot = norm(slab_dipole[i])/nframes;
677 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
678 ((i+0.5)*box[idim][idim])/nslice,
679 slab_dipole[i][XX]/nframes,
680 slab_dipole[i][YY]/nframes,
681 slab_dipole[i][ZZ]/nframes,
685 do_view(oenv, fn, "-autoscale xy -nxy");
688 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
691 double dc, d, ddc1, ddc2, ddc3;
692 rvec xxx = { 1, 0, 0 };
693 rvec yyy = { 0, 1, 0 };
694 rvec zzz = { 0, 0, 1 };
697 ddc1 = ddc2 = ddc3 = 0;
698 for (i = k = 0; (i < n); i++)
700 ddc1 += fabs(cos_angle(dip[i], xxx));
701 ddc2 += fabs(cos_angle(dip[i], yyy));
702 ddc3 += fabs(cos_angle(dip[i], zzz));
705 for (j = i+1; (j < n); j++, k++)
707 dc = cos_angle(dip[i], dip[j]);
718 static void do_dip(t_topology *top, int ePBC, real volume,
720 const char *out_mtot, const char *out_eps,
721 const char *out_aver, const char *dipdist,
722 const char *cosaver, const char *fndip3d,
723 const char *fnadip, gmx_bool bPairs,
724 const char *corrtype, const char *corf,
725 gmx_bool bGkr, const char *gkrfn,
726 gmx_bool bPhi, int *nlevels, int ndegrees,
728 const char *cmap, real rcmax,
730 gmx_bool bMU, const char *mufn,
731 int *gnx, int *molindex[],
732 real mu_max, real mu_aver,
733 real epsilonRF, real temp,
734 int *gkatom, int skip,
735 gmx_bool bSlab, int nslices,
736 const char *axtitle, const char *slabfn,
737 const output_env_t oenv)
739 const char *leg_mtot[] = {
745 #define NLEGMTOT asize(leg_mtot)
746 const char *leg_eps[] = {
751 #define NLEGEPS asize(leg_eps)
752 const char *leg_aver[] = {
755 "< |M|\\S2\\N > - < |M| >\\S2\\N",
756 "< |M| >\\S2\\N / < |M|\\S2\\N >"
758 #define NLEGAVER asize(leg_aver)
759 const char *leg_cosaver[] = {
760 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
762 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
763 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
764 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
766 #define NLEGCOSAVER asize(leg_cosaver)
767 const char *leg_adip[] = {
772 #define NLEGADIP asize(leg_adip)
774 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
775 FILE *dip3d = NULL, *adip = NULL;
776 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
777 t_gkrbin *gkrbin = NULL;
778 gmx_enxnm_t *enm = NULL;
780 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
781 ener_file_t fmu = NULL;
782 int i, n, m, natom = 0, gnx_tot, teller, tel3;
784 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
786 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
789 gmx_bool bCorr, bTotal, bCont;
790 double M_diff = 0, epsilon, invtel, vol_aver;
791 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
792 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
793 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
796 rvec *slab_dipoles = NULL;
798 t_block *mols = NULL;
799 gmx_rmpbc_t gpbc = NULL;
811 /* initialize to a negative value so we can see that it wasn't picked up */
812 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
815 fmu = open_enx(mufn, "r");
816 do_enxnms(fmu, &nre, &enm);
818 /* Determine the indexes of the energy grps we need */
819 for (i = 0; (i < nre); i++)
821 if (strstr(enm[i].name, "Volume"))
825 else if (strstr(enm[i].name, "Mu-X"))
829 else if (strstr(enm[i].name, "Mu-Y"))
833 else if (strstr(enm[i].name, "Mu-Z"))
838 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
840 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
845 atom = top->atoms.atom;
848 if ((iVol == -1) && bMU)
850 printf("Using Volume from topology: %g nm^3\n", volume);
853 /* Correlation stuff */
854 bCorr = (corrtype[0] != 'n');
855 bTotal = (corrtype[0] == 't');
861 snew(muall[0], nframes*DIM);
866 for (i = 0; (i < gnx[0]); i++)
868 snew(muall[i], nframes*DIM);
873 /* Allocate array which contains for every molecule in a frame the
878 snew(dipole, gnx_tot);
883 for (i = 0; (i < DIM); i++)
885 Qlsq[i] = gmx_stats_init();
887 mulsq = gmx_stats_init();
889 /* Open all the files */
890 outmtot = xvgropen(out_mtot,
891 "Total dipole moment of the simulation box vs. time",
892 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
893 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
894 "Time (ps)", "", oenv);
895 outaver = xvgropen(out_aver, "Total dipole moment",
896 "Time (ps)", "D", oenv);
899 idim = axtitle[0] - 'X';
900 if ((idim < 0) || (idim >= DIM))
902 idim = axtitle[0] - 'x';
904 if ((idim < 0) || (idim >= DIM))
912 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
913 axtitle, nslices, idim);
916 snew(slab_dipoles, nslices);
917 fprintf(stderr, "Doing slab analysis\n");
923 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
924 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
929 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
930 "Average absolute dipole orientation", "Time (ps)", "", oenv);
931 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
937 snew(dipsp, gnx_tot);
939 /* we need a dummy file for gnuplot */
940 dip3d = (FILE *)ffopen("dummy.dat", "w");
941 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
944 dip3d = (FILE *)ffopen(fndip3d, "w");
945 fprintf(dip3d, "# This file was created by %s\n", Program());
946 fprintf(dip3d, "# which is part of G R O M A C S:\n");
947 fprintf(dip3d, "#\n");
950 /* Write legends to all the files */
951 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
952 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
954 if (bMU && (mu_aver == -1))
956 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
960 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
966 /* Read the first frame from energy or traj file */
971 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
974 timecheck = check_times(t);
979 if ((teller % 10) == 0)
981 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
986 printf("End of %s reached\n", mufn);
990 while (bCont && (timecheck < 0));
994 natom = read_first_x(oenv, &status, fn, &t, &x, box);
997 /* Calculate spacing for dipole bin (simple histogram) */
998 ndipbin = 1 + static_cast<int>(mu_max/0.01);
999 snew(dipole_bin, ndipbin);
1001 for (m = 0; (m < DIM); m++)
1003 M[m] = M2[m] = M4[m] = 0.0;
1008 /* Use 0.7 iso 0.5 to account for pressure scaling */
1009 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1010 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1012 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1014 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1016 /* Start while loop over frames */
1021 if (bCorr && (teller >= nframes))
1026 srenew(muall[0], nframes*DIM);
1030 for (i = 0; (i < gnx_tot); i++)
1032 srenew(muall[i], nframes*DIM);
1038 muframelsq = gmx_stats_init();
1041 for (m = 0; (m < DIM); m++)
1048 /* Copy rvec into double precision local variable */
1049 for (m = 0; (m < DIM); m++)
1057 for (m = 0; (m < DIM); m++)
1062 gmx_rmpbc(gpbc, natom, box, x);
1064 /* Begin loop of all molecules in frame */
1065 for (n = 0; (n < ncos); n++)
1067 for (i = 0; (i < gnx[n]); i++)
1071 ind0 = mols->index[molindex[n][i]];
1072 ind1 = mols->index[molindex[n][i]+1];
1074 mol_dip(ind0, ind1, x, atom, dipole[i]);
1075 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1076 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1079 update_slab_dipoles(ind0, ind1, x,
1080 dipole[i], idim, nslices, slab_dipoles, box);
1084 mol_quad(ind0, ind1, x, atom, quad);
1085 for (m = 0; (m < DIM); m++)
1087 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1090 if (bCorr && !bTotal)
1093 muall[i][tel3+XX] = dipole[i][XX];
1094 muall[i][tel3+YY] = dipole[i][YY];
1095 muall[i][tel3+ZZ] = dipole[i][ZZ];
1098 for (m = 0; (m < DIM); m++)
1100 M_av[m] += dipole[i][m]; /* M per frame */
1101 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1103 mu_mol = sqrt(mu_mol);
1105 mu_ave += mu_mol; /* calc. the average mu */
1107 /* Update the dipole distribution */
1108 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1116 rvec2sprvec(dipole[i], dipsp[i]);
1118 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1120 if (dipsp[i][ZZ] < 0.5 * M_PI)
1129 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1131 if (dipsp[i][ZZ] < 0.5 * M_PI)
1140 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1142 if (dipsp[i][ZZ] < 0.5 * M_PI)
1151 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1153 if (dipsp[i][ZZ] < 0.5 * M_PI)
1164 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1169 x[ind0][XX]+dipole[i][XX]/25,
1170 x[ind0][YY]+dipole[i][YY]/25,
1171 x[ind0][ZZ]+dipole[i][ZZ]/25,
1175 } /* End loop of all molecules in frame */
1179 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1180 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1181 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1182 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1183 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1184 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1188 /* Compute square of total dipole */
1189 for (m = 0; (m < DIM); m++)
1191 M_av2[m] = M_av[m]*M_av[m];
1196 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1197 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1198 sqr(dipaxis[YY]-0.5)+
1199 sqr(dipaxis[ZZ]-0.5));
1202 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1203 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1207 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1208 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1214 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1221 muall[0][tel3+XX] = M_av[XX];
1222 muall[0][tel3+YY] = M_av[YY];
1223 muall[0][tel3+ZZ] = M_av[ZZ];
1226 /* Write to file the total dipole moment of the box, and its components
1229 if ((skip == 0) || ((teller % skip) == 0))
1231 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1232 t, M_av[XX], M_av[YY], M_av[ZZ],
1233 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1236 for (m = 0; (m < DIM); m++)
1240 M4[m] += sqr(M_av2[m]);
1242 /* Increment loop counter */
1245 /* Calculate for output the running averages */
1246 invtel = 1.0/teller;
1247 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1248 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1249 M_diff = M2_ave - M_ave2;
1251 /* Compute volume from box in traj, else we use the one from above */
1258 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1260 /* Calculate running average for dipole */
1263 mu_aver = (mu_ave/gnx_tot)*invtel;
1266 if ((skip == 0) || ((teller % skip) == 0))
1268 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1269 * the two. Here M is sum mu_i. Further write the finite system
1270 * Kirkwood G factor and epsilon.
1272 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1273 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1278 gmx_stats_get_average(muframelsq, &aver);
1279 fprintf(adip, "%10g %f \n", t, aver);
1282 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1284 if (!bMU || (mu_aver != -1))
1286 /* Finite system Kirkwood G-factor */
1287 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1288 /* Infinite system Kirkwood G-factor */
1289 if (epsilonRF == 0.0)
1291 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1295 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1296 Gk/(3*epsilon*(2*epsilonRF+1)));
1299 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1304 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1307 gmx_stats_done(muframelsq);
1311 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1315 bCont = read_next_x(oenv, status, &t, x, box);
1317 timecheck = check_times(t);
1319 while (bCont && (timecheck == 0) );
1321 gmx_rmpbc_done(gpbc);
1344 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1345 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1346 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1347 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1348 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1354 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1355 sfree(slab_dipoles);
1359 printf("Average volume over run is %g\n", vol_aver);
1362 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1363 print_cmap(cmap, gkrbin, nlevels);
1365 /* Autocorrelation function */
1370 printf("Not enough frames for autocorrelation\n");
1374 dt = (t1 - t0)/(teller-1);
1375 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1381 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1382 teller, 1, muall, dt, mode, TRUE);
1386 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1387 teller, gnx_tot, muall, dt,
1388 mode, strcmp(corrtype, "molsep"));
1394 real aver, sigma, error;
1396 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1397 printf("\nDipole moment (Debye)\n");
1398 printf("---------------------\n");
1399 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1400 aver, sigma, error);
1404 for (m = 0; (m < DIM); m++)
1406 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1409 printf("\nQuadrupole moment (Debye-Ang)\n");
1410 printf("-----------------------------\n");
1411 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1412 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1413 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1417 printf("The following averages for the complete trajectory have been calculated:\n\n");
1418 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1419 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1420 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1422 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1423 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1424 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1426 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1427 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1429 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1431 if (!bMU || (mu_aver != -1))
1433 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1434 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1436 printf("Epsilon = %g\n", epsilon);
1440 /* Write to file the dipole moment distibution during the simulation.
1442 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1443 for (i = 0; (i < ndipbin); i++)
1445 fprintf(outdd, "%10g %10f\n",
1446 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1453 done_gkrbin(&gkrbin);
1457 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1466 while (m < mols->nr && index[i] != mols->index[m])
1472 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1474 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1476 if (i >= *n || index[i] != j)
1478 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1482 /* Modify the index in place */
1485 printf("There are %d molecules in the selection\n", nmol);
1489 int gmx_dipoles(int argc, char *argv[])
1491 const char *desc[] = {
1492 "[TT]g_dipoles[tt] computes the total dipole plus fluctuations of a simulation",
1493 "system. From this you can compute e.g. the dielectric constant for",
1494 "low-dielectric media.",
1495 "For molecules with a net charge, the net charge is subtracted at",
1496 "center of mass of the molecule.[PAR]",
1497 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1498 "components as well as the norm of the vector.",
1499 "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",
1501 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1503 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1504 "Furthermore, the dipole autocorrelation function will be computed when",
1505 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1507 "The correlation functions can be averaged over all molecules",
1508 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1509 "or it can be computed over the total dipole moment of the simulation box",
1510 "([TT]total[tt]).[PAR]",
1511 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1512 "G-factor, as well as the average cosine of the angle between the dipoles",
1513 "as a function of the distance. The plot also includes gOO and hOO",
1514 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1515 "we also include the energy per scale computed by taking the inner product of",
1516 "the dipoles divided by the distance to the third power.[PAR]",
1519 "[TT]g_dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1520 "This will calculate the autocorrelation function of the molecular",
1521 "dipoles using a first order Legendre polynomial of the angle of the",
1522 "dipole vector and itself a time t later. For this calculation 1001",
1523 "frames will be used. Further, the dielectric constant will be calculated",
1524 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1525 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1526 "distribution function a maximum of 5.0 will be used."
1528 real mu_max = 5, mu_aver = -1, rcmax = 0;
1529 real epsilonRF = 0.0, temp = 300;
1530 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1531 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1532 const char *axtitle = "Z";
1533 int nslices = 10; /* nr of slices defined */
1534 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1535 int nlevels = 20, ndegrees = 90;
1538 { "-mu", FALSE, etREAL, {&mu_aver},
1539 "dipole of a single molecule (in Debye)" },
1540 { "-mumax", FALSE, etREAL, {&mu_max},
1541 "max dipole in Debye (for histogram)" },
1542 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1543 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1544 { "-skip", FALSE, etINT, {&skip},
1545 "Skip steps in the output (but not in the computations)" },
1546 { "-temp", FALSE, etREAL, {&temp},
1547 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1548 { "-corr", FALSE, etENUM, {corrtype},
1549 "Correlation function to calculate" },
1550 { "-pairs", FALSE, etBOOL, {&bPairs},
1551 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1552 { "-quad", FALSE, etBOOL, {&bQuad},
1553 "Take quadrupole into account"},
1554 { "-ncos", FALSE, etINT, {&ncos},
1555 "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." },
1556 { "-axis", FALSE, etSTR, {&axtitle},
1557 "Take the normal on the computational box in direction X, Y or Z." },
1558 { "-sl", FALSE, etINT, {&nslices},
1559 "Divide the box into this number of slices." },
1560 { "-gkratom", FALSE, etINT, {&nFA},
1561 "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" },
1562 { "-gkratom2", FALSE, etINT, {&nFB},
1563 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1564 { "-rcmax", FALSE, etREAL, {&rcmax},
1565 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1566 { "-phi", FALSE, etBOOL, {&bPhi},
1567 "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." },
1568 { "-nlevels", FALSE, etINT, {&nlevels},
1569 "Number of colors in the cmap output" },
1570 { "-ndegrees", FALSE, etINT, {&ndegrees},
1571 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1576 char **grpname = NULL;
1577 gmx_bool bGkr, bMU, bSlab;
1579 { efEDR, "-en", NULL, ffOPTRD },
1580 { efTRX, "-f", NULL, ffREAD },
1581 { efTPX, NULL, NULL, ffREAD },
1582 { efNDX, NULL, NULL, ffOPTRD },
1583 { efXVG, "-o", "Mtot", ffWRITE },
1584 { efXVG, "-eps", "epsilon", ffWRITE },
1585 { efXVG, "-a", "aver", ffWRITE },
1586 { efXVG, "-d", "dipdist", ffWRITE },
1587 { efXVG, "-c", "dipcorr", ffOPTWR },
1588 { efXVG, "-g", "gkr", ffOPTWR },
1589 { efXVG, "-adip", "adip", ffOPTWR },
1590 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1591 { efXVG, "-cos", "cosaver", ffOPTWR },
1592 { efXPM, "-cmap", "cmap", ffOPTWR },
1593 { efXVG, "-slab", "slab", ffOPTWR }
1595 #define NFILE asize(fnm)
1604 ppa = add_acf_pargs(&npargs, pa);
1605 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1606 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv);
1608 printf("Using %g as mu_max and %g as the dipole moment.\n",
1610 if (epsilonRF == 0.0)
1612 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1615 bMU = opt2bSet("-en", NFILE, fnm);
1618 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.");
1620 bGkr = opt2bSet("-g", NFILE, fnm);
1621 if (opt2parg_bSet("-ncos", asize(pa), pa))
1623 if ((ncos != 1) && (ncos != 2))
1625 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1629 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1630 opt2parg_bSet("-axis", asize(pa), pa));
1635 printf("WARNING: Can not determine quadrupoles from energy file\n");
1640 printf("WARNING: Can not determine Gk(r) from energy file\n");
1646 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1647 " not enter a valid dipole for the molecules\n");
1652 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1653 &natoms, NULL, NULL, NULL, top);
1656 snew(grpname, ncos);
1657 snew(grpindex, ncos);
1658 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1659 ncos, gnx, grpindex, grpname);
1660 for (k = 0; (k < ncos); k++)
1662 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1663 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1667 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1668 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1669 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1670 opt2fn_null("-cos", NFILE, fnm),
1671 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1672 bPairs, corrtype[0],
1673 opt2fn("-c", NFILE, fnm),
1674 bGkr, opt2fn("-g", NFILE, fnm),
1675 bPhi, &nlevels, ndegrees,
1677 opt2fn("-cmap", NFILE, fnm), rcmax,
1678 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1679 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1680 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1682 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1683 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1684 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1685 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1686 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");