src/gromacs/gmxana/gmx_spol.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/commandline/viewit.h"
  43 #include "gromacs/fileio/tpxio.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/gmxana/gstat.h"
  48 #include "gromacs/math/functions.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdtypes/inputrec.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/pbcutil/rmpbc.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 static void calc_com_pbc(int nrefat, const t_topology* top, rvec x[], t_pbc* pbc, const int index[], rvec xref, int ePBC)
  62 {
  63     const real tol = 1e-4;
  64     gmx_bool   bChanged;
  65     int        m, j, ai, iter;
  66     real       mass, mtot;
  67     rvec       dx, xtest;
  68
  69     /* First simple calculation */
  70     clear_rvec(xref);
  71     mtot = 0;
  72     for (m = 0; (m < nrefat); m++)
  73     {
  74         ai   = index[m];
  75         mass = top->atoms.atom[ai].m;
  76         for (j = 0; (j < DIM); j++)
  77         {
  78             xref[j] += mass * x[ai][j];
  79         }
  80         mtot += mass;
  81     }
  82     svmul(1 / mtot, xref, xref);
  83     /* Now check if any atom is more than half the box from the COM */
  84     if (ePBC != epbcNONE)
  85     {
  86         iter = 0;
  87         do
  88         {
  89             bChanged = FALSE;
  90             for (m = 0; (m < nrefat); m++)
  91             {
  92                 ai   = index[m];
  93                 mass = top->atoms.atom[ai].m / mtot;
  94                 pbc_dx(pbc, x[ai], xref, dx);
  95                 rvec_add(xref, dx, xtest);
  96                 for (j = 0; (j < DIM); j++)
  97                 {
  98                     if (std::abs(xtest[j] - x[ai][j]) > tol)
  99                     {
 100                         /* Here we have used the wrong image for contributing to the COM */
 101                         xref[j] += mass * (xtest[j] - x[ai][j]);
 102                         x[ai][j] = xtest[j];
 103                         bChanged = TRUE;
 104                     }
 105                 }
 106             }
 107             if (bChanged)
 108             {
 109                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 110             }
 111             iter++;
 112         } while (bChanged);
 113     }
 114 }
 115
 116 static void spol_atom2molindex(int* n, int* index, const t_block* mols)
 117 {
 118     int nmol, i, j, m;
 119
 120     nmol = 0;
 121     i    = 0;
 122     while (i < *n)
 123     {
 124         m = 0;
 125         while (m < mols->nr && index[i] != mols->index[m])
 126         {
 127             m++;
 128         }
 129         if (m == mols->nr)
 130         {
 131             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule",
 132                       i + 1, index[i] + 1);
 133         }
 134         for (j = mols->index[m]; j < mols->index[m + 1]; j++)
 135         {
 136             if (i >= *n || index[i] != j)
 137             {
 138                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
 139             }
 140             i++;
 141         }
 142         /* Modify the index in place */
 143         index[nmol++] = m;
 144     }
 145     printf("There are %d molecules in the selection\n", nmol);
 146
 147     *n = nmol;
 148 }
 149
 150 int gmx_spol(int argc, char* argv[])
 151 {
 152     t_topology*  top;
 153     t_atom*      atom;
 154     t_trxstatus* status;
 155     int          nrefat, natoms, nf, ntot;
 156     real         t;
 157     rvec *       x, xref, trial, dx = { 0 }, dip, dir;
 158     matrix       box;
 159
 160     FILE*       fp;
 161     int *       isize, nrefgrp;
 162     int **      index, *molindex;
 163     char**      grpname;
 164     real        rmin2, rmax2, rcut, rcut2, rdx2 = 0, rtry2, qav, q, dip2, invbw;
 165     int         nbin, i, m, mol, a0, a1, a, d;
 166     double      sdip, sdip2, sinp, sdinp, nmol;
 167     int*        hist;
 168     t_pbc       pbc;
 169     gmx_rmpbc_t gpbc = nullptr;
 170
 171
 172     const char* desc[] = {
 173         "[THISMODULE] analyzes dipoles around a solute; it is especially useful",
 174         "for polarizable water. A group of reference atoms, or a center",
 175         "of mass reference (option [TT]-com[tt]) and a group of solvent",
 176         "atoms is required. The program splits the group of solvent atoms",
 177         "into molecules. For each solvent molecule the distance to the",
 178         "closest atom in reference group or to the COM is determined.",
 179         "A cumulative distribution of these distances is plotted.",
 180         "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
 181         "the inner product of the distance vector",
 182         "and the dipole of the solvent molecule is determined.",
 183         "For solvent molecules with net charge (ions), the net charge of the ion",
 184         "is subtracted evenly from all atoms in the selection of each ion.",
 185         "The average of these dipole components is printed.",
 186         "The same is done for the polarization, where the average dipole is",
 187         "subtracted from the instantaneous dipole. The magnitude of the average",
 188         "dipole is set with the option [TT]-dip[tt], the direction is defined",
 189         "by the vector from the first atom in the selected solvent group",
 190         "to the midpoint between the second and the third atom."
 191     };
 192
 193     gmx_output_env_t* oenv;
 194     static gmx_bool   bCom   = FALSE;
 195     static int        srefat = 1;
 196     static real       rmin = 0.0, rmax = 0.32, refdip = 0, bw = 0.01;
 197     t_pargs           pa[] = {
 198         { "-com", FALSE, etBOOL, { &bCom }, "Use the center of mass as the reference position" },
 199         { "-refat", FALSE, etINT, { &srefat }, "The reference atom of the solvent molecule" },
 200         { "-rmin", FALSE, etREAL, { &rmin }, "Maximum distance (nm)" },
 201         { "-rmax", FALSE, etREAL, { &rmax }, "Maximum distance (nm)" },
 202         { "-dip", FALSE, etREAL, { &refdip }, "The average dipole (D)" },
 203         { "-bw", FALSE, etREAL, { &bw }, "The bin width" }
 204     };
 205
 206     t_filenm fnm[] = { { efTRX, nullptr, nullptr, ffREAD },
 207                        { efTPR, nullptr, nullptr, ffREAD },
 208                        { efNDX, nullptr, nullptr, ffOPTRD },
 209                        { efXVG, nullptr, "scdist", ffWRITE } };
 210 #define NFILE asize(fnm)
 211
 212     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
 213                            asize(desc), desc, 0, nullptr, &oenv))
 214     {
 215         return 0;
 216     }
 217
 218     snew(top, 1);
 219     // TODO: Only ePBC is used, not the full inputrec.
 220     t_inputrec  irInstance;
 221     t_inputrec* ir = &irInstance;
 222     read_tpx_top(ftp2fn(efTPR, NFILE, fnm), ir, box, &natoms, nullptr, nullptr, top);
 223
 224     /* get index groups */
 225     printf("Select a group of reference particles and a solvent group:\n");
 226     snew(grpname, 2);
 227     snew(index, 2);
 228     snew(isize, 2);
 229     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 230
 231     if (bCom)
 232     {
 233         nrefgrp = 1;
 234         nrefat  = isize[0];
 235     }
 236     else
 237     {
 238         nrefgrp = isize[0];
 239         nrefat  = 1;
 240     }
 241
 242     spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
 243     srefat--;
 244
 245     /* initialize reading trajectory:                         */
 246     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 247
 248     rcut = 0.99 * std::sqrt(max_cutoff2(ir->ePBC, box));
 249     if (rcut == 0)
 250     {
 251         rcut = 10 * rmax;
 252     }
 253     rcut2 = gmx::square(rcut);
 254     invbw = 1 / bw;
 255     nbin  = static_cast<int>(rcut * invbw) + 2;
 256     snew(hist, nbin);
 257
 258     rmin2 = gmx::square(rmin);
 259     rmax2 = gmx::square(rmax);
 260
 261     nf    = 0;
 262     ntot  = 0;
 263     sdip  = 0;
 264     sdip2 = 0;
 265     sinp  = 0;
 266     sdinp = 0;
 267
 268     molindex = top->mols.index;
 269     atom     = top->atoms.atom;
 270
 271     gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
 272
 273     /* start analysis of trajectory */
 274     do
 275     {
 276         /* make molecules whole again */
 277         gmx_rmpbc(gpbc, natoms, box, x);
 278
 279         set_pbc(&pbc, ir->ePBC, box);
 280         if (bCom)
 281         {
 282             calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC);
 283         }
 284
 285         for (m = 0; m < isize[1]; m++)
 286         {
 287             mol = index[1][m];
 288             a0  = molindex[mol];
 289             a1  = molindex[mol + 1];
 290             for (i = 0; i < nrefgrp; i++)
 291             {
 292                 pbc_dx(&pbc, x[a0 + srefat], bCom ? xref : x[index[0][i]], trial);
 293                 rtry2 = norm2(trial);
 294                 if (i == 0 || rtry2 < rdx2)
 295                 {
 296                     copy_rvec(trial, dx);
 297                     rdx2 = rtry2;
 298                 }
 299             }
 300             if (rdx2 < rcut2)
 301             {
 302                 hist[static_cast<int>(std::sqrt(rdx2) * invbw) + 1]++;
 303             }
 304             if (rdx2 >= rmin2 && rdx2 < rmax2)
 305             {
 306                 unitv(dx, dx);
 307                 clear_rvec(dip);
 308                 qav = 0;
 309                 for (a = a0; a < a1; a++)
 310                 {
 311                     qav += atom[a].q;
 312                 }
 313                 qav /= (a1 - a0);
 314                 for (a = a0; a < a1; a++)
 315                 {
 316                     q = atom[a].q - qav;
 317                     for (d = 0; d < DIM; d++)
 318                     {
 319                         dip[d] += q * x[a][d];
 320                     }
 321                 }
 322                 for (d = 0; d < DIM; d++)
 323                 {
 324                     dir[d] = -x[a0][d];
 325                 }
 326                 for (a = a0 + 1; a < a0 + 3; a++)
 327                 {
 328                     for (d = 0; d < DIM; d++)
 329                     {
 330                         dir[d] += 0.5 * x[a][d];
 331                     }
 332                 }
 333                 unitv(dir, dir);
 334
 335                 svmul(ENM2DEBYE, dip, dip);
 336                 dip2 = norm2(dip);
 337                 sdip += std::sqrt(dip2);
 338                 sdip2 += dip2;
 339                 for (d = 0; d < DIM; d++)
 340                 {
 341                     sinp += dx[d] * dip[d];
 342                     sdinp += dx[d] * (dip[d] - refdip * dir[d]);
 343                 }
 344
 345                 ntot++;
 346             }
 347         }
 348         nf++;
 349
 350     } while (read_next_x(oenv, status, &t, x, box));
 351
 352     gmx_rmpbc_done(gpbc);
 353
 354     /* clean up */
 355     sfree(x);
 356     close_trx(status);
 357
 358     fprintf(stderr, "Average number of molecules within %g nm is %.1f\n", rmax,
 359             static_cast<real>(ntot) / nf);
 360     if (ntot > 0)
 361     {
 362         sdip /= ntot;
 363         sdip2 /= ntot;
 364         sinp /= ntot;
 365         sdinp /= ntot;
 366         fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n", sdip,
 367                 std::sqrt(sdip2 - gmx::square(sdip)));
 368         fprintf(stderr, "Average radial component of the dipole:       %f (D)\n", sinp);
 369         fprintf(stderr, "Average radial component of the polarization: %f (D)\n", sdinp);
 370     }
 371
 372     fp   = xvgropen(opt2fn("-o", NFILE, fnm), "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
 373     nmol = 0;
 374     for (i = 0; i <= nbin; i++)
 375     {
 376         nmol += hist[i];
 377         fprintf(fp, "%g %g\n", i * bw, nmol / nf);
 378     }
 379     xvgrclose(fp);
 380
 381     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 382
 383     return 0;
 384 }