Remove unused function generate_excls and make clean_excls static

[gromacs.git] / src / gromacs / gmxana / gmx_dipoles.cpp

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team.
 * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 * and including many others, as listed in the AUTHORS file in the
 * top-level source directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */
#include "gmxpre.h"

#include <cmath>
#include <cstring>

#include <algorithm>

#include "gromacs/commandline/pargs.h"
#include "gromacs/commandline/viewit.h"
#include "gromacs/correlationfunctions/autocorr.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/linearalgebra/nrjac.h"
#include "gromacs/listed_forces/bonded.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/rmpbc.h"
#include "gromacs/statistics/statistics.h"
#include "gromacs/topology/index.h"
#include "gromacs/topology/topology.h"
#include "gromacs/trajectory/energyframe.h"
#include "gromacs/utility/arraysize.h"
#include "gromacs/utility/binaryinformation.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"

#define e2d(x) ENM2DEBYE*(x)
#define EANG2CM  (E_CHARGE*1.0e-10)       /* e Angstrom to Coulomb meter */
#define CM2D  (SPEED_OF_LIGHT*1.0e+24)    /* Coulomb meter to Debye */

typedef struct {
    int      nelem;
    real     spacing, radius;
    real    *elem;
    int     *count;
    gmx_bool bPhi;
    int      nx, ny;
    real   **cmap;
} t_gkrbin;

static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
{
    t_gkrbin *gb;
    char     *ptr;
    int       i;

    snew(gb, 1);

    if ((ptr = getenv("GMX_DIPOLE_SPACING")) != nullptr)
    {
        double bw = strtod(ptr, nullptr);
        gb->spacing = bw;
    }
    else
    {
        gb->spacing = 0.01; /* nm */
    }
    gb->nelem   = 1 + static_cast<int>(radius/gb->spacing);
    if (rcmax == 0)
    {
        gb->nx = gb->nelem;
    }
    else
    {
        gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
    }
    gb->radius  = radius;
    snew(gb->elem, gb->nelem);
    snew(gb->count, gb->nelem);

    snew(gb->cmap, gb->nx);
    gb->ny = std::max(2, ndegrees);
    for (i = 0; (i < gb->nx); i++)
    {
        snew(gb->cmap[i], gb->ny);
    }
    gb->bPhi = bPhi;

    return gb;
}

static void done_gkrbin(t_gkrbin **gb)
{
    sfree((*gb)->elem);
    sfree((*gb)->count);
    sfree((*gb));
    *gb = nullptr;
}

static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
{
    int  cy, index = gmx::roundToInt(r/gb->spacing);
    real alpha;

    if (index < gb->nelem)
    {
        gb->elem[index]  += cosa;
        gb->count[index]++;
    }
    if (index < gb->nx)
    {
        alpha = std::acos(cosa);
        if (gb->bPhi)
        {
            cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
        }
        else
        {
            cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
        }
        cy = std::min(gb->ny-1, std::max(0, cy));
        if (debug)
        {
            fprintf(debug, "CY: %10f  %5d\n", alpha, cy);
        }
        gb->cmap[index][cy] += 1;
    }
}

static void rvec2sprvec(rvec dipcart, rvec dipsp)
{
    clear_rvec(dipsp);
    dipsp[0] = std::sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]);  /* R */
    dipsp[1] = std::atan2(dipcart[YY], dipcart[XX]);                                                /* Theta */
    dipsp[2] = std::atan2(std::sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
}



static void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
                   const int mindex[], rvec x[], rvec mu[],
                   int ePBC, const matrix box, const t_atom *atom, const int *nAtom)
{
    static rvec *xcm[2] = { nullptr, nullptr};
    int          gi, gj, j0, j1, i, j, k, n, grp0, grp1;
    real         qtot, q, cosa, rr, phi;
    rvec         dx;
    t_pbc        pbc;

    for (n = 0; (n < ncos); n++)
    {
        if (!xcm[n])
        {
            snew(xcm[n], ngrp[n]);
        }
        for (i = 0; (i < ngrp[n]); i++)
        {
            /* Calculate center of mass of molecule */
            gi = molindex[n][i];
            j0 = mindex[gi];

            if (nAtom[n] > 0)
            {
                copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
            }
            else
            {
                j1 = mindex[gi+1];
                clear_rvec(xcm[n][i]);
                qtot = 0;
                for (j = j0; j < j1; j++)
                {
                    q     = std::abs(atom[j].q);
                    qtot += q;
                    for (k = 0; k < DIM; k++)
                    {
                        xcm[n][i][k] += q*x[j][k];
                    }
                }
                svmul(1/qtot, xcm[n][i], xcm[n][i]);
            }
        }
    }
    set_pbc(&pbc, ePBC, box);
    grp0 = 0;
    grp1 = ncos-1;
    for (i = 0; i < ngrp[grp0]; i++)
    {
        gi = molindex[grp0][i];
        j0 = (ncos == 2) ? 0 : i+1;
        for (j = j0; j < ngrp[grp1]; j++)
        {
            gj   = molindex[grp1][j];
            if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
            {
                /* Compute distance between molecules including PBC */
                pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
                rr = norm(dx);

                if (gb->bPhi)
                {
                    rvec xi, xj, xk, xl;
                    rvec r_ij, r_kj, r_kl, mm, nn;
                    int  t1, t2, t3;

                    copy_rvec(xcm[grp0][i], xj);
                    copy_rvec(xcm[grp1][j], xk);
                    rvec_add(xj, mu[gi], xi);
                    rvec_add(xk, mu[gj], xl);
                    phi = dih_angle(xi, xj, xk, xl, &pbc,
                                    r_ij, r_kj, r_kl, mm, nn, /* out */
                                    &t1, &t2, &t3);
                    cosa = std::cos(phi);
                }
                else
                {
                    cosa = cos_angle(mu[gi], mu[gj]);
                    phi  = 0;
                }
                if (debug || std::isnan(cosa))
                {
                    fprintf(debug ? debug : stderr,
                            "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",
                            gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
                            gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
                            rr, cosa);
                }

                add2gkr(gb, rr, cosa, phi);
            }
        }
    }
}

