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