src/gromacs/gmxana/gmx_spol.c

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