static real normalize_cmap(t_gkrbin *gb)
{
    int    i, j;
    real   hi;
    double vol;

    hi = 0;
    for (i = 0; (i < gb->nx); i++)
    {
        vol = 4*M_PI*gmx::square(gb->spacing*i)*gb->spacing;
        for (j = 0; (j < gb->ny); j++)
        {
            gb->cmap[i][j] /= vol;
            hi              = std::max(hi, gb->cmap[i][j]);
        }
    }
    if (hi <= 0)
    {
        gmx_fatal(FARGS, "No data in the cmap");
    }
    return hi;
}

static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
{
    FILE   *out;
    int     i, j;
    real    hi;

    real   *xaxis, *yaxis;
    t_rgb   rlo = { 1, 1, 1 };
    t_rgb   rhi = { 0, 0, 0 };

    hi = normalize_cmap(gb);
    snew(xaxis, gb->nx+1);
    for (i = 0; (i < gb->nx+1); i++)
    {
        xaxis[i] = i*gb->spacing;
    }
    snew(yaxis, gb->ny);
    for (j = 0; (j < gb->ny); j++)
    {
        if (gb->bPhi)
        {
            yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
        }
        else
        {
            yaxis[j] = (180.0*j)/(gb->ny-1.0);
        }
        /*2.0*j/(gb->ny-1.0)-1.0;*/
    }
    out = gmx_ffopen(cmap, "w");
    write_xpm(out, 0,
              "Dipole Orientation Distribution", "Fraction", "r (nm)",
              gb->bPhi ? "Phi" : "Alpha",
              gb->nx, gb->ny, xaxis, yaxis,
              gb->cmap, 0, hi, rlo, rhi, nlevels);
    gmx_ffclose(out);
    sfree(xaxis);
    sfree(yaxis);
}

static void print_gkrbin(const char *fn, t_gkrbin *gb,
                         int ngrp, int nframes, real volume,
                         const gmx_output_env_t *oenv)
{
    /* We compute Gk(r), gOO and hOO according to
     * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
     * In this implementation the angle between dipoles is stored
     * rather than their inner product. This allows to take polarizible
     * models into account. The RDF is calculated as well, almost for free!
     */
    FILE       *fp;
    const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
    int         i, last;
    real        x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
    double      fac;

    fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
    xvgr_legend(fp, asize(leg), leg, oenv);

    Gkr = 1;  /* Self-dipole inproduct = 1 */
    rho = ngrp/volume;

    if (debug)
    {
        fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
        fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
    }

    last = gb->nelem-1;
    while (last > 1 && gb->elem[last-1] == 0)
    {
        last--;
    }

    /* Divide by dipole squared, by number of frames, by number of origins.
     * Multiply by 2 because we only take half the matrix of interactions
     * into account.
     */
    fac  = 2.0/(ngrp * nframes);

    x0 = 0;
    for (i = 0; i < last; i++)
    {
        /* Centre of the coordinate in the spherical layer */
        x1    = x0+gb->spacing;

        /* Volume of the layer */
        vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);

        /* gOO */
        gOO   = gb->count[i]*fac/(rho*vol_s);

        /* Dipole correlation hOO, normalized by the relative number density, like
         * in a Radial distribution function.
         */
        ggg  = gb->elem[i]*fac;
        hOO  = 3.0*ggg/(rho*vol_s);
        Gkr += ggg;
        if (gb->count[i])
        {
            cosav = gb->elem[i]/gb->count[i];
        }
        else
        {
            cosav = 0;
        }
        ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);

        fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e  %12.5e\n",
                x1, Gkr, cosav, hOO, gOO, ener);

        /* Swap x0 and x1 */
        x0 = x1;
    }
    xvgrclose(fp);
}

static gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, const ivec iMu, rvec mu, real *vol,
                                 real *t, int nre, t_enxframe *fr)
{
    int          i;
    gmx_bool     bCont;
    char         buf[22];

    bCont = do_enx(fmu, fr);
    if (fr->nre != nre)
    {
        fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
                nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
    }

    if (bCont)
    {
        if (Vol != -1)          /* we've got Volume in the energy file */
        {
            *vol = fr->ener[Vol].e;
        }
        for (i = 0; i < DIM; i++)
        {
            mu[i] = fr->ener[iMu[i]].e;
        }
        *t = fr->t;
    }

    return bCont;
}

static void neutralize_mols(int n, const int *index, const t_block *mols, t_atom *atom)
{
    double mtot, qtot;
    int    ncharged, m, a0, a1, a;

    ncharged = 0;
    for (m = 0; m < n; m++)
    {
        a0   = mols->index[index[m]];
        a1   = mols->index[index[m]+1];
        mtot = 0;
        qtot = 0;
        for (a = a0; a < a1; a++)
        {
            mtot += atom[a].m;
            qtot += atom[a].q;
        }
        /* This check is only for the count print */
        if (std::abs(qtot) > 0.01)
        {
            ncharged++;
        }
        if (mtot > 0)
        {
            /* Remove the net charge at the center of mass */
            for (a = a0; a < a1; a++)
            {
                atom[a].q -= qtot*atom[a].m/mtot;
            }
        }
    }

    if (ncharged)
    {
        printf("There are %d charged molecules in the selection,\n"
               "will subtract their charge at their center of mass\n", ncharged);
    }
}

