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
51 #include "gmx_statistics.h"
64 #define e2d(x) ENM2DEBYE*(x)
65 #define EANG2CM E_CHARGE*1.0e-10 /* e Angstrom to Coulomb meter */
66 #define CM2D SPEED_OF_LIGHT*1.0e+24 /* Coulomb meter to Debye */
78 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
86 if ((ptr = getenv("GKRWIDTH")) != NULL)
90 sscanf(ptr, "%lf", &bw);
95 gb->spacing = 0.01; /* nm */
97 gb->nelem = 1 + radius/gb->spacing;
104 gb->nx = 1 + rcmax/gb->spacing;
107 snew(gb->elem, gb->nelem);
108 snew(gb->count, gb->nelem);
110 snew(gb->cmap, gb->nx);
111 gb->ny = 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 = (M_PI+phi)*gb->ny/(2*M_PI);
148 cy = (alpha*gb->ny)/M_PI; /*((1+cosa)*0.5*(gb->ny));*/
150 cy = min(gb->ny-1, 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, index, grp0, grp1;
175 real qtot, q, r2, 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 || (cosa != 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)
271 for (i = 0; (i < gb->nx); i++)
273 vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
274 for (j = 0; (j < gb->ny); j++)
276 gb->cmap[i][j] /= vol;
277 hi = max(hi, gb->cmap[i][j]);
282 gmx_fatal(FARGS, "No data in the cmap");
287 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
294 t_rgb rlo = { 1, 1, 1 };
295 t_rgb rhi = { 0, 0, 0 };
297 hi = normalize_cmap(gb);
298 snew(xaxis, gb->nx+1);
299 for (i = 0; (i < gb->nx+1); i++)
301 xaxis[i] = i*gb->spacing;
304 for (j = 0; (j < gb->ny); j++)
308 yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
312 yaxis[j] = (180.0*j)/(gb->ny-1.0);
314 /*2.0*j/(gb->ny-1.0)-1.0;*/
316 out = ffopen(cmap, "w");
318 "Dipole Orientation Distribution", "Fraction", "r (nm)",
319 gb->bPhi ? "Phi" : "Alpha",
320 gb->nx, gb->ny, xaxis, yaxis,
321 gb->cmap, 0, hi, rlo, rhi, nlevels);
327 static void print_gkrbin(const char *fn, t_gkrbin *gb,
328 int ngrp, int nframes, real volume,
329 const output_env_t oenv)
331 /* We compute Gk(r), gOO and hOO according to
332 * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
333 * In this implementation the angle between dipoles is stored
334 * rather than their inner product. This allows to take polarizible
335 * models into account. The RDF is calculated as well, almost for free!
338 const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
340 real x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
343 fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
344 xvgr_legend(fp, asize(leg), leg, oenv);
346 Gkr = 1; /* Self-dipole inproduct = 1 */
351 fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
352 fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
356 while (last > 1 && gb->elem[last-1] == 0)
361 /* Divide by dipole squared, by number of frames, by number of origins.
362 * Multiply by 2 because we only take half the matrix of interactions
365 fac = 2.0/((double) ngrp * (double) nframes);
368 for (i = 0; i < last; i++)
370 /* Centre of the coordinate in the spherical layer */
373 /* Volume of the layer */
374 vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
377 gOO = gb->count[i]*fac/(rho*vol_s);
379 /* Dipole correlation hOO, normalized by the relative number density, like
380 * in a Radial distribution function.
382 ggg = gb->elem[i]*fac;
383 hOO = 3.0*ggg/(rho*vol_s);
387 cosav = gb->elem[i]/gb->count[i];
393 ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
395 fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
396 x1, Gkr, cosav, hOO, gOO, ener);
404 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
405 real *t, int nre, t_enxframe *fr)
411 bCont = do_enx(fmu, fr);
414 fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
415 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
420 if (Vol != -1) /* we've got Volume in the energy file */
422 *vol = fr->ener[Vol].e;
424 for (i = 0; i < DIM; i++)
426 mu[i] = fr->ener[iMu[i]].e;
434 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
437 int ncharged, m, a0, a1, a;
440 for (m = 0; m < n; m++)
442 a0 = mols->index[index[m]];
443 a1 = mols->index[index[m]+1];
446 for (a = a0; a < a1; a++)
451 /* This check is only for the count print */
452 if (fabs(qtot) > 0.01)
458 /* Remove the net charge at the center of mass */
459 for (a = a0; a < a1; a++)
461 atom[a].q -= qtot*atom[a].m/mtot;
468 printf("There are %d charged molecules in the selection,\n"
469 "will subtract their charge at their center of mass\n", ncharged);
473 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
479 for (k = k0; (k < k1); k++)
482 for (m = 0; (m < DIM); m++)
489 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
491 int i, k, m, n, niter;
492 real q, r2, mass, masstot;
493 rvec com; /* center of mass */
494 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, dc1, d, n5, 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, j, k, n, m, natom = 0, nmol, gnx_tot, teller, tel3;
784 int *dipole_bin, ndipbin, ibin, iVol, step, idim = -1;
787 real rcut = 0, t, t0, t1, dt, lambda, dd, rms_cos;
790 gmx_bool bCorr, bTotal, bCont;
791 double M_diff = 0, epsilon, invtel, vol_aver;
792 double mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
793 double M[3], M2[3], M4[3], Gk = 0, g_k = 0;
794 gmx_stats_t Mx, My, Mz, Msq, Vol, *Qlsq, mulsq, muframelsq = NULL;
797 rvec *slab_dipoles = NULL;
799 t_block *mols = NULL;
800 gmx_rmpbc_t gpbc = NULL;
812 /* initialize to a negative value so we can see that it wasn't picked up */
813 iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
816 fmu = open_enx(mufn, "r");
817 do_enxnms(fmu, &nre, &enm);
819 /* Determine the indexes of the energy grps we need */
820 for (i = 0; (i < nre); i++)
822 if (strstr(enm[i].name, "Volume"))
826 else if (strstr(enm[i].name, "Mu-X"))
830 else if (strstr(enm[i].name, "Mu-Y"))
834 else if (strstr(enm[i].name, "Mu-Z"))
839 if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
841 gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
846 atom = top->atoms.atom;
849 if ((iVol == -1) && bMU)
851 printf("Using Volume from topology: %g nm^3\n", volume);
854 /* Correlation stuff */
855 bCorr = (corrtype[0] != 'n');
856 bTotal = (corrtype[0] == 't');
862 snew(muall[0], nframes*DIM);
867 for (i = 0; (i < gnx[0]); i++)
869 snew(muall[i], nframes*DIM);
874 /* Allocate array which contains for every molecule in a frame the
879 snew(dipole, gnx_tot);
884 for (i = 0; (i < DIM); i++)
886 Qlsq[i] = gmx_stats_init();
888 mulsq = gmx_stats_init();
890 /* Open all the files */
891 outmtot = xvgropen(out_mtot,
892 "Total dipole moment of the simulation box vs. time",
893 "Time (ps)", "Total Dipole Moment (Debye)", oenv);
894 outeps = xvgropen(out_eps, "Epsilon and Kirkwood factors",
895 "Time (ps)", "", oenv);
896 outaver = xvgropen(out_aver, "Total dipole moment",
897 "Time (ps)", "D", oenv);
900 idim = axtitle[0] - 'X';
901 if ((idim < 0) || (idim >= DIM))
903 idim = axtitle[0] - 'x';
905 if ((idim < 0) || (idim >= DIM))
913 fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
914 axtitle, nslices, idim);
917 snew(slab_dipoles, nslices);
918 fprintf(stderr, "Doing slab analysis\n");
924 adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
925 xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
930 caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
931 "Average absolute dipole orientation", "Time (ps)", "", oenv);
932 xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
938 snew(dipsp, gnx_tot);
940 /* we need a dummy file for gnuplot */
941 dip3d = (FILE *)ffopen("dummy.dat", "w");
942 fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
945 dip3d = (FILE *)ffopen(fndip3d, "w");
946 fprintf(dip3d, "# This file was created by %s\n", Program());
947 fprintf(dip3d, "# which is part of G R O M A C S:\n");
948 fprintf(dip3d, "#\n");
951 /* Write legends to all the files */
952 xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
953 xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
955 if (bMU && (mu_aver == -1))
957 xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
961 xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
967 /* Read the first frame from energy or traj file */
972 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
975 timecheck = check_times(t);
980 if ((teller % 10) == 0)
982 fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
987 printf("End of %s reached\n", mufn);
991 while (bCont && (timecheck < 0));
995 natom = read_first_x(oenv, &status, fn, &t, &x, box);
998 /* Calculate spacing for dipole bin (simple histogram) */
999 ndipbin = 1+(mu_max/0.01);
1000 snew(dipole_bin, ndipbin);
1003 for (m = 0; (m < DIM); m++)
1005 M[m] = M2[m] = M4[m] = 0.0;
1010 /* Use 0.7 iso 0.5 to account for pressure scaling */
1011 /* rcut = 0.7*sqrt(max_cutoff2(box)); */
1012 rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1014 gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1016 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom, box);
1018 /* Start while loop over frames */
1023 if (bCorr && (teller >= nframes))
1028 srenew(muall[0], nframes*DIM);
1032 for (i = 0; (i < gnx_tot); i++)
1034 srenew(muall[i], nframes*DIM);
1040 muframelsq = gmx_stats_init();
1043 for (m = 0; (m < DIM); m++)
1050 /* Copy rvec into double precision local variable */
1051 for (m = 0; (m < DIM); m++)
1059 for (m = 0; (m < DIM); m++)
1064 gmx_rmpbc(gpbc, natom, box, x);
1066 /* Begin loop of all molecules in frame */
1067 for (n = 0; (n < ncos); n++)
1069 for (i = 0; (i < gnx[n]); i++)
1073 ind0 = mols->index[molindex[n][i]];
1074 ind1 = mols->index[molindex[n][i]+1];
1076 mol_dip(ind0, ind1, x, atom, dipole[i]);
1077 gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1078 gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1081 update_slab_dipoles(ind0, ind1, x,
1082 dipole[i], idim, nslices, slab_dipoles, box);
1086 mol_quad(ind0, ind1, x, atom, quad);
1087 for (m = 0; (m < DIM); m++)
1089 gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1092 if (bCorr && !bTotal)
1095 muall[i][tel3+XX] = dipole[i][XX];
1096 muall[i][tel3+YY] = dipole[i][YY];
1097 muall[i][tel3+ZZ] = dipole[i][ZZ];
1100 for (m = 0; (m < DIM); m++)
1102 M_av[m] += dipole[i][m]; /* M per frame */
1103 mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
1105 mu_mol = sqrt(mu_mol);
1107 mu_ave += mu_mol; /* calc. the average mu */
1109 /* Update the dipole distribution */
1110 ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1118 rvec2sprvec(dipole[i], dipsp[i]);
1120 if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1122 if (dipsp[i][ZZ] < 0.5 * M_PI)
1131 else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1133 if (dipsp[i][ZZ] < 0.5 * M_PI)
1142 else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1144 if (dipsp[i][ZZ] < 0.5 * M_PI)
1153 else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1155 if (dipsp[i][ZZ] < 0.5 * M_PI)
1166 fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d # %d %d\n",
1171 x[ind0][XX]+dipole[i][XX]/25,
1172 x[ind0][YY]+dipole[i][YY]/25,
1173 x[ind0][ZZ]+dipole[i][ZZ]/25,
1177 } /* End loop of all molecules in frame */
1181 fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1182 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1183 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1184 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1185 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1186 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1190 /* Compute square of total dipole */
1191 for (m = 0; (m < DIM); m++)
1193 M_av2[m] = M_av[m]*M_av[m];
1198 compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1199 rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1200 sqr(dipaxis[YY]-0.5)+
1201 sqr(dipaxis[ZZ]-0.5));
1204 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
1205 t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1209 fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e\n",
1210 t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1216 do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1223 muall[0][tel3+XX] = M_av[XX];
1224 muall[0][tel3+YY] = M_av[YY];
1225 muall[0][tel3+ZZ] = M_av[ZZ];
1228 /* Write to file the total dipole moment of the box, and its components
1231 if ((skip == 0) || ((teller % skip) == 0))
1233 fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
1234 t, M_av[XX], M_av[YY], M_av[ZZ],
1235 sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1238 for (m = 0; (m < DIM); m++)
1242 M4[m] += sqr(M_av2[m]);
1244 /* Increment loop counter */
1247 /* Calculate for output the running averages */
1248 invtel = 1.0/teller;
1249 M2_ave = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1250 M_ave2 = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1251 M_diff = M2_ave - M_ave2;
1253 /* Compute volume from box in traj, else we use the one from above */
1260 epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1262 /* Calculate running average for dipole */
1265 mu_aver = (mu_ave/gnx_tot)*invtel;
1268 if ((skip == 0) || ((teller % skip) == 0))
1270 /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1271 * the two. Here M is sum mu_i. Further write the finite system
1272 * Kirkwood G factor and epsilon.
1274 fprintf(outaver, "%10g %10.3e %10.3e %10.3e %10.3e\n",
1275 t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1280 gmx_stats_get_average(muframelsq, &aver);
1281 fprintf(adip, "%10g %f \n", t, aver);
1284 printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1286 if (!bMU || (mu_aver != -1))
1288 /* Finite system Kirkwood G-factor */
1289 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1290 /* Infinite system Kirkwood G-factor */
1291 if (epsilonRF == 0.0)
1293 g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1297 g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1298 Gk/(3*epsilon*(2*epsilonRF+1)));
1301 fprintf(outeps, "%10g %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1306 fprintf(outeps, "%10g %12.8e\n", t, epsilon);
1309 gmx_stats_done(muframelsq);
1313 bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1317 bCont = read_next_x(oenv, status, &t, natom, x, box);
1319 timecheck = check_times(t);
1321 while (bCont && (timecheck == 0) );
1323 gmx_rmpbc_done(gpbc);
1346 fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1347 fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1348 fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1349 fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1350 fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1356 dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1357 sfree(slab_dipoles);
1361 printf("Average volume over run is %g\n", vol_aver);
1364 print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1365 print_cmap(cmap, gkrbin, nlevels);
1367 /* Autocorrelation function */
1372 printf("Not enough frames for autocorrelation\n");
1376 dt = (t1 - t0)/(teller-1);
1377 printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1383 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1384 teller, 1, muall, dt, mode, TRUE);
1388 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1389 teller, gnx_tot, muall, dt,
1390 mode, strcmp(corrtype, "molsep"));
1396 real aver, sigma, error, lsq;
1398 gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1399 printf("\nDipole moment (Debye)\n");
1400 printf("---------------------\n");
1401 printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n",
1402 aver, sigma, error);
1407 for (m = 0; (m < DIM); m++)
1409 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1412 printf("\nQuadrupole moment (Debye-Ang)\n");
1413 printf("-----------------------------\n");
1414 printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
1415 printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
1416 printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
1420 printf("The following averages for the complete trajectory have been calculated:\n\n");
1421 printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1422 printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1423 printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1425 printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1426 printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1427 printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1429 printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1430 printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1432 printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1434 if (!bMU || (mu_aver != -1))
1436 printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1437 printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1439 printf("Epsilon = %g\n", epsilon);
1443 /* Write to file the dipole moment distibution during the simulation.
1445 outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1446 for (i = 0; (i < ndipbin); i++)
1448 fprintf(outdd, "%10g %10f\n",
1449 (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1456 done_gkrbin(&gkrbin);
1460 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1469 while (m < mols->nr && index[i] != mols->index[m])
1475 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1477 for (j = mols->index[m]; j < mols->index[m+1]; j++)
1479 if (i >= *n || index[i] != j)
1481 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1485 /* Modify the index in place */
1488 printf("There are %d molecules in the selection\n", nmol);
1492 int gmx_dipoles(int argc, char *argv[])
1494 const char *desc[] = {
1495 "[TT]g_dipoles[tt] computes the total dipole plus fluctuations of a simulation",
1496 "system. From this you can compute e.g. the dielectric constant for",
1497 "low-dielectric media.",
1498 "For molecules with a net charge, the net charge is subtracted at",
1499 "center of mass of the molecule.[PAR]",
1500 "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1501 "components as well as the norm of the vector.",
1502 "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",
1504 "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1506 "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1507 "Furthermore, the dipole autocorrelation function will be computed when",
1508 "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1510 "The correlation functions can be averaged over all molecules",
1511 "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1512 "or it can be computed over the total dipole moment of the simulation box",
1513 "([TT]total[tt]).[PAR]",
1514 "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1515 "G-factor, as well as the average cosine of the angle between the dipoles",
1516 "as a function of the distance. The plot also includes gOO and hOO",
1517 "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1518 "we also include the energy per scale computed by taking the inner product of",
1519 "the dipoles divided by the distance to the third power.[PAR]",
1522 "[TT]g_dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1523 "This will calculate the autocorrelation function of the molecular",
1524 "dipoles using a first order Legendre polynomial of the angle of the",
1525 "dipole vector and itself a time t later. For this calculation 1001",
1526 "frames will be used. Further, the dielectric constant will be calculated",
1527 "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1528 "an average dipole moment of the molecule of 2.273 (SPC). For the",
1529 "distribution function a maximum of 5.0 will be used."
1531 real mu_max = 5, mu_aver = -1, rcmax = 0;
1532 real epsilonRF = 0.0, temp = 300;
1533 gmx_bool bAverCorr = FALSE, bMolCorr = FALSE, bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1534 const char *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1535 const char *axtitle = "Z";
1536 int nslices = 10; /* nr of slices defined */
1537 int skip = 0, nFA = 0, nFB = 0, ncos = 1;
1538 int nlevels = 20, ndegrees = 90;
1541 { "-mu", FALSE, etREAL, {&mu_aver},
1542 "dipole of a single molecule (in Debye)" },
1543 { "-mumax", FALSE, etREAL, {&mu_max},
1544 "max dipole in Debye (for histogram)" },
1545 { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1546 "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1547 { "-skip", FALSE, etINT, {&skip},
1548 "Skip steps in the output (but not in the computations)" },
1549 { "-temp", FALSE, etREAL, {&temp},
1550 "Average temperature of the simulation (needed for dielectric constant calculation)" },
1551 { "-corr", FALSE, etENUM, {corrtype},
1552 "Correlation function to calculate" },
1553 { "-pairs", FALSE, etBOOL, {&bPairs},
1554 "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1555 { "-quad", FALSE, etBOOL, {&bQuad},
1556 "Take quadrupole into account"},
1557 { "-ncos", FALSE, etINT, {&ncos},
1558 "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." },
1559 { "-axis", FALSE, etSTR, {&axtitle},
1560 "Take the normal on the computational box in direction X, Y or Z." },
1561 { "-sl", FALSE, etINT, {&nslices},
1562 "Divide the box into this number of slices." },
1563 { "-gkratom", FALSE, etINT, {&nFA},
1564 "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" },
1565 { "-gkratom2", FALSE, etINT, {&nFB},
1566 "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1567 { "-rcmax", FALSE, etREAL, {&rcmax},
1568 "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1569 { "-phi", FALSE, etBOOL, {&bPhi},
1570 "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." },
1571 { "-nlevels", FALSE, etINT, {&nlevels},
1572 "Number of colors in the cmap output" },
1573 { "-ndegrees", FALSE, etINT, {&ndegrees},
1574 "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1579 char **grpname = NULL;
1580 gmx_bool bCorr, bGkr, bMU, bSlab;
1582 { efEDR, "-en", NULL, ffOPTRD },
1583 { efTRX, "-f", NULL, ffREAD },
1584 { efTPX, NULL, NULL, ffREAD },
1585 { efNDX, NULL, NULL, ffOPTRD },
1586 { efXVG, "-o", "Mtot", ffWRITE },
1587 { efXVG, "-eps", "epsilon", ffWRITE },
1588 { efXVG, "-a", "aver", ffWRITE },
1589 { efXVG, "-d", "dipdist", ffWRITE },
1590 { efXVG, "-c", "dipcorr", ffOPTWR },
1591 { efXVG, "-g", "gkr", ffOPTWR },
1592 { efXVG, "-adip", "adip", ffOPTWR },
1593 { efXVG, "-dip3d", "dip3d", ffOPTWR },
1594 { efXVG, "-cos", "cosaver", ffOPTWR },
1595 { efXPM, "-cmap", "cmap", ffOPTWR },
1596 { efXVG, "-slab", "slab", ffOPTWR }
1598 #define NFILE asize(fnm)
1607 ppa = add_acf_pargs(&npargs, pa);
1608 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1609 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv);
1611 printf("Using %g as mu_max and %g as the dipole moment.\n",
1613 if (epsilonRF == 0.0)
1615 printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1618 bMU = opt2bSet("-en", NFILE, fnm);
1621 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.");
1623 bGkr = opt2bSet("-g", NFILE, fnm);
1624 if (opt2parg_bSet("-ncos", asize(pa), pa))
1626 if ((ncos != 1) && (ncos != 2))
1628 gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1632 bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1633 opt2parg_bSet("-axis", asize(pa), pa));
1639 printf("WARNING: Can not determine quadrupoles from energy file\n");
1644 printf("WARNING: Can not determine Gk(r) from energy file\n");
1650 printf("WARNING: Can not calculate Gk and gk, since you did\n"
1651 " not enter a valid dipole for the molecules\n");
1656 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1657 &natoms, NULL, NULL, NULL, top);
1660 snew(grpname, ncos);
1661 snew(grpindex, ncos);
1662 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1663 ncos, gnx, grpindex, grpname);
1664 for (k = 0; (k < ncos); k++)
1666 dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1667 neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1671 do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1672 opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1673 opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1674 opt2fn_null("-cos", NFILE, fnm),
1675 opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1676 bPairs, corrtype[0],
1677 opt2fn("-c", NFILE, fnm),
1678 bGkr, opt2fn("-g", NFILE, fnm),
1679 bPhi, &nlevels, ndegrees,
1681 opt2fn("-cmap", NFILE, fnm), rcmax,
1682 bQuad, bMU, opt2fn("-en", NFILE, fnm),
1683 gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1684 bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1686 do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1687 do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1688 do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1689 do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1690 do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");