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 "config.h"
  40
  41 #include "gromacs/legacyheaders/macros.h"
  42 #include "gstat.h"
  43 #include "gromacs/legacyheaders/viewit.h"
  44 #include "gromacs/pbcutil/pbc.h"
  45 #include "gromacs/topology/index.h"
  46 #include "gromacs/fileio/tpxio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/math/units.h"
  49 #include "gmx_ana.h"
  50
  51 #include "gromacs/commandline/pargs.h"
  52 #include "gromacs/fileio/xvgr.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/pbcutil/rmpbc.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
  59                          atom_id index[], rvec xref, int ePBC)
  60 {
  61     const real tol = 1e-4;
  62     gmx_bool   bChanged;
  63     int        m, j, ai, iter;
  64     real       mass, mtot;
  65     rvec       dx, xtest;
  66
  67     /* First simple calculation */
  68     clear_rvec(xref);
  69     mtot = 0;
  70     for (m = 0; (m < nrefat); m++)
  71     {
  72         ai   = index[m];
  73         mass = top->atoms.atom[ai].m;
  74         for (j = 0; (j < DIM); j++)
  75         {
  76             xref[j] += mass*x[ai][j];
  77         }
  78         mtot += mass;
  79     }
  80     svmul(1/mtot, xref, xref);
  81     /* Now check if any atom is more than half the box from the COM */
  82     if (ePBC != epbcNONE)
  83     {
  84         iter = 0;
  85         do
  86         {
  87             bChanged = FALSE;
  88             for (m = 0; (m < nrefat); m++)
  89             {
  90                 ai   = index[m];
  91                 mass = top->atoms.atom[ai].m/mtot;
  92                 pbc_dx(pbc, x[ai], xref, dx);
  93                 rvec_add(xref, dx, xtest);
  94                 for (j = 0; (j < DIM); j++)
  95                 {
  96                     if (fabs(xtest[j]-x[ai][j]) > tol)
  97                     {
  98                         /* Here we have used the wrong image for contributing to the COM */
  99                         xref[j] += mass*(xtest[j]-x[ai][j]);
 100                         x[ai][j] = xtest[j];
 101                         bChanged = TRUE;
 102                     }
 103                 }
 104             }
 105             if (bChanged)
 106             {
 107                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 108             }
 109             iter++;
 110         }
 111         while (bChanged);
 112     }
 113 }
 114
 115 void spol_atom2molindex(int *n, int *index, t_block *mols)
 116 {
 117     int nmol, i, j, m;
 118
 119     nmol = 0;
 120     i    = 0;
 121     while (i < *n)
 122     {
 123         m = 0;
 124         while (m < mols->nr && index[i] != mols->index[m])
 125         {
 126             m++;
 127         }
 128         if (m == mols->nr)
 129         {
 130             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
 131         }
 132         for (j = mols->index[m]; j < mols->index[m+1]; j++)
 133         {
 134             if (i >= *n || index[i] != j)
 135             {
 136                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
 137             }
 138             i++;
 139         }
 140         /* Modify the index in place */
 141         index[nmol++] = m;
 142     }
 143     printf("There are %d molecules in the selection\n", nmol);
 144
 145     *n = nmol;
 146 }
 147
 148 int gmx_spol(int argc, char *argv[])
 149 {
 150     t_topology  *top;
 151     t_inputrec  *ir;
 152     t_atom      *atom;
 153     char         title[STRLEN];
 154     t_trxstatus *status;
 155     int          nrefat, natoms, nf, ntot;
 156     real         t;
 157     rvec        *xtop, *x, xref, trial, dx = {0}, dip, dir;
 158     matrix       box;
 159
 160     FILE        *fp;
 161     int         *isize, nrefgrp;
 162     atom_id    **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 = NULL;
 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     output_env_t    oenv;
 194     static gmx_bool bCom   = FALSE, bPBC = 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},
 199           "Use the center of mass as the reference postion" },
 200         { "-refat",  FALSE, etINT, {&srefat},
 201           "The reference atom of the solvent molecule" },
 202         { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
 203         { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 204         { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
 205         { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
 206     };
 207
 208     t_filenm        fnm[] = {
 209         { efTRX, NULL,  NULL,  ffREAD },
 210         { efTPX, NULL,  NULL,  ffREAD },
 211         { efNDX, NULL,  NULL,  ffOPTRD },
 212         { efXVG, NULL,  "scdist.xvg",  ffWRITE }
 213     };
 214 #define NFILE asize(fnm)
 215
 216     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 217                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 218     {
 219         return 0;
 220     }
 221
 222     snew(top, 1);
 223     snew(ir, 1);
 224     read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
 225                  ir, box, &natoms, NULL, NULL, NULL, 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*sqrt(max_cutoff2(ir->ePBC, box));
 252     if (rcut == 0)
 253     {
 254         rcut = 10*rmax;
 255     }
 256     rcut2 = sqr(rcut);
 257     invbw = 1/bw;
 258     nbin  = (int)(rcut*invbw)+2;
 259     snew(hist, nbin);
 260
 261     rmin2 = sqr(rmin);
 262     rmax2 = sqr(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->ePBC, 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->ePBC, box);
 283         if (bCom)
 284         {
 285             calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC);
 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[(int)(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  += 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     }
 354     while (read_next_x(oenv, status, &t, x, box));
 355
 356     gmx_rmpbc_done(gpbc);
 357
 358     /* clean up */
 359     sfree(x);
 360     close_trj(status);
 361
 362     fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
 363             rmax, (real)ntot/(real)nf);
 364     if (ntot > 0)
 365     {
 366         sdip  /= ntot;
 367         sdip2 /= ntot;
 368         sinp  /= ntot;
 369         sdinp /= ntot;
 370         fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n",
 371                 sdip, sqrt(sdip2-sqr(sdip)));
 372         fprintf(stderr, "Average radial component of the dipole:       %f (D)\n",
 373                 sinp);
 374         fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
 375                 sdinp);
 376     }
 377
 378     fp = xvgropen(opt2fn("-o", NFILE, fnm),
 379                   "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
 380     nmol = 0;
 381     for (i = 0; i <= nbin; i++)
 382     {
 383         nmol += hist[i];
 384         fprintf(fp, "%g %g\n", i*bw, nmol/nf);
 385     }
 386     gmx_ffclose(fp);
 387
 388     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 389
 390     return 0;
 391 }