static void mol_dip(int k0, int k1, rvec x[], const t_atom atom[], rvec mu)
{
    int  k, m;
    real q;

    clear_rvec(mu);
    for (k = k0; (k < k1); k++)
    {
        q  = e2d(atom[k].q);
        for (m = 0; (m < DIM); m++)
        {
            mu[m] += q*x[k][m];
        }
    }
}

static void mol_quad(int k0, int k1, rvec x[], const t_atom atom[], rvec quad)
{
    int      i, k, m, n, niter;
    real     q, r2, mass, masstot;
    rvec     com;        /* center of mass */
    rvec     r;          /* distance of atoms to center of mass */
    double **inten;
    double   dd[DIM], **ev;

    snew(inten, DIM);
    snew(ev, DIM);
    for (i = 0; (i < DIM); i++)
    {
        snew(inten[i], DIM);
        snew(ev[i], DIM);
        dd[i] = 0.0;
    }

    /* Compute center of mass */
    clear_rvec(com);
    masstot = 0.0;
    for (k = k0; (k < k1); k++)
    {
        mass     = atom[k].m;
        masstot += mass;
        for (i = 0; (i < DIM); i++)
        {
            com[i] += mass*x[k][i];
        }
    }
    svmul((1.0/masstot), com, com);

    /* We want traceless quadrupole moments, so let us calculate the complete
     * quadrupole moment tensor and diagonalize this tensor to get
     * the individual components on the diagonal.
     */

#define delta(a, b) (( (a) == (b) ) ? 1.0 : 0.0)

    for (m = 0; (m < DIM); m++)
    {
        for (n = 0; (n < DIM); n++)
        {
            inten[m][n] = 0;
        }
    }
    for (k = k0; (k < k1); k++)         /* loop over atoms in a molecule */
    {
        q  = (atom[k].q)*100.0;
        rvec_sub(x[k], com, r);
        r2 = iprod(r, r);
        for (m = 0; (m < DIM); m++)
        {
            for (n = 0; (n < DIM); n++)
            {
                inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
            }
        }
    }
    if (debug)
    {
        for (i = 0; (i < DIM); i++)
        {
            fprintf(debug, "Q[%d] = %8.3f  %8.3f  %8.3f\n",
                    i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
        }
    }

    /* We've got the quadrupole tensor, now diagonalize the sucker */
    jacobi(inten, 3, dd, ev, &niter);

    if (debug)
    {
        for (i = 0; (i < DIM); i++)
        {
            fprintf(debug, "ev[%d] = %8.3f  %8.3f  %8.3f\n",
                    i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
        }
        for (i = 0; (i < DIM); i++)
        {
            fprintf(debug, "Q'[%d] = %8.3f  %8.3f  %8.3f\n",
                    i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
        }
    }
    /* Sort the eigenvalues, for water we know that the order is as follows:
     *
     * Q_yy, Q_zz, Q_xx
     *
     * At the moment I have no idea how this will work out for other molecules...
     */

    if (dd[1] > dd[0])
    {
        std::swap(dd[0], dd[1]);
    }
    if (dd[2] > dd[1])
    {
        std::swap(dd[1], dd[2]);
    }
    if (dd[1] > dd[0])
    {
        std::swap(dd[0], dd[1]);
    }

    for (m = 0; (m < DIM); m++)
    {
        quad[0] = dd[2];  /* yy */
        quad[1] = dd[0];  /* zz */
        quad[2] = dd[1];  /* xx */
    }

    if (debug)
    {
        pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
    }

    /* clean-up */
    for (i = 0; (i < DIM); i++)
    {
        sfree(inten[i]);
        sfree(ev[i]);
    }
    sfree(inten);
    sfree(ev);
}

/*
 * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
 */
static real calc_eps(double M_diff, double volume, double epsRF, double temp)
{
    double eps, A, teller, noemer;
    double eps_0 = 8.854187817e-12;   /* epsilon_0 in C^2 J^-1 m^-1 */
    double fac   = 1.112650021e-59;   /* converts Debye^2 to C^2 m^2 */

    A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);

    if (epsRF == 0.0)
    {
        teller = 1 + A;
        noemer = 1;
    }
    else
    {
        teller = 1 + (A*2*epsRF/(2*epsRF+1));
        noemer = 1 - (A/(2*epsRF+1));
    }
    eps = teller / noemer;

    return eps;
}

static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
                                int idim, int nslice, rvec slab_dipole[],
                                matrix box)
{
    int  k;
    real xdim = 0;

    for (k = k0; (k < k1); k++)
    {
        xdim += x[k][idim];
    }
    xdim /= (k1-k0);
    k     = (static_cast<int>(xdim*nslice/box[idim][idim] + nslice)) % nslice;
    rvec_inc(slab_dipole[k], mu);
}

static void dump_slab_dipoles(const char *fn, int idim, int nslice,
                              rvec slab_dipole[], matrix box, int nframes,
                              const gmx_output_env_t *oenv)
{
    FILE       *fp;
    char        buf[STRLEN];
    int         i;
    real        mutot;
    const char *leg_dim[4] = {
        "\\f{12}m\\f{4}\\sX\\N",
        "\\f{12}m\\f{4}\\sY\\N",
        "\\f{12}m\\f{4}\\sZ\\N",
        "\\f{12}m\\f{4}\\stot\\N"
    };

    sprintf(buf, "Box-%c (nm)", 'X'+idim);
    fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
                  oenv);
    xvgr_legend(fp, DIM, leg_dim, oenv);
    for (i = 0; (i < nslice); i++)
    {
        mutot = norm(slab_dipole[i])/nframes;
        fprintf(fp, "%10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
                ((i+0.5)*box[idim][idim])/nslice,
                slab_dipole[i][XX]/nframes,
                slab_dipole[i][YY]/nframes,
                slab_dipole[i][ZZ]/nframes,
                mutot);
    }
    xvgrclose(fp);
    do_view(oenv, fn, "-autoscale xy -nxy");
}

