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, 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/mdrunutility/mdmodules.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,
  63                          int index[], rvec xref, int ePBC)
  64 {
  65     const real tol = 1e-4;
  66     gmx_bool   bChanged;
  67     int        m, j, ai, iter;
  68     real       mass, mtot;
  69     rvec       dx, xtest;
  70
  71     /* First simple calculation */
  72     clear_rvec(xref);
  73     mtot = 0;
  74     for (m = 0; (m < nrefat); m++)
  75     {
  76         ai   = index[m];
  77         mass = top->atoms.atom[ai].m;
  78         for (j = 0; (j < DIM); j++)
  79         {
  80             xref[j] += mass*x[ai][j];
  81         }
  82         mtot += mass;
  83     }
  84     svmul(1/mtot, xref, xref);
  85     /* Now check if any atom is more than half the box from the COM */
  86     if (ePBC != epbcNONE)
  87     {
  88         iter = 0;
  89         do
  90         {
  91             bChanged = FALSE;
  92             for (m = 0; (m < nrefat); m++)
  93             {
  94                 ai   = index[m];
  95                 mass = top->atoms.atom[ai].m/mtot;
  96                 pbc_dx(pbc, x[ai], xref, dx);
  97                 rvec_add(xref, dx, xtest);
  98                 for (j = 0; (j < DIM); j++)
  99                 {
 100                     if (std::abs(xtest[j]-x[ai][j]) > tol)
 101                     {
 102                         /* Here we have used the wrong image for contributing to the COM */
 103                         xref[j] += mass*(xtest[j]-x[ai][j]);
 104                         x[ai][j] = xtest[j];
 105                         bChanged = TRUE;
 106                     }
 107                 }
 108             }
 109             if (bChanged)
 110             {
 111                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 112             }
 113             iter++;
 114         }
 115         while (bChanged);
 116     }
 117 }
 118
 119 void spol_atom2molindex(int *n, int *index, const t_block *mols)
 120 {
 121     int nmol, i, j, m;
 122
 123     nmol = 0;
 124     i    = 0;
 125     while (i < *n)
 126     {
 127         m = 0;
 128         while (m < mols->nr && index[i] != mols->index[m])
 129         {
 130             m++;
 131         }
 132         if (m == mols->nr)
 133         {
 134             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
 135         }
 136         for (j = mols->index[m]; j < mols->index[m+1]; j++)
 137         {
 138             if (i >= *n || index[i] != j)
 139             {
 140                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
 141             }
 142             i++;
 143         }
 144         /* Modify the index in place */
 145         index[nmol++] = m;
 146     }
 147     printf("There are %d molecules in the selection\n", nmol);
 148
 149     *n = nmol;
 150 }
 151
 152 int gmx_spol(int argc, char *argv[])
 153 {
 154     t_topology  *top;
 155     t_inputrec  *ir;
 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},
 202           "Use the center of mass as the reference position" },
 203         { "-refat",  FALSE, etINT, {&srefat},
 204           "The reference atom of the solvent molecule" },
 205         { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
 206         { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 207         { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
 208         { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
 209     };
 210
 211     t_filenm          fnm[] = {
 212         { efTRX, nullptr,  nullptr,  ffREAD },
 213         { efTPR, nullptr,  nullptr,  ffREAD },
 214         { efNDX, nullptr,  nullptr,  ffOPTRD },
 215         { efXVG, nullptr,  "scdist",  ffWRITE }
 216     };
 217 #define NFILE asize(fnm)
 218
 219     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 220                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 221     {
 222         return 0;
 223     }
 224
 225     snew(top, 1);
 226     // TODO: Only ePBC is used, not the full inputrec.
 227     gmx::MDModules mdModules;
 228     ir = mdModules.inputrec();
 229     read_tpx_top(ftp2fn(efTPR, NFILE, fnm),
 230                  ir, box, &natoms, nullptr, nullptr, top);
 231
 232     /* get index groups */
 233     printf("Select a group of reference particles and a solvent group:\n");
 234     snew(grpname, 2);
 235     snew(index, 2);
 236     snew(isize, 2);
 237     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 238
 239     if (bCom)
 240     {
 241         nrefgrp = 1;
 242         nrefat  = isize[0];
 243     }
 244     else
 245     {
 246         nrefgrp = isize[0];
 247         nrefat  = 1;
 248     }
 249
 250     spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
 251     srefat--;
 252
 253     /* initialize reading trajectory:                         */
 254     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 255
 256     rcut  = 0.99*std::sqrt(max_cutoff2(ir->ePBC, box));
 257     if (rcut == 0)
 258     {
 259         rcut = 10*rmax;
 260     }
 261     rcut2 = gmx::square(rcut);
 262     invbw = 1/bw;
 263     nbin  = static_cast<int>(rcut*invbw)+2;
 264     snew(hist, nbin);
 265
 266     rmin2 = gmx::square(rmin);
 267     rmax2 = gmx::square(rmax);
 268
 269     nf    = 0;
 270     ntot  = 0;
 271     sdip  = 0;
 272     sdip2 = 0;
 273     sinp  = 0;
 274     sdinp = 0;
 275
 276     molindex = top->mols.index;
 277     atom     = top->atoms.atom;
 278
 279     gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
 280
 281     /* start analysis of trajectory */
 282     do
 283     {
 284         /* make molecules whole again */
 285         gmx_rmpbc(gpbc, natoms, box, x);
 286
 287         set_pbc(&pbc, ir->ePBC, box);
 288         if (bCom)
 289         {
 290             calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC);
 291         }
 292
 293         for (m = 0; m < isize[1]; m++)
 294         {
 295             mol = index[1][m];
 296             a0  = molindex[mol];
 297             a1  = molindex[mol+1];
 298             for (i = 0; i < nrefgrp; i++)
 299             {
 300                 pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
 301                 rtry2 = norm2(trial);
 302                 if (i == 0 || rtry2 < rdx2)
 303                 {
 304                     copy_rvec(trial, dx);
 305                     rdx2 = rtry2;
 306                 }
 307             }
 308             if (rdx2 < rcut2)
 309             {
 310                 hist[static_cast<int>(std::sqrt(rdx2)*invbw)+1]++;
 311             }
 312             if (rdx2 >= rmin2 && rdx2 < rmax2)
 313             {
 314                 unitv(dx, dx);
 315                 clear_rvec(dip);
 316                 qav = 0;
 317                 for (a = a0; a < a1; a++)
 318                 {
 319                     qav += atom[a].q;
 320                 }
 321                 qav /= (a1 - a0);
 322                 for (a = a0; a < a1; a++)
 323                 {
 324                     q = atom[a].q - qav;
 325                     for (d = 0; d < DIM; d++)
 326                     {
 327                         dip[d] += q*x[a][d];
 328                     }
 329                 }
 330                 for (d = 0; d < DIM; d++)
 331                 {
 332                     dir[d] = -x[a0][d];
 333                 }
 334                 for (a = a0+1; a < a0+3; a++)
 335                 {
 336                     for (d = 0; d < DIM; d++)
 337                     {
 338                         dir[d] += 0.5*x[a][d];
 339                     }
 340                 }
 341                 unitv(dir, dir);
 342
 343                 svmul(ENM2DEBYE, dip, dip);
 344                 dip2   = norm2(dip);
 345                 sdip  += std::sqrt(dip2);
 346                 sdip2 += dip2;
 347                 for (d = 0; d < DIM; d++)
 348                 {
 349                     sinp  += dx[d]*dip[d];
 350                     sdinp += dx[d]*(dip[d] - refdip*dir[d]);
 351                 }
 352
 353                 ntot++;
 354             }
 355         }
 356         nf++;
 357
 358     }
 359     while (read_next_x(oenv, status, &t, x, box));
 360
 361     gmx_rmpbc_done(gpbc);
 362
 363     /* clean up */
 364     sfree(x);
 365     close_trj(status);
 366
 367     fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
 368             rmax, static_cast<real>(ntot)/nf);
 369     if (ntot > 0)
 370     {
 371         sdip  /= ntot;
 372         sdip2 /= ntot;
 373         sinp  /= ntot;
 374         sdinp /= ntot;
 375         fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n",
 376                 sdip, std::sqrt(sdip2-gmx::square(sdip)));
 377         fprintf(stderr, "Average radial component of the dipole:       %f (D)\n",
 378                 sinp);
 379         fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
 380                 sdinp);
 381     }
 382
 383     fp = xvgropen(opt2fn("-o", NFILE, fnm),
 384                   "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
 385     nmol = 0;
 386     for (i = 0; i <= nbin; i++)
 387     {
 388         nmol += hist[i];
 389         fprintf(fp, "%g %g\n", i*bw, nmol/nf);
 390     }
 391     xvgrclose(fp);
 392
 393     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 394
 395     return 0;
 396 }