src/gromacs/gmxana/gmx_spol.c

   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, 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 "gromacs/commandline/pargs.h"
  40 #include "gromacs/fileio/tpxio.h"
  41 #include "gromacs/fileio/trxio.h"
  42 #include "gromacs/fileio/xvgr.h"
  43 #include "gromacs/gmxana/gmx_ana.h"
  44 #include "gromacs/gmxana/gstat.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/legacyheaders/viewit.h"
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
  56                          atom_id index[], rvec xref, int ePBC)
  57 {
  58     const real tol = 1e-4;
  59     gmx_bool   bChanged;
  60     int        m, j, ai, iter;
  61     real       mass, mtot;
  62     rvec       dx, xtest;
  63
  64     /* First simple calculation */
  65     clear_rvec(xref);
  66     mtot = 0;
  67     for (m = 0; (m < nrefat); m++)
  68     {
  69         ai   = index[m];
  70         mass = top->atoms.atom[ai].m;
  71         for (j = 0; (j < DIM); j++)
  72         {
  73             xref[j] += mass*x[ai][j];
  74         }
  75         mtot += mass;
  76     }
  77     svmul(1/mtot, xref, xref);
  78     /* Now check if any atom is more than half the box from the COM */
  79     if (ePBC != epbcNONE)
  80     {
  81         iter = 0;
  82         do
  83         {
  84             bChanged = FALSE;
  85             for (m = 0; (m < nrefat); m++)
  86             {
  87                 ai   = index[m];
  88                 mass = top->atoms.atom[ai].m/mtot;
  89                 pbc_dx(pbc, x[ai], xref, dx);
  90                 rvec_add(xref, dx, xtest);
  91                 for (j = 0; (j < DIM); j++)
  92                 {
  93                     if (fabs(xtest[j]-x[ai][j]) > tol)
  94                     {
  95                         /* Here we have used the wrong image for contributing to the COM */
  96                         xref[j] += mass*(xtest[j]-x[ai][j]);
  97                         x[ai][j] = xtest[j];
  98                         bChanged = TRUE;
  99                     }
 100                 }
 101             }
 102             if (bChanged)
 103             {
 104                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 105             }
 106             iter++;
 107         }
 108         while (bChanged);
 109     }
 110 }
 111
 112 void spol_atom2molindex(int *n, int *index, t_block *mols)
 113 {
 114     int nmol, i, j, m;
 115
 116     nmol = 0;
 117     i    = 0;
 118     while (i < *n)
 119     {
 120         m = 0;
 121         while (m < mols->nr && index[i] != mols->index[m])
 122         {
 123             m++;
 124         }
 125         if (m == mols->nr)
 126         {
 127             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
 128         }
 129         for (j = mols->index[m]; j < mols->index[m+1]; j++)
 130         {
 131             if (i >= *n || index[i] != j)
 132             {
 133                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
 134             }
 135             i++;
 136         }
 137         /* Modify the index in place */
 138         index[nmol++] = m;
 139     }
 140     printf("There are %d molecules in the selection\n", nmol);
 141
 142     *n = nmol;
 143 }
 144
 145 int gmx_spol(int argc, char *argv[])
 146 {
 147     t_topology  *top;
 148     t_inputrec  *ir;
 149     t_atom      *atom;
 150     char         title[STRLEN];
 151     t_trxstatus *status;
 152     int          nrefat, natoms, nf, ntot;
 153     real         t;
 154     rvec        *xtop, *x, xref, trial, dx = {0}, dip, dir;
 155     matrix       box;
 156
 157     FILE        *fp;
 158     int         *isize, nrefgrp;
 159     atom_id    **index, *molindex;
 160     char       **grpname;
 161     real         rmin2, rmax2, rcut, rcut2, rdx2 = 0, rtry2, qav, q, dip2, invbw;
 162     int          nbin, i, m, mol, a0, a1, a, d;
 163     double       sdip, sdip2, sinp, sdinp, nmol;
 164     int         *hist;
 165     t_pbc        pbc;
 166     gmx_rmpbc_t  gpbc = NULL;
 167
 168
 169     const char     *desc[] = {
 170         "[THISMODULE] analyzes dipoles around a solute; it is especially useful",
 171         "for polarizable water. A group of reference atoms, or a center",
 172         "of mass reference (option [TT]-com[tt]) and a group of solvent",
 173         "atoms is required. The program splits the group of solvent atoms",
 174         "into molecules. For each solvent molecule the distance to the",
 175         "closest atom in reference group or to the COM is determined.",
 176         "A cumulative distribution of these distances is plotted.",
 177         "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
 178         "the inner product of the distance vector",
 179         "and the dipole of the solvent molecule is determined.",
 180         "For solvent molecules with net charge (ions), the net charge of the ion",
 181         "is subtracted evenly from all atoms in the selection of each ion.",
 182         "The average of these dipole components is printed.",
 183         "The same is done for the polarization, where the average dipole is",
 184         "subtracted from the instantaneous dipole. The magnitude of the average",
 185         "dipole is set with the option [TT]-dip[tt], the direction is defined",
 186         "by the vector from the first atom in the selected solvent group",
 187         "to the midpoint between the second and the third atom."
 188     };
 189
 190     output_env_t    oenv;
 191     static gmx_bool bCom   = FALSE, bPBC = FALSE;
 192     static int      srefat = 1;
 193     static real     rmin   = 0.0, rmax = 0.32, refdip = 0, bw = 0.01;
 194     t_pargs         pa[]   = {
 195         { "-com",  FALSE, etBOOL,  {&bCom},
 196           "Use the center of mass as the reference postion" },
 197         { "-refat",  FALSE, etINT, {&srefat},
 198           "The reference atom of the solvent molecule" },
 199         { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
 200         { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 201         { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
 202         { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
 203     };
 204
 205     t_filenm        fnm[] = {
 206         { efTRX, NULL,  NULL,  ffREAD },
 207         { efTPR, NULL,  NULL,  ffREAD },
 208         { efNDX, NULL,  NULL,  ffOPTRD },
 209         { efXVG, NULL,  "scdist",  ffWRITE }
 210     };
 211 #define NFILE asize(fnm)
 212
 213     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 214                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 215     {
 216         return 0;
 217     }
 218
 219     snew(top, 1);
 220     snew(ir, 1);
 221     read_tpx_top(ftp2fn(efTPR, NFILE, fnm),
 222                  ir, box, &natoms, NULL, NULL, NULL, 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*sqrt(max_cutoff2(ir->ePBC, box));
 249     if (rcut == 0)
 250     {
 251         rcut = 10*rmax;
 252     }
 253     rcut2 = sqr(rcut);
 254     invbw = 1/bw;
 255     nbin  = (int)(rcut*invbw)+2;
 256     snew(hist, nbin);
 257
 258     rmin2 = sqr(rmin);
 259     rmax2 = sqr(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[(int)(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  += 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     }
 351     while (read_next_x(oenv, status, &t, x, box));
 352
 353     gmx_rmpbc_done(gpbc);
 354
 355     /* clean up */
 356     sfree(x);
 357     close_trj(status);
 358
 359     fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
 360             rmax, (real)ntot/(real)nf);
 361     if (ntot > 0)
 362     {
 363         sdip  /= ntot;
 364         sdip2 /= ntot;
 365         sinp  /= ntot;
 366         sdinp /= ntot;
 367         fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n",
 368                 sdip, sqrt(sdip2-sqr(sdip)));
 369         fprintf(stderr, "Average radial component of the dipole:       %f (D)\n",
 370                 sinp);
 371         fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
 372                 sdinp);
 373     }
 374
 375     fp = xvgropen(opt2fn("-o", NFILE, fnm),
 376                   "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     gmx_ffclose(fp);
 384
 385     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 386
 387     return 0;
 388 }