static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
{
    int    i, j, k;
    double dc, d, ddc1, ddc2, ddc3;
    rvec   xxx = { 1, 0, 0 };
    rvec   yyy = { 0, 1, 0 };
    rvec   zzz = { 0, 0, 1 };

    d    = 0;
    ddc1 = ddc2 = ddc3 = 0;
    for (i = k = 0; (i < n); i++)
    {
        ddc1 += std::abs(cos_angle(dip[i], xxx));
        ddc2 += std::abs(cos_angle(dip[i], yyy));
        ddc3 += std::abs(cos_angle(dip[i], zzz));
        if (bPairs)
        {
            for (j = i+1; (j < n); j++, k++)
            {
                dc  = cos_angle(dip[i], dip[j]);
                d  += std::abs(dc);
            }
        }
    }
    *dd      = d/k;
    axis[XX] = ddc1/n;
    axis[YY] = ddc2/n;
    axis[ZZ] = ddc3/n;
}

static void do_dip(const t_topology *top, int ePBC, real volume,
                   const char *fn,
                   const char *out_mtot, const char *out_eps,
                   const char *out_aver, const char *dipdist,
                   const char *cosaver, const char *fndip3d,
                   const char *fnadip,  gmx_bool bPairs,
                   const char *corrtype, const char *corf,
                   gmx_bool bGkr,     const char *gkrfn,
                   gmx_bool bPhi,     int  *nlevels,  int ndegrees,
                   int  ncos,
                   const char *cmap,    real rcmax,
                   gmx_bool bQuad,
                   gmx_bool bMU,      const char *mufn,
                   int  *gnx,     int  *molindex[],
                   real mu_max,   real mu_aver,
                   real epsilonRF, real temp,
                   int  *gkatom,  int skip,
                   gmx_bool bSlab,    int nslices,
                   const char *axtitle, const char *slabfn,
                   const gmx_output_env_t *oenv)
{
    const char *leg_mtot[] = {
        "M\\sx \\N",
        "M\\sy \\N",
        "M\\sz \\N",
        "|M\\stot \\N|"
    };
#define NLEGMTOT asize(leg_mtot)
    const char *leg_eps[] = {
        "epsilon",
        "G\\sk",
        "g\\sk"
    };
#define NLEGEPS asize(leg_eps)
    const char *leg_aver[] = {
        "< |M|\\S2\\N >",
        "< |M| >\\S2\\N",
        "< |M|\\S2\\N > - < |M| >\\S2\\N",
        "< |M| >\\S2\\N / < |M|\\S2\\N >"
    };
#define NLEGAVER asize(leg_aver)
    const char *leg_cosaver[] = {
        "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
        "RMSD cos",
        "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
        "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
        "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
    };
#define NLEGCOSAVER asize(leg_cosaver)
    const char *leg_adip[] = {
        "<mu>",
        "Std. Dev.",
        "Error"
    };
#define NLEGADIP asize(leg_adip)

    FILE          *outdd, *outmtot, *outaver, *outeps, *caver = nullptr;
    FILE          *dip3d = nullptr, *adip = nullptr;
    rvec          *x, *dipole = nullptr, mu_t, quad, *dipsp = nullptr;
    t_gkrbin      *gkrbin = nullptr;
    gmx_enxnm_t   *enm    = nullptr;
    t_enxframe    *fr;
    int            nframes = 1000, nre, timecheck = 0, ncolour = 0;
    ener_file_t    fmu     = nullptr;
    int            i, n, m, natom = 0, gnx_tot, teller, tel3;
    t_trxstatus   *status;
    int           *dipole_bin, ndipbin, ibin, iVol, idim = -1;
    unsigned long  mode;
    real           rcut = 0, t, t0, t1, dt, dd, rms_cos;
    rvec           dipaxis;
    matrix         box;
    gmx_bool       bCorr, bTotal, bCont;
    double         M_diff = 0, epsilon, invtel, vol_aver;
    double         mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
    double         M[3], M2[3], M4[3], Gk = 0, g_k = 0;
    gmx_stats_t   *Qlsq, mulsq, muframelsq = nullptr;
    ivec           iMu;
    real         **muall        = nullptr;
    rvec          *slab_dipoles = nullptr;
    const t_atom  *atom         = nullptr;
    const t_block *mols         = nullptr;
    gmx_rmpbc_t    gpbc         = nullptr;

    gnx_tot = gnx[0];
    if (ncos == 2)
    {
        gnx_tot += gnx[1];
    }
    GMX_RELEASE_ASSERT(ncos == 1 || ncos == 2, "Invalid number of groups used with -ncos");

    vol_aver = 0.0;

    iVol = -1;

    /* initialize to a negative value so we can see that it wasn't picked up */
    iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
    if (bMU)
    {
        fmu = open_enx(mufn, "r");
        do_enxnms(fmu, &nre, &enm);

        /* Determine the indexes of the energy grps we need */
        for (i = 0; (i < nre); i++)
        {
            if (std::strstr(enm[i].name, "Volume"))
            {
                iVol = i;
            }
            else if (std::strstr(enm[i].name, "Mu-X"))
            {
                iMu[XX] = i;
            }
            else if (std::strstr(enm[i].name, "Mu-Y"))
            {
                iMu[YY] = i;
            }
            else if (std::strstr(enm[i].name, "Mu-Z"))
            {
                iMu[ZZ] = i;
            }
        }
        if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
        {
            gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
        }
    }
    else
    {
        atom = top->atoms.atom;
        mols = &(top->mols);
    }
    if ((iVol == -1) && bMU)
    {
        printf("Using Volume from topology: %g nm^3\n", volume);
    }

    /* Correlation stuff */
    bCorr  = (corrtype[0] != 'n');
    bTotal = (corrtype[0] == 't');
    if (bCorr)
    {
        if (bTotal)
        {
            snew(muall, 1);
            snew(muall[0], nframes*DIM);
        }
        else
        {
            snew(muall, gnx[0]);
            for (i = 0; (i < gnx[0]); i++)
            {
                snew(muall[i], nframes*DIM);
            }
        }
    }

    /* Allocate array which contains for every molecule in a frame the
     * dipole moment.
     */
    if (!bMU)
    {
        snew(dipole, gnx_tot);
    }

    /* Statistics */
    snew(Qlsq, DIM);
    for (i = 0; (i < DIM); i++)
    {
        Qlsq[i] = gmx_stats_init();
    }
    mulsq = gmx_stats_init();

    /* Open all the files */
    outmtot = xvgropen(out_mtot,
                       "Total dipole moment of the simulation box vs. time",
                       "Time (ps)", "Total Dipole Moment (Debye)", oenv);
    outeps  = xvgropen(out_eps, "Epsilon and Kirkwood factors",
                       "Time (ps)", "", oenv);
    outaver = xvgropen(out_aver, "Total dipole moment",
                       "Time (ps)", "D", oenv);
    if (bSlab)
    {
        idim = axtitle[0] - 'X';
        if ((idim < 0) || (idim >= DIM))
        {
            idim = axtitle[0] - 'x';
        }
        if ((idim < 0) || (idim >= DIM))
        {
            bSlab = FALSE;
        }
        if (nslices < 2)
        {
            bSlab = FALSE;
        }
        fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
                axtitle, nslices, idim);
        if (bSlab)
        {
            snew(slab_dipoles, nslices);
            fprintf(stderr, "Doing slab analysis\n");
        }
    }

    if (fnadip)
    {
        adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
        xvgr_legend(adip, NLEGADIP, leg_adip, oenv);

    }
    if (cosaver)
    {
        caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
                         "Average absolute dipole orientation", "Time (ps)", "", oenv);
        xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
                    oenv);
    }

    if (fndip3d)
    {
        snew(dipsp, gnx_tot);

        /* we need a dummy file for gnuplot */
        dip3d = gmx_ffopen("dummy.dat", "w");
        fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
        gmx_ffclose(dip3d);

        dip3d = gmx_ffopen(fndip3d, "w");
        try
        {
            gmx::BinaryInformationSettings settings;
            settings.generatedByHeader(true);
            settings.linePrefix("# ");
            gmx::printBinaryInformation(dip3d, output_env_get_program_context(oenv),
                                        settings);
        }
        GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
    }

    /* Write legends to all the files */
    xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
    xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);

    if (bMU && (mu_aver == -1))
    {
        xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
    }
    else
    {
        xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
    }

    snew(fr, 1);
    clear_rvec(mu_t);
    teller = 0;
    /* Read the first frame from energy or traj file */
    if (bMU)
    {
        do
        {
            bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
            if (bCont)
            {
                timecheck = check_times(t);
                if (timecheck < 0)
                {
                    teller++;
                }
                if ((teller % 10) == 0)
                {
                    fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
                    fflush(stderr);
                }
            }
            else
            {
                printf("End of %s reached\n", mufn);
                break;
            }
        }
        while (bCont && (timecheck < 0));
    }
    else
    {
        natom  = read_first_x(oenv, &status, fn, &t, &x, box);
    }

    /* Calculate spacing for dipole bin (simple histogram) */
    ndipbin = 1 + static_cast<int>(mu_max/0.01);
    snew(dipole_bin, ndipbin);
    mu_ave     = 0.0;
    for (m = 0; (m < DIM); m++)
    {
        M[m] = M2[m] = M4[m] = 0.0;
    }

    if (bGkr)
    {
        /* Use 0.7 iso 0.5 to account for pressure scaling */
        /*  rcut   = 0.7*sqrt(max_cutoff2(box)); */
        rcut   = 0.7*std::sqrt(gmx::square(box[XX][XX])+gmx::square(box[YY][YY])+gmx::square(box[ZZ][ZZ]));

        gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
    }
    gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);

    /* Start while loop over frames */
    t0     = t;
    teller = 0;
    do
    {
        if (bCorr && (teller >= nframes))
        {
            nframes += 1000;
            if (bTotal)
            {
                srenew(muall[0], nframes*DIM);
            }
            else
            {
                for (i = 0; (i < gnx_tot); i++)
                {
                    srenew(muall[i], nframes*DIM);
                }
            }
        }
        t1 = t;

        muframelsq = gmx_stats_init();

        /* Initialise */
        for (m = 0; (m < DIM); m++)
        {
            M_av2[m] = 0;
        }

        if (bMU)
        {
            /* Copy rvec into double precision local variable */
            for (m = 0; (m < DIM); m++)
            {
                M_av[m]  = mu_t[m];
            }
        }
        else
        {
            /* Initialise */
            for (m = 0; (m < DIM); m++)
            {
                M_av[m] = 0;
            }

            gmx_rmpbc(gpbc, natom, box, x);

            /* Begin loop of all molecules in frame */
            for (n = 0; (n < ncos); n++)
            {
                for (i = 0; (i < gnx[n]); i++)
                {
                    int ind0, ind1;

                    ind0  = mols->index[molindex[n][i]];
                    ind1  = mols->index[molindex[n][i]+1];

                    mol_dip(ind0, ind1, x, atom, dipole[i]);
                    gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
                    gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
                    if (bSlab)
                    {
                        update_slab_dipoles(ind0, ind1, x,
                                            dipole[i], idim, nslices, slab_dipoles, box);
                    }
                    if (bQuad)
                    {
                        mol_quad(ind0, ind1, x, atom, quad);
                        for (m = 0; (m < DIM); m++)
                        {
                            gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
                        }
                    }
                    if (bCorr && !bTotal)
                    {
                        tel3              = DIM*teller;
                        muall[i][tel3+XX] = dipole[i][XX];
                        muall[i][tel3+YY] = dipole[i][YY];
                        muall[i][tel3+ZZ] = dipole[i][ZZ];
                    }
                    mu_mol = 0.0;
                    for (m = 0; (m < DIM); m++)
                    {
                        M_av[m]  += dipole[i][m];               /* M per frame */
                        mu_mol   += dipole[i][m]*dipole[i][m];  /* calc. mu for distribution */
                    }
                    mu_mol = std::sqrt(mu_mol);

                    mu_ave += mu_mol;                         /* calc. the average mu */

                    /* Update the dipole distribution */
                    ibin = gmx::roundToInt(ndipbin*mu_mol/mu_max);
                    if (ibin < ndipbin)
                    {
                        dipole_bin[ibin]++;
                    }

                    if (fndip3d)
                    {
                        rvec2sprvec(dipole[i], dipsp[i]);

                        if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
                        {
                            if (dipsp[i][ZZ] < 0.5 * M_PI)
                            {
                                ncolour = 1;
                            }
                            else
                            {
                                ncolour = 2;
                            }
                        }
                        else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
                        {
                            if (dipsp[i][ZZ] < 0.5 * M_PI)
                            {
                                ncolour = 3;
                            }
                            else
                            {
                                ncolour = 4;
                            }
                        }
                        else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
                        {
                            if (dipsp[i][ZZ] < 0.5 * M_PI)
                            {
                                ncolour = 5;
                            }
                            else
                            {
                                ncolour = 6;
                            }
                        }
                        else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
                        {
                            if (dipsp[i][ZZ] < 0.5 * M_PI)
                            {
                                ncolour = 7;
                            }
                            else
                            {
                                ncolour = 8;
                            }
                        }
                        if (dip3d)
                        {
                            fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d  # %d %d\n",
                                    i+1,
                                    x[ind0][XX],
                                    x[ind0][YY],
                                    x[ind0][ZZ],
                                    x[ind0][XX]+dipole[i][XX]/25,
                                    x[ind0][YY]+dipole[i][YY]/25,
                                    x[ind0][ZZ]+dipole[i][ZZ]/25,
                                    ncolour, ind0, i);
                        }
                    }
                } /* End loop of all molecules in frame */

                if (dip3d)
                {
                    fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
                    fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
                    fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
                    fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
                    fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
                    fprintf(dip3d, "pause -1 'Hit return to continue'\n");
                }
            }
        }
        /* Compute square of total dipole */
        for (m = 0; (m < DIM); m++)
        {
            M_av2[m] = M_av[m]*M_av[m];
        }

        if (cosaver)
        {
            compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
            rms_cos = std::sqrt(gmx::square(dipaxis[XX]-0.5)+
                                gmx::square(dipaxis[YY]-0.5)+
                                gmx::square(dipaxis[ZZ]-0.5));
            if (bPairs)
            {
                fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
                        t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
            }
            else
            {
                fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
                        t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
            }
        }

        if (bGkr)
        {
            do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
                   atom, gkatom);
        }

        if (bTotal)
        {
            tel3              = DIM*teller;
            muall[0][tel3+XX] = M_av[XX];
            muall[0][tel3+YY] = M_av[YY];
            muall[0][tel3+ZZ] = M_av[ZZ];
        }

        /* Write to file the total dipole moment of the box, and its components
         * for this frame.
         */
        if ((skip == 0) || ((teller % skip) == 0))
        {
            fprintf(outmtot, "%10g  %12.8e %12.8e %12.8e %12.8e\n",
                    t, M_av[XX], M_av[YY], M_av[ZZ],
                    std::sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
        }

        for (m = 0; (m < DIM); m++)
        {
            M[m]  += M_av[m];
            M2[m] += M_av2[m];
            M4[m] += gmx::square(M_av2[m]);
        }
        /* Increment loop counter */
        teller++;

        /* Calculate for output the running averages */
        invtel  = 1.0/teller;
        M2_ave  = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
        M_ave2  = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
        M_diff  = M2_ave - M_ave2;

        /* Compute volume from box in traj, else we use the one from above */
        if (!bMU)
        {
            volume  = det(box);
        }
        vol_aver += volume;

        epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);

        /* Calculate running average for dipole */
        if (mu_ave != 0)
        {
            mu_aver = (mu_ave/gnx_tot)*invtel;
        }

        if ((skip == 0) || ((teller % skip) == 0))
        {
            /* Write to file < |M|^2 >, |< M >|^2. And the difference between
             * the two. Here M is sum mu_i. Further write the finite system
             * Kirkwood G factor and epsilon.
             */
            fprintf(outaver, "%10g  %10.3e %10.3e %10.3e %10.3e\n",
                    t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);

            if (fnadip)
            {
                real aver;
                gmx_stats_get_average(muframelsq, &aver);
                fprintf(adip, "%10g %f \n", t, aver);
            }
            /*if (dipole)
               printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
             */
            if (!bMU || (mu_aver != -1))
            {
                /* Finite system Kirkwood G-factor */
                Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
                /* Infinite system Kirkwood G-factor */
                if (epsilonRF == 0.0)
                {
                    g_k = ((2*epsilon+1)*Gk/(3*epsilon));
                }
                else
                {
                    g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
                           Gk/(3*epsilon*(2*epsilonRF+1)));
                }

                fprintf(outeps, "%10g  %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);

            }
            else
            {
                fprintf(outeps, "%10g  %12.8e\n", t, epsilon);
            }
        }
        gmx_stats_free(muframelsq);

        if (bMU)
        {
            bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
        }
        else
        {
            bCont = read_next_x(oenv, status, &t, x, box);
        }
        timecheck = check_times(t);
    }
    while (bCont && (timecheck == 0) );

    gmx_rmpbc_done(gpbc);

    if (!bMU)
    {
        close_trx(status);
    }

    xvgrclose(outmtot);
    xvgrclose(outaver);
    xvgrclose(outeps);

    if (fnadip)
    {
        xvgrclose(adip);
    }

    if (cosaver)
    {
        xvgrclose(caver);
    }

    if (dip3d)
    {
        fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
        fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
        fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
        fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
        fprintf(dip3d, "pause -1 'Hit return to continue'\n");
        gmx_ffclose(dip3d);
    }

    if (bSlab)
    {
        dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
        sfree(slab_dipoles);
    }

    vol_aver /= teller;
    printf("Average volume over run is %g\n", vol_aver);
    if (bGkr)
    {
        print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
        print_cmap(cmap, gkrbin, nlevels);
    }
    /* Autocorrelation function */
    if (bCorr)
    {
        if (teller < 2)
        {
            printf("Not enough frames for autocorrelation\n");
        }
        else
        {
            dt = (t1 - t0)/(teller-1);
            printf("t0 %g, t %g, teller %d\n", t0, t, teller);

            mode = eacVector;

            if (bTotal)
            {
                do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
                            teller, 1, muall, dt, mode, TRUE);
            }
            else
            {
                do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
                            teller, gnx_tot, muall, dt,
                            mode, std::strcmp(corrtype, "molsep") != 0);
            }
        }
    }
    if (!bMU)
    {
        real aver, sigma, error;

        gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
        printf("\nDipole moment (Debye)\n");
        printf("---------------------\n");
        printf("Average  = %8.4f  Std. Dev. = %8.4f  Error = %8.4f\n",
               aver, sigma, error);
        if (bQuad)
        {
            rvec a, s, e;
            for (m = 0; (m < DIM); m++)
            {
                gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
            }

            printf("\nQuadrupole moment (Debye-Ang)\n");
            printf("-----------------------------\n");
            printf("Averages  = %8.4f  %8.4f  %8.4f\n", a[XX], a[YY], a[ZZ]);
            printf("Std. Dev. = %8.4f  %8.4f  %8.4f\n", s[XX], s[YY], s[ZZ]);
            printf("Error     = %8.4f  %8.4f  %8.4f\n", e[XX], e[YY], e[ZZ]);
        }
        printf("\n");
    }
    printf("The following averages for the complete trajectory have been calculated:\n\n");
    printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
    printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
    printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);

    printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
    printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
    printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);

    printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
    printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);

    printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);

    if (!bMU || (mu_aver != -1))
    {
        printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
        printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
    }
    printf("Epsilon = %g\n", epsilon);

    if (!bMU)
    {
        /* Write to file the dipole moment distibution during the simulation.
         */
        outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
        for (i = 0; (i < ndipbin); i++)
        {
            fprintf(outdd, "%10g  %10f\n",
                    (i*mu_max)/ndipbin, static_cast<real>(dipole_bin[i])/teller);
        }
        xvgrclose(outdd);
        sfree(dipole_bin);
    }
    if (bGkr)
    {
        done_gkrbin(&gkrbin);
    }
}

