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
50 #include "gromacs/fileio/futil.h"
53 #include "gmx_statistics.h"
60 #include "gromacs/fileio/enxio.h"
62 #include "gromacs/fileio/matio.h"
65 #include "gromacs/fileio/trxio.h"
67 #include "gromacs/utility/exceptions.h"
68 #include "gromacs/utility/programinfo.h"
70 #define e2d(x) ENM2DEBYE*(x)
71 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
72 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
84 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
92 if ((ptr = getenv("GKRWIDTH")) != NULL)
94 double bw = strtod(ptr, NULL);
99 gb->spacing = 0.01; /* nm */
101 gb->nelem = 1 + static_cast<int>(radius/gb->spacing);
108 gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
111 snew(gb->elem, gb->nelem);
112 snew(gb->count, gb->nelem);
114 snew(gb->cmap, gb->nx);
115 gb->ny = std::max(2, ndegrees);
116 for (i = 0; (i < gb->nx); i++)
118 snew(gb->cmap[i], gb->ny);
125 static void done_gkrbin(t_gkrbin **gb)
133 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
135 int cy, index = gmx_nint(r/gb->spacing);
138 if (index < gb->nelem)
140 gb->elem[index] += cosa;
148 cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
152 cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
154 cy = std::min(gb->ny-1, std::max(0, cy));
157 fprintf(debug, "CY: %10f %5d\n", alpha, cy);
159 gb->cmap[index][cy] += 1;
163 static void rvec2sprvec(rvec dipcart, rvec dipsp)
166 dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
167 dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
168 dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
173 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
174 int mindex[], rvec x[], rvec mu[],
175 int ePBC, matrix box, t_atom *atom, int *nAtom)
177 static rvec *xcm[2] = { NULL, NULL};
178 int gi, gj, j0, j1, i, j, k, n, grp0, grp1;
179 real qtot, q, cosa, rr, phi;
183 for (n = 0; (n < ncos); n++)
187 snew(xcm[n], ngrp[n]);
189 for (i = 0; (i < ngrp[n]); i++)
191 /* Calculate center of mass of molecule */
197 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
202 clear_rvec(xcm[n][i]);
204 for (j = j0; j < j1; j++)
208 for (k = 0; k < DIM; k++)
210 xcm[n][i][k] += q*x[j][k];
213 svmul(1/qtot, xcm[n][i], xcm[n][i]);
217 set_pbc(&pbc, ePBC, box);
220 for (i = 0; i < ngrp[grp0]; i++)
222 gi = molindex[grp0][i];
223 j0 = (ncos == 2) ? 0 : i+1;
224 for (j = j0; j < ngrp[grp1]; j++)
226 gj = molindex[grp1][j];
227 if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
229 /* Compute distance between molecules including PBC */
230 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
236 rvec r_ij, r_kj, r_kl, mm, nn;
240 copy_rvec(xcm[grp0][i], xj);
241 copy_rvec(xcm[grp1][j], xk);
242 rvec_add(xj, mu[gi], xi);
243 rvec_add(xk, mu[gj], xl);
244 phi = dih_angle(xi, xj, xk, xl, &pbc,
245 r_ij, r_kj, r_kl, mm, nn, /* out */
246 &sign, &t1, &t2, &t3);
251 cosa = cos_angle(mu[gi], mu[gj]);
254 if (debug || gmx_isnan(cosa))
256 fprintf(debug ? debug : stderr,
257 "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",
258 gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
259 gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
263 add2gkr(gb, rr, cosa, phi);
269 static real normalize_cmap(t_gkrbin *gb)
276 for (i = 0; (i < gb->nx); i++)
278 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
279 for (j = 0; (j < gb->ny); j++)
281 gb->cmap[i][j] /= vol;
282 hi = std::max(hi, gb->cmap[i][j]);
287 gmx_fatal(FARGS, "No data in the cmap");
292 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
299 t_rgb rlo = { 1, 1, 1 };
300 t_rgb rhi = { 0, 0, 0 };
302 hi = normalize_cmap(gb);
303 snew(xaxis, gb->nx+1);
304 for (i = 0; (i < gb->nx+1); i++)
306 xaxis[i] = i*gb->spacing;
309 for (j = 0; (j < gb->ny); j++)
313 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
317 yaxis[j] = (180.0*j)/(gb->ny-1.0);
319 /*2.0*j/(gb->ny-1.0)-1.0;*/
321 out = ffopen(cmap, "w");
323 "Dipole Orientation Distribution", "Fraction", "r (nm)",
324 gb->bPhi ? "Phi" : "Alpha",
325 gb->nx, gb->ny, xaxis, yaxis,
326 gb->cmap, 0, hi, rlo, rhi, nlevels);
332 static void print_gkrbin(const char *fn, t_gkrbin *gb,
333 int ngrp, int nframes, real volume,
334 const output_env_t oenv)
336 /* We compute Gk(r), gOO and hOO according to
337 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
338 * In this implementation the angle between dipoles is stored
339 * rather than their inner product. This allows to take polarizible
340 * models into account. The RDF is calculated as well, almost for free!
343 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
345 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
348 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
349 xvgr_legend(fp, asize(leg), leg, oenv);
351 Gkr = 1; /* Self-dipole inproduct = 1 */
356 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
357 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
361 while (last > 1 && gb->elem[last-1] == 0)
366 /* Divide by dipole squared, by number of frames, by number of origins.
367 * Multiply by 2 because we only take half the matrix of interactions
370 fac = 2.0/((double) ngrp * (double) nframes);
373 for (i = 0; i < last; i++)
375 /* Centre of the coordinate in the spherical layer */
378 /* Volume of the layer */
379 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
382 gOO = gb->count[i]*fac/(rho*vol_s);
384 /* Dipole correlation hOO, normalized by the relative number density, like
385 * in a Radial distribution function.
387 ggg = gb->elem[i]*fac;
388 hOO = 3.0*ggg/(rho*vol_s);
392 cosav = gb->elem[i]/gb->count[i];
398 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
400 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
401 x1, Gkr, cosav, hOO, gOO, ener);
409 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
410 real *t, int nre, t_enxframe *fr)
416 bCont = do_enx(fmu, fr);
419 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
420 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
425 if (Vol != -1) /* we've got Volume in the energy file */
427 *vol = fr->ener[Vol].e;
429 for (i = 0; i < DIM; i++)
431 mu[i] = fr->ener[iMu[i]].e;
439 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
442 int ncharged, m, a0, a1, a;
445 for (m = 0; m < n; m++)
447 a0 = mols->index[index[m]];
448 a1 = mols->index[index[m]+1];
451 for (a = a0; a < a1; a++)
456 /* This check is only for the count print */
457 if (fabs(qtot) > 0.01)
463 /* Remove the net charge at the center of mass */
464 for (a = a0; a < a1; a++)
466 atom[a].q -= qtot*atom[a].m/mtot;
473 printf("There are %d charged molecules in the selection,\n"
474 "will subtract their charge at their center of mass\n", ncharged);
478 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
484 for (k = k0; (k < k1); k++)
487 for (m = 0; (m < DIM); m++)
494 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
496 int i, k, m, n, niter;
497 real q, r2, mass, masstot;
498 rvec com; /* center of mass */
499 rvec r; /* distance of atoms to center of mass */
501 double dd[DIM], **ev, tmp;
505 for (i = 0; (i < DIM); i++)
512 /* Compute center of mass */
515 for (k = k0; (k < k1); k++)
519 for (i = 0; (i < DIM); i++)
521 com[i] += mass*x[k][i];
524 svmul((1.0/masstot), com, com);
526 /* We want traceless quadrupole moments, so let us calculate the complete
527 * quadrupole moment tensor and diagonalize this tensor to get
528 * the individual components on the diagonal.
531 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
533 for (m = 0; (m < DIM); m++)
535 for (n = 0; (n < DIM); n++)
540 for (k = k0; (k < k1); k++) /* loop over atoms in a molecule */
542 q = (atom[k].q)*100.0;
543 rvec_sub(x[k], com, r);
545 for (m = 0; (m < DIM); m++)
547 for (n = 0; (n < DIM); n++)
549 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
555 for (i = 0; (i < DIM); i++)
557 fprintf(debug, "Q[%d] = %8.3f %8.3f %8.3f\n",
558 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
562 /* We've got the quadrupole tensor, now diagonalize the sucker */
563 jacobi(inten, 3, dd, ev, &niter);
567 for (i = 0; (i < DIM); i++)
569 fprintf(debug, "ev[%d] = %8.3f %8.3f %8.3f\n",
570 i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
572 for (i = 0; (i < DIM); i++)
574 fprintf(debug, "Q'[%d] = %8.3f %8.3f %8.3f\n",
575 i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
578 /* Sort the eigenvalues, for water we know that the order is as follows:
582 * At the moment I have no idea how this will work out for other molecules...
586 if (dd[i+1] > dd[i]) { \
595 for (m = 0; (m < DIM); m++)
597 quad[0] = dd[2]; /* yy */
598 quad[1] = dd[0]; /* zz */
599 quad[2] = dd[1]; /* xx */
604 pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
608 for (i = 0; (i < DIM); i++)
618 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
620 real calc_eps(double M_diff, double volume, double epsRF, double temp)
622 double eps, A, teller, noemer;
623 double eps_0 = 8.854187817e-12; /* epsilon_0 in C^2 J^-1 m^-1 */
624 double fac = 1.112650021e-59; /* converts Debye^2 to C^2 m^2 */
626 A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
635 teller = 1 + (A*2*epsRF/(2*epsRF+1));
636 noemer = 1 - (A/(2*epsRF+1));
638 eps = teller / noemer;
643 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
644 int idim, int nslice, rvec slab_dipole[],
650 for (k = k0; (k < k1); k++)
655 k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
656 rvec_inc(slab_dipole[k], mu);
659 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
660 rvec slab_dipole[], matrix box, int nframes,
661 const output_env_t oenv)
667 const char *leg_dim[4] = {
668 "\\f{12}m\\f{4}\\sX\\N",
669 "\\f{12}m\\f{4}\\sY\\N",
670 "\\f{12}m\\f{4}\\sZ\\N",
671 "\\f{12}m\\f{4}\\stot\\N"
674 sprintf(buf, "Box-%c (nm)", 'X'+idim);
675 fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
677 xvgr_legend(fp, DIM, leg_dim, oenv);
678 for (i = 0; (i < nslice); i++)
680 mutot = norm(slab_dipole[i])/nframes;
681 fprintf(fp, "%10.3f %10.3f %10.3f %10.3f %10.3f\n",
682 ((i+0.5)*box[idim][idim])/nslice,
683 slab_dipole[i][XX]/nframes,
684 slab_dipole[i][YY]/nframes,
685 slab_dipole[i][ZZ]/nframes,
689 do_view(oenv, fn, "-autoscale xy -nxy");
692 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
695 double dc, d, ddc1, ddc2, ddc3;
696 rvec xxx = { 1, 0, 0 };
697 rvec yyy = { 0, 1, 0 };
698 rvec zzz = { 0, 0, 1 };
701 ddc1 = ddc2 = ddc3 = 0;
702 for (i = k = 0; (i < n); i++)
704 ddc1 += fabs(cos_angle(dip[i], xxx));
705 ddc2 += fabs(cos_angle(dip[i], yyy));
706 ddc3 += fabs(cos_angle(dip[i], zzz));
709 for (j = i+1; (j < n); j++, k++)
711 dc = cos_angle(dip[i], dip[j]);
722 static void do_dip(t_topology *top, int ePBC, real volume,
724 const char *out_mtot, const char *out_eps,
725 const char *out_aver, const char *dipdist,
726 const char *cosaver, const char *fndip3d,
727 const char *fnadip, gmx_bool bPairs,
728 const char *corrtype, const char *corf,
729 gmx_bool bGkr, const char *gkrfn,
730 gmx_bool bPhi, int *nlevels, int ndegrees,
732 const char *cmap, real rcmax,
734 gmx_bool bMU, const char *mufn,
735 int *gnx, int *molindex[],
736 real mu_max, real mu_aver,
737 real epsilonRF, real temp,
738 int *gkatom, int skip,
739 gmx_bool bSlab, int nslices,
740 const char *axtitle, const char *slabfn,
741 const output_env_t oenv)
743 const char *leg_mtot[] = {
749 #define NLEGMTOT asize(leg_mtot)
750 const char *leg_eps[] = {
755 #define NLEGEPS asize(leg_eps)
756 const char *leg_aver[] = {
759 "< |M|\\S2\\N > - < |M| >\\S2\\N",
760 "< |M| >\\S2\\N / < |M|\\S2\\N >"
762 #define NLEGAVER asize(leg_aver)
763 const char *leg_cosaver[] = {
764 "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
766 "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
767 "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
768 "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
770 #define NLEGCOSAVER asize(leg_cosaver)
771 const char *leg_adip[] = {
776 #define NLEGADIP asize(leg_adip)
778 FILE *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
779 FILE *dip3d = NULL, *adip = NULL;
780 rvec *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
781 t_gkrbin *gkrbin = NULL;
782 gmx_enxnm_t *enm = NULL;
784 int nframes = 1000, nre, timecheck = 0, ncolour = 0;
785 ener_file_t fmu = NULL;
786 int i, n, m, natom = 0, gnx_tot, teller, tel3;
788 int *dipole_bin, ndipbin, ibin, iVol, idim = -1;
790 real rcut = 0, t, t0, t1, dt, dd, rms_cos;
793 gmx_bool bCorr, bTotal, bCont;
794 double M_diff = 0, epsilon, invtel, vol_aver;
795 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
796 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
797 gmx_stats_t *Qlsq, mulsq, muframelsq = NULL;
800 rvec *slab_dipoles = NULL;
802 t_block *mols = NULL;
803 gmx_rmpbc_t gpbc = NULL;
815 /* initialize to a negative value so we can see that it wasn't picked up */
816 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
819 fmu = open_enx(mufn, "r");
820 do_enxnms(fmu, &nre, &enm);
822 /* Determine the indexes of the energy grps we need */
823 for (i = 0; (i < nre); i++)
825 if (strstr(enm[i].name, "Volume"))
829 else if (strstr(enm[i].name, "Mu-X"))
833 else if (strstr(enm[i].name, "Mu-Y"))
837 else if (strstr(enm[i].name, "Mu-Z"))
842 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
844 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
849 atom = top->atoms.atom;
852 if ((iVol == -1) && bMU)
854 printf("Using Volume from topology: %g nm^3\n", volume);
857 /* Correlation stuff */
858 bCorr = (corrtype[0] != 'n');
859 bTotal = (corrtype[0] == 't');
865 snew(muall[0], nframes*DIM);
870 for (i = 0; (i < gnx[0]); i++)
872 snew(muall[i], nframes*DIM);
877 /* Allocate array which contains for every molecule in a frame the
882 snew(dipole, gnx_tot);
887 for (i = 0; (i < DIM); i++)
889 Qlsq[i] = gmx_stats_init();
891 mulsq = gmx_stats_init();
893 /* Open all the files */
894 outmtot = xvgropen(out_mtot,
895 "Total dipole moment of the simulation box vs. time",
896 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
897 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
898 "Time (ps)", "", oenv);
899 outaver = xvgropen(out_aver, "Total dipole moment",
900 "Time (ps)", "D", oenv);
903 idim = axtitle[0] - 'X';
904 if ((idim < 0) || (idim >= DIM))
906 idim = axtitle[0] - 'x';
908 if ((idim < 0) || (idim >= DIM))
916 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
917 axtitle, nslices, idim);
920 snew(slab_dipoles, nslices);
921 fprintf(stderr, "Doing slab analysis\n");
927 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
928 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
933 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
934 "Average absolute dipole orientation", "Time (ps)", "", oenv);
935 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
941 snew(dipsp, gnx_tot);
943 /* we need a dummy file for gnuplot */
944 dip3d = (FILE *)ffopen("dummy.dat", "w");
945 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
948 dip3d = (FILE *)ffopen(fndip3d, "w");
951 gmx::BinaryInformationSettings settings;
952 settings.generatedByHeader(true);
953 settings.linePrefix("# ");
954 gmx::printBinaryInformation(dip3d, gmx::ProgramInfo::getInstance(),
957 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
960 /* Write legends to all the files */
961 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
962 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
964 if (bMU && (mu_aver == -1))
966 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
970 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
976 /* Read the first frame from energy or traj file */
981 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
984 timecheck = check_times(t);
989 if ((teller % 10) == 0)
991 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
996 printf("End of %s reached\n", mufn);
1000 while (bCont && (timecheck < 0));
1004 natom = read_first_x(oenv, &status, fn, &t, &x, box);
1007 /* Calculate spacing for dipole bin (simple histogram) */
1008 ndipbin = 1 + static_cast<int>(mu_max/0.01);
1009 snew(dipole_bin, ndipbin);
1011 for (m = 0; (m < DIM); m++)
1013 M[m] = M2[m] = M4[m] = 0.0;
1018 /* Use 0.7 iso 0.5 to account for pressure scaling */
1019 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1020 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1022 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1024 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1026 /* Start while loop over frames */
1031 if (bCorr && (teller >= nframes))
1036 srenew(muall[0], nframes*DIM);
1040 for (i = 0; (i < gnx_tot); i++)
1042 srenew(muall[i], nframes*DIM);
1048 muframelsq = gmx_stats_init();
1051 for (m = 0; (m < DIM); m++)
1058 /* Copy rvec into double precision local variable */
1059 for (m = 0; (m < DIM); m++)
1067 for (m = 0; (m < DIM); m++)
1072 gmx_rmpbc(gpbc, natom, box, x);
1074 /* Begin loop of all molecules in frame */
1075 for (n = 0; (n < ncos); n++)
1077 for (i = 0; (i < gnx[n]); i++)
1081 ind0 = mols->index[molindex[n][i]];
1082 ind1 = mols->index[molindex[n][i]+1];
1084 mol_dip(ind0, ind1, x, atom, dipole[i]);
1085 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1086 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1089 update_slab_dipoles(ind0, ind1, x,
1090 dipole[i], idim, nslices, slab_dipoles, box);
1094 mol_quad(ind0, ind1, x, atom, quad);
1095 for (m = 0; (m < DIM); m++)
1097 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1100 if (bCorr && !bTotal)
1103 muall[i][tel3+XX] = dipole[i][XX];
1104 muall[i][tel3+YY] = dipole[i][YY];
1105 muall[i][tel3+ZZ] = dipole[i][ZZ];
1108 for (m = 0; (m < DIM); m++)
1110 M_av[m] += dipole[i][m]; /* M per frame */
1111 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1113 mu_mol = sqrt(mu_mol);
1115 mu_ave += mu_mol; /* calc. the average mu */
1117 /* Update the dipole distribution */
1118 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1126 rvec2sprvec(dipole[i], dipsp[i]);
1128 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1130 if (dipsp[i][ZZ] < 0.5 * M_PI)
1139 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1141 if (dipsp[i][ZZ] < 0.5 * M_PI)
1150 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1152 if (dipsp[i][ZZ] < 0.5 * M_PI)
1161 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1163 if (dipsp[i][ZZ] < 0.5 * M_PI)
1174 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1179 x[ind0][XX]+dipole[i][XX]/25,
1180 x[ind0][YY]+dipole[i][YY]/25,
1181 x[ind0][ZZ]+dipole[i][ZZ]/25,
1185 } /* End loop of all molecules in frame */
1189 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1190 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1191 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1192 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1193 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1194 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1198 /* Compute square of total dipole */
1199 for (m = 0; (m < DIM); m++)
1201 M_av2[m] = M_av[m]*M_av[m];
1206 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1207 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1208 sqr(dipaxis[YY]-0.5)+
1209 sqr(dipaxis[ZZ]-0.5));
1212 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1213 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1217 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1218 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1224 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1231 muall[0][tel3+XX] = M_av[XX];
1232 muall[0][tel3+YY] = M_av[YY];
1233 muall[0][tel3+ZZ] = M_av[ZZ];
1236 /* Write to file the total dipole moment of the box, and its components
1239 if ((skip == 0) || ((teller % skip) == 0))
1241 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1242 t, M_av[XX], M_av[YY], M_av[ZZ],
1243 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1246 for (m = 0; (m < DIM); m++)
1250 M4[m] += sqr(M_av2[m]);
1252 /* Increment loop counter */
1255 /* Calculate for output the running averages */
1256 invtel = 1.0/teller;
1257 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1258 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1259 M_diff = M2_ave - M_ave2;
1261 /* Compute volume from box in traj, else we use the one from above */
1268 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1270 /* Calculate running average for dipole */
1273 mu_aver = (mu_ave/gnx_tot)*invtel;
1276 if ((skip == 0) || ((teller % skip) == 0))
1278 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1279 * the two. Here M is sum mu_i. Further write the finite system
1280 * Kirkwood G factor and epsilon.
1282 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1283 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1288 gmx_stats_get_average(muframelsq, &aver);
1289 fprintf(adip, "%10g %f \n", t, aver);
1292 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1294 if (!bMU || (mu_aver != -1))
1296 /* Finite system Kirkwood G-factor */
1297 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1298 /* Infinite system Kirkwood G-factor */
1299 if (epsilonRF == 0.0)
1301 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1305 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1306 Gk/(3*epsilon*(2*epsilonRF+1)));
1309 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1314 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1317 gmx_stats_done(muframelsq);
1321 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1325 bCont = read_next_x(oenv, status, &t, x, box);
1327 timecheck = check_times(t);
1329 while (bCont && (timecheck == 0) );
1331 gmx_rmpbc_done(gpbc);
1354 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1355 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1356 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1357 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1358 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1364 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1365 sfree(slab_dipoles);
1369 printf("Average volume over run is %g\n", vol_aver);
1372 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1373 print_cmap(cmap, gkrbin, nlevels);
1375 /* Autocorrelation function */
1380 printf("Not enough frames for autocorrelation\n");
1384 dt = (t1 - t0)/(teller-1);
1385 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1391 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1392 teller, 1, muall, dt, mode, TRUE);
1396 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1397 teller, gnx_tot, muall, dt,
1398 mode, strcmp(corrtype, "molsep"));
1404 real aver, sigma, error;
1406 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1407 printf("\nDipole moment (Debye)\n");
1408 printf("---------------------\n");
1409 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1410 aver, sigma, error);
1414 for (m = 0; (m < DIM); m++)
1416 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1419 printf("\nQuadrupole moment (Debye-Ang)\n");
1420 printf("-----------------------------\n");
1421 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1422 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1423 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1427 printf("The following averages for the complete trajectory have been calculated:\n\n");
1428 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1429 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1430 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1432 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1433 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1434 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1436 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1437 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1439 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1441 if (!bMU || (mu_aver != -1))
1443 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1444 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1446 printf("Epsilon = %g\n", epsilon);
1450 /* Write to file the dipole moment distibution during the simulation.
1452 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1453 for (i = 0; (i < ndipbin); i++)
1455 fprintf(outdd, "%10g %10f\n",
1456 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1463 done_gkrbin(&gkrbin);
1467 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1476 while (m < mols->nr && index[i] != mols->index[m])
1482 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1484 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1486 if (i >= *n || index[i] != j)
1488 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1492 /* Modify the index in place */
1495 printf("There are %d molecules in the selection\n", nmol);
1499 int gmx_dipoles(int argc, char *argv[])
1501 const char *desc[] = {
1502 "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1503 "system. From this you can compute e.g. the dielectric constant for",
1504 "low-dielectric media.",
1505 "For molecules with a net charge, the net charge is subtracted at",
1506 "center of mass of the molecule.[PAR]",
1507 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1508 "components as well as the norm of the vector.",
1509 "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",
1511 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1513 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1514 "Furthermore, the dipole autocorrelation function will be computed when",
1515 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1517 "The correlation functions can be averaged over all molecules",
1518 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1519 "or it can be computed over the total dipole moment of the simulation box",
1520 "([TT]total[tt]).[PAR]",
1521 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1522 "G-factor, as well as the average cosine of the angle between the dipoles",
1523 "as a function of the distance. The plot also includes gOO and hOO",
1524 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1525 "we also include the energy per scale computed by taking the inner product of",
1526 "the dipoles divided by the distance to the third power.[PAR]",
1529 "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1530 "This will calculate the autocorrelation function of the molecular",
1531 "dipoles using a first order Legendre polynomial of the angle of the",
1532 "dipole vector and itself a time t later. For this calculation 1001",
1533 "frames will be used. Further, the dielectric constant will be calculated",
1534 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1535 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1536 "distribution function a maximum of 5.0 will be used."
1538 real mu_max = 5, mu_aver = -1, rcmax = 0;
1539 real epsilonRF = 0.0, temp = 300;
1540 gmx_bool bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1541 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1542 const char *axtitle = "Z";
1543 int nslices = 10; /* nr of slices defined */
1544 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1545 int nlevels = 20, ndegrees = 90;
1548 { "-mu", FALSE, etREAL, {&mu_aver},
1549 "dipole of a single molecule (in Debye)" },
1550 { "-mumax", FALSE, etREAL, {&mu_max},
1551 "max dipole in Debye (for histogram)" },
1552 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1553 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1554 { "-skip", FALSE, etINT, {&skip},
1555 "Skip steps in the output (but not in the computations)" },
1556 { "-temp", FALSE, etREAL, {&temp},
1557 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1558 { "-corr", FALSE, etENUM, {corrtype},
1559 "Correlation function to calculate" },
1560 { "-pairs", FALSE, etBOOL, {&bPairs},
1561 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1562 { "-quad", FALSE, etBOOL, {&bQuad},
1563 "Take quadrupole into account"},
1564 { "-ncos", FALSE, etINT, {&ncos},
1565 "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." },
1566 { "-axis", FALSE, etSTR, {&axtitle},
1567 "Take the normal on the computational box in direction X, Y or Z." },
1568 { "-sl", FALSE, etINT, {&nslices},
1569 "Divide the box into this number of slices." },
1570 { "-gkratom", FALSE, etINT, {&nFA},
1571 "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" },
1572 { "-gkratom2", FALSE, etINT, {&nFB},
1573 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1574 { "-rcmax", FALSE, etREAL, {&rcmax},
1575 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1576 { "-phi", FALSE, etBOOL, {&bPhi},
1577 "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." },
1578 { "-nlevels", FALSE, etINT, {&nlevels},
1579 "Number of colors in the cmap output" },
1580 { "-ndegrees", FALSE, etINT, {&ndegrees},
1581 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1586 char **grpname = NULL;
1587 gmx_bool bGkr, bMU, bSlab;
1589 { efEDR, "-en", NULL, ffOPTRD },
1590 { efTRX, "-f", NULL, ffREAD },
1591 { efTPX, NULL, NULL, ffREAD },
1592 { efNDX, NULL, NULL, ffOPTRD },
1593 { efXVG, "-o", "Mtot", ffWRITE },
1594 { efXVG, "-eps", "epsilon", ffWRITE },
1595 { efXVG, "-a", "aver", ffWRITE },
1596 { efXVG, "-d", "dipdist", ffWRITE },
1597 { efXVG, "-c", "dipcorr", ffOPTWR },
1598 { efXVG, "-g", "gkr", ffOPTWR },
1599 { efXVG, "-adip", "adip", ffOPTWR },
1600 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1601 { efXVG, "-cos", "cosaver", ffOPTWR },
1602 { efXPM, "-cmap", "cmap", ffOPTWR },
1603 { efXVG, "-slab", "slab", ffOPTWR }
1605 #define NFILE asize(fnm)
1614 ppa = add_acf_pargs(&npargs, pa);
1615 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1616 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1621 printf("Using %g as mu_max and %g as the dipole moment.\n",
1623 if (epsilonRF == 0.0)
1625 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1628 bMU = opt2bSet("-en", NFILE, fnm);
1631 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.");
1633 bGkr = opt2bSet("-g", NFILE, fnm);
1634 if (opt2parg_bSet("-ncos", asize(pa), pa))
1636 if ((ncos != 1) && (ncos != 2))
1638 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1642 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1643 opt2parg_bSet("-axis", asize(pa), pa));
1648 printf("WARNING: Can not determine quadrupoles from energy file\n");
1653 printf("WARNING: Can not determine Gk(r) from energy file\n");
1659 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1660 " not enter a valid dipole for the molecules\n");
1665 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1666 &natoms, NULL, NULL, NULL, top);
1669 snew(grpname, ncos);
1670 snew(grpindex, ncos);
1671 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1672 ncos, gnx, grpindex, grpname);
1673 for (k = 0; (k < ncos); k++)
1675 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1676 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1680 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1681 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1682 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1683 opt2fn_null("-cos", NFILE, fnm),
1684 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1685 bPairs, corrtype[0],
1686 opt2fn("-c", NFILE, fnm),
1687 bGkr, opt2fn("-g", NFILE, fnm),
1688 bPhi, &nlevels, ndegrees,
1690 opt2fn("-cmap", NFILE, fnm), rcmax,
1691 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1692 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1693 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1695 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1696 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1697 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1698 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1699 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");