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 "copyrite.h"
  43 #include "gstat.h"
  44 #include "vec.h"
  45 #include "xvgr.h"
  46 #include "pbc.h"
  47 #include "index.h"
  48 #include "tpxio.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, matrix box)
  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         "[TT]g_spol[tt] 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     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     snew(top, 1);
 215     snew(ir, 1);
 216     read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
 217                  ir, box, &natoms, NULL, NULL, NULL, top);
 218
 219     /* get index groups */
 220     printf("Select a group of reference particles and a solvent group:\n");
 221     snew(grpname, 2);
 222     snew(index, 2);
 223     snew(isize, 2);
 224     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 225
 226     if (bCom)
 227     {
 228         nrefgrp = 1;
 229         nrefat  = isize[0];
 230     }
 231     else
 232     {
 233         nrefgrp = isize[0];
 234         nrefat  = 1;
 235     }
 236
 237     spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
 238     srefat--;
 239
 240     /* initialize reading trajectory:                         */
 241     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 242
 243     rcut  = 0.99*sqrt(max_cutoff2(ir->ePBC, box));
 244     if (rcut == 0)
 245     {
 246         rcut = 10*rmax;
 247     }
 248     rcut2 = sqr(rcut);
 249     invbw = 1/bw;
 250     nbin  = (int)(rcut*invbw)+2;
 251     snew(hist, nbin);
 252
 253     rmin2 = sqr(rmin);
 254     rmax2 = sqr(rmax);
 255
 256     nf    = 0;
 257     ntot  = 0;
 258     sdip  = 0;
 259     sdip2 = 0;
 260     sinp  = 0;
 261     sdinp = 0;
 262
 263     molindex = top->mols.index;
 264     atom     = top->atoms.atom;
 265
 266     gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms, box);
 267
 268     /* start analysis of trajectory */
 269     do
 270     {
 271         /* make molecules whole again */
 272         gmx_rmpbc(gpbc, natoms, box, x);
 273
 274         set_pbc(&pbc, ir->ePBC, box);
 275         if (bCom)
 276         {
 277             calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC, box);
 278         }
 279
 280         for (m = 0; m < isize[1]; m++)
 281         {
 282             mol = index[1][m];
 283             a0  = molindex[mol];
 284             a1  = molindex[mol+1];
 285             for (i = 0; i < nrefgrp; i++)
 286             {
 287                 pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
 288                 rtry2 = norm2(trial);
 289                 if (i == 0 || rtry2 < rdx2)
 290                 {
 291                     copy_rvec(trial, dx);
 292                     rdx2 = rtry2;
 293                 }
 294             }
 295             if (rdx2 < rcut2)
 296             {
 297                 hist[(int)(sqrt(rdx2)*invbw)+1]++;
 298             }
 299             if (rdx2 >= rmin2 && rdx2 < rmax2)
 300             {
 301                 unitv(dx, dx);
 302                 clear_rvec(dip);
 303                 qav = 0;
 304                 for (a = a0; a < a1; a++)
 305                 {
 306                     qav += atom[a].q;
 307                 }
 308                 qav /= (a1 - a0);
 309                 for (a = a0; a < a1; a++)
 310                 {
 311                     q = atom[a].q - qav;
 312                     for (d = 0; d < DIM; d++)
 313                     {
 314                         dip[d] += q*x[a][d];
 315                     }
 316                 }
 317                 for (d = 0; d < DIM; d++)
 318                 {
 319                     dir[d] = -x[a0][d];
 320                 }
 321                 for (a = a0+1; a < a0+3; a++)
 322                 {
 323                     for (d = 0; d < DIM; d++)
 324                     {
 325                         dir[d] += 0.5*x[a][d];
 326                     }
 327                 }
 328                 unitv(dir, dir);
 329
 330                 svmul(ENM2DEBYE, dip, dip);
 331                 dip2   = norm2(dip);
 332                 sdip  += sqrt(dip2);
 333                 sdip2 += dip2;
 334                 for (d = 0; d < DIM; d++)
 335                 {
 336                     sinp  += dx[d]*dip[d];
 337                     sdinp += dx[d]*(dip[d] - refdip*dir[d]);
 338                 }
 339
 340                 ntot++;
 341             }
 342         }
 343         nf++;
 344
 345     }
 346     while (read_next_x(oenv, status, &t, natoms, x, box));
 347
 348     gmx_rmpbc_done(gpbc);
 349
 350     /* clean up */
 351     sfree(x);
 352     close_trj(status);
 353
 354     fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
 355             rmax, (real)ntot/(real)nf);
 356     if (ntot > 0)
 357     {
 358         sdip  /= ntot;
 359         sdip2 /= ntot;
 360         sinp  /= ntot;
 361         sdinp /= ntot;
 362         fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n",
 363                 sdip, sqrt(sdip2-sqr(sdip)));
 364         fprintf(stderr, "Average radial component of the dipole:       %f (D)\n",
 365                 sinp);
 366         fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
 367                 sdinp);
 368     }
 369
 370     fp = xvgropen(opt2fn("-o", NFILE, fnm),
 371                   "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
 372     nmol = 0;
 373     for (i = 0; i <= nbin; i++)
 374     {
 375         nmol += hist[i];
 376         fprintf(fp, "%g %g\n", i*bw, nmol/nf);
 377     }
 378     ffclose(fp);
 379
 380     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 381
 382     thanx(stderr);
 383
 384     return 0;
 385 }