static void dipole_atom2molindex(int *n, int *index, const t_block *mols)
{
    int nmol, i, j, m;

    nmol = 0;
    i    = 0;
    while (i < *n)
    {
        m = 0;
        while (m < mols->nr && index[i] != mols->index[m])
        {
            m++;
        }
        if (m == mols->nr)
        {
            gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
        }
        for (j = mols->index[m]; j < mols->index[m+1]; j++)
        {
            if (i >= *n || index[i] != j)
            {
                gmx_fatal(FARGS, "The index group is not a set of whole molecules");
            }
            i++;
        }
        /* Modify the index in place */
        index[nmol++] = m;
    }
    printf("There are %d molecules in the selection\n", nmol);

    *n = nmol;
}
int gmx_dipoles(int argc, char *argv[])
{
    const char       *desc[] = {
        "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
        "system. From this you can compute e.g. the dielectric constant for",
        "low-dielectric media.",
        "For molecules with a net charge, the net charge is subtracted at",
        "center of mass of the molecule.[PAR]",
        "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
        "components as well as the norm of the vector.",
        "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",
        "simulation.",
        "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
        "the simulation",
        "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
        "Furthermore, the dipole autocorrelation function will be computed when",
        "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
        "option.",
        "The correlation functions can be averaged over all molecules",
        "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
        "or it can be computed over the total dipole moment of the simulation box",
        "([TT]total[tt]).[PAR]",
        "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
        "G-factor, as well as the average cosine of the angle between the dipoles",
        "as a function of the distance. The plot also includes gOO and hOO",
        "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
        "we also include the energy per scale computed by taking the inner product of",
        "the dipoles divided by the distance to the third power.[PAR]",
        "[PAR]",
        "EXAMPLES[PAR]",
        "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
        "This will calculate the autocorrelation function of the molecular",
        "dipoles using a first order Legendre polynomial of the angle of the",
        "dipole vector and itself a time t later. For this calculation 1001",
        "frames will be used. Further, the dielectric constant will be calculated",
        "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
        "an average dipole moment of the molecule of 2.273 (SPC). For the",
        "distribution function a maximum of 5.0 will be used."
    };
    real              mu_max     = 5, mu_aver = -1, rcmax = 0;
    real              epsilonRF  = 0.0, temp = 300;
    gmx_bool          bPairs     = TRUE, bPhi = FALSE, bQuad = FALSE;
    const char       *corrtype[] = {nullptr, "none", "mol", "molsep", "total", nullptr};
    const char       *axtitle    = "Z";
    int               nslices    = 10; /* nr of slices defined       */
    int               skip       = 0, nFA = 0, nFB = 0, ncos = 1;
    int               nlevels    = 20, ndegrees = 90;
    gmx_output_env_t *oenv;
    t_pargs           pa[] = {
        { "-mu",       FALSE, etREAL, {&mu_aver},
          "dipole of a single molecule (in Debye)" },
        { "-mumax",    FALSE, etREAL, {&mu_max},
          "max dipole in Debye (for histogram)" },
        { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
          "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
        { "-skip",     FALSE, etINT, {&skip},
          "Skip steps in the output (but not in the computations)" },
        { "-temp",     FALSE, etREAL, {&temp},
          "Average temperature of the simulation (needed for dielectric constant calculation)" },
        { "-corr",     FALSE, etENUM, {corrtype},
          "Correlation function to calculate" },
        { "-pairs",    FALSE, etBOOL, {&bPairs},
          "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
        { "-quad",     FALSE, etBOOL, {&bQuad},
          "Take quadrupole into account"},
        { "-ncos",     FALSE, etINT, {&ncos},
          "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." },
        { "-axis",     FALSE, etSTR, {&axtitle},
          "Take the normal on the computational box in direction X, Y or Z." },
        { "-sl",       FALSE, etINT, {&nslices},
          "Divide the box into this number of slices." },
        { "-gkratom",  FALSE, etINT, {&nFA},
          "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" },
        { "-gkratom2", FALSE, etINT, {&nFB},
          "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
        { "-rcmax",    FALSE, etREAL, {&rcmax},
          "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
        { "-phi",      FALSE, etBOOL, {&bPhi},
          "Plot the 'torsion angle' defined as the rotation of the two dipole vectors around the distance vector between the two molecules in the [REF].xpm[ref] file from the [TT]-cmap[tt] option. By default the cosine of the angle between the dipoles is plotted." },
        { "-nlevels",  FALSE, etINT, {&nlevels},
          "Number of colors in the cmap output" },
        { "-ndegrees", FALSE, etINT, {&ndegrees},
          "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
    };
    int              *gnx;
    int               nFF[2];
    int             **grpindex;
    char            **grpname = nullptr;
    gmx_bool          bGkr, bMU, bSlab;
    t_filenm          fnm[] = {
        { efEDR, "-en", nullptr,         ffOPTRD },
        { efTRX, "-f", nullptr,           ffREAD },
        { efTPR, nullptr, nullptr,           ffREAD },
        { efNDX, nullptr, nullptr,           ffOPTRD },
        { efXVG, "-o",   "Mtot",       ffWRITE },
        { efXVG, "-eps", "epsilon",    ffWRITE },
        { efXVG, "-a",   "aver",       ffWRITE },
        { efXVG, "-d",   "dipdist",    ffWRITE },
        { efXVG, "-c",   "dipcorr",    ffOPTWR },
        { efXVG, "-g",   "gkr",        ffOPTWR },
        { efXVG, "-adip", "adip",       ffOPTWR },
        { efXVG, "-dip3d", "dip3d",    ffOPTWR },
        { efXVG, "-cos", "cosaver",    ffOPTWR },
        { efXPM, "-cmap", "cmap",       ffOPTWR },
        { efXVG, "-slab", "slab",       ffOPTWR }
    };
#define NFILE asize(fnm)
    int               npargs;
    t_pargs          *ppa;
    t_topology       *top;
    int               ePBC;
    int               k, natoms;
    matrix            box;

    npargs = asize(pa);
    ppa    = add_acf_pargs(&npargs, pa);
    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
                           NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
    {
        sfree(ppa);
        return 0;
    }

    printf("Using %g as mu_max and %g as the dipole moment.\n",
           mu_max, mu_aver);
    if (epsilonRF == 0.0)
    {
        printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
    }

    bMU   = opt2bSet("-en", NFILE, fnm);
    if (bMU)
    {
        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.");
    }
    bGkr  = opt2bSet("-g", NFILE, fnm);
    if (opt2parg_bSet("-ncos", asize(pa), pa))
    {
        if ((ncos != 1) && (ncos != 2))
        {
            gmx_fatal(FARGS, "ncos has to be either 1 or 2");
        }
        bGkr = TRUE;
    }
    bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
             opt2parg_bSet("-axis", asize(pa), pa));
    if (bMU)
    {
        if (bQuad)
        {
            printf("WARNING: Can not determine quadrupoles from energy file\n");
            bQuad = FALSE;
        }
        if (bGkr)
        {
            printf("WARNING: Can not determine Gk(r) from energy file\n");
            bGkr  = FALSE;
            ncos  = 1;
        }
        if (mu_aver == -1)
        {
            printf("WARNING: Can not calculate Gk and gk, since you did\n"
                   "         not enter a valid dipole for the molecules\n");
        }
    }

    snew(top, 1);
    ePBC = read_tpx_top(ftp2fn(efTPR, NFILE, fnm), nullptr, box,
                        &natoms, nullptr, nullptr, top);

    snew(gnx, ncos);
    snew(grpname, ncos);
    snew(grpindex, ncos);
    get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
              ncos, gnx, grpindex, grpname);
    for (k = 0; (k < ncos); k++)
    {
        dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
        neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
    }
    nFF[0] = nFA;
    nFF[1] = nFB;
    do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
           opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
           opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
           opt2fn_null("-cos", NFILE, fnm),
           opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
           bPairs, corrtype[0],
           opt2fn("-c", NFILE, fnm),
           bGkr,    opt2fn("-g", NFILE, fnm),
           bPhi,    &nlevels,  ndegrees,
           ncos,
           opt2fn("-cmap", NFILE, fnm), rcmax,
           bQuad, bMU,     opt2fn("-en", NFILE, fnm),
           gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
           bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);

    do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
    do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
    do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
    do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
    do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");

    return 0;
}