src/gromacs/gmxana/gmx_sorient.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, 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/confio.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/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/pbcutil/rmpbc.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
  60                          int index[], rvec xref, gmx_bool bPBC)
  61 {
  62     const real tol = 1e-4;
  63     gmx_bool   bChanged;
  64     int        m, j, ai, iter;
  65     real       mass, mtot;
  66     rvec       dx, xtest;
  67
  68     /* First simple calculation */
  69     clear_rvec(xref);
  70     mtot = 0;
  71     for (m = 0; (m < nrefat); m++)
  72     {
  73         ai   = index[m];
  74         mass = top->atoms.atom[ai].m;
  75         for (j = 0; (j < DIM); j++)
  76         {
  77             xref[j] += mass*x[ai][j];
  78         }
  79         mtot += mass;
  80     }
  81     svmul(1/mtot, xref, xref);
  82     /* Now check if any atom is more than half the box from the COM */
  83     if (bPBC)
  84     {
  85         iter = 0;
  86         do
  87         {
  88             bChanged = FALSE;
  89             for (m = 0; (m < nrefat); m++)
  90             {
  91                 ai   = index[m];
  92                 mass = top->atoms.atom[ai].m/mtot;
  93                 pbc_dx(pbc, x[ai], xref, dx);
  94                 rvec_add(xref, dx, xtest);
  95                 for (j = 0; (j < DIM); j++)
  96                 {
  97                     if (std::abs(xtest[j]-x[ai][j]) > tol)
  98                     {
  99                         /* Here we have used the wrong image for contributing to the COM */
 100                         xref[j] += mass*(xtest[j]-x[ai][j]);
 101                         x[ai][j] = xtest[j];
 102                         bChanged = TRUE;
 103                     }
 104                 }
 105             }
 106             if (bChanged)
 107             {
 108                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 109             }
 110             iter++;
 111         }
 112         while (bChanged);
 113     }
 114 }
 115
 116 int gmx_sorient(int argc, char *argv[])
 117 {
 118     t_topology        top;
 119     int               ePBC = -1;
 120     t_trxstatus      *status;
 121     int               natoms;
 122     real              t;
 123     rvec             *xtop, *x;
 124     matrix            box;
 125
 126     FILE             *fp;
 127     int               i, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
 128     real             *histi1, *histi2, invbw, invrbw;
 129     double            sum1, sum2;
 130     int              *isize, nrefgrp, nrefat;
 131     int             **index;
 132     char            **grpname;
 133     real              inp, outp, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r;
 134     real              c1, c2;
 135     char              str[STRLEN];
 136     gmx_bool          bTPS;
 137     rvec              xref, dx, dxh1, dxh2, outer;
 138     gmx_rmpbc_t       gpbc = NULL;
 139     t_pbc             pbc;
 140     const char       *legr[] = {
 141         "<cos(\\8q\\4\\s1\\N)>",
 142         "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
 143     };
 144     const char       *legc[] = {
 145         "cos(\\8q\\4\\s1\\N)",
 146         "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
 147     };
 148
 149     const char       *desc[] = {
 150         "[THISMODULE] analyzes solvent orientation around solutes.",
 151         "It calculates two angles between the vector from one or more",
 152         "reference positions to the first atom of each solvent molecule:",
 153         "",
 154         " * [GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
 155         "   molecule to the midpoint between atoms 2 and 3.",
 156         " * [GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
 157         "   same three atoms, or, when the option [TT]-v23[tt] is set, ",
 158         "   the angle with the vector between atoms 2 and 3.",
 159         "",
 160         "The reference can be a set of atoms or",
 161         "the center of mass of a set of atoms. The group of solvent atoms should",
 162         "consist of 3 atoms per solvent molecule.",
 163         "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
 164         "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
 165         "[TT]-o[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
 166         "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
 167         "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
 168         "distance.[PAR]",
 169         "[TT]-co[tt]: the sum over all solvent molecules within distance r",
 170         "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] as a function of r.[PAR]",
 171         "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
 172     };
 173
 174     gmx_output_env_t *oenv;
 175     static gmx_bool   bCom = FALSE, bVec23 = FALSE, bPBC = FALSE;
 176     static real       rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02;
 177     t_pargs           pa[] = {
 178         { "-com",  FALSE, etBOOL,  {&bCom},
 179           "Use the center of mass as the reference position" },
 180         { "-v23",  FALSE, etBOOL,  {&bVec23},
 181           "Use the vector between atoms 2 and 3" },
 182         { "-rmin",  FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
 183         { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 184         { "-cbin",  FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
 185         { "-rbin",  FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
 186         { "-pbc",   FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
 187     };
 188
 189     t_filenm          fnm[] = {
 190         { efTRX, NULL,  NULL,  ffREAD },
 191         { efTPS, NULL,  NULL,  ffREAD },
 192         { efNDX, NULL,  NULL,  ffOPTRD },
 193         { efXVG, NULL,  "sori",   ffWRITE },
 194         { efXVG, "-no", "snor",   ffWRITE },
 195         { efXVG, "-ro", "sord",   ffWRITE },
 196         { efXVG, "-co", "scum",   ffWRITE },
 197         { efXVG, "-rc", "scount", ffWRITE }
 198     };
 199 #define NFILE asize(fnm)
 200
 201     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 202                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 203     {
 204         return 0;
 205     }
 206
 207     bTPS = (opt2bSet("-s", NFILE, fnm) || !opt2bSet("-n", NFILE, fnm) || bCom);
 208     if (bTPS)
 209     {
 210         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, NULL, box,
 211                       bCom);
 212     }
 213
 214     /* get index groups */
 215     printf("Select a group of reference particles and a solvent group:\n");
 216     snew(grpname, 2);
 217     snew(index, 2);
 218     snew(isize, 2);
 219     if (bTPS)
 220     {
 221         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 222     }
 223     else
 224     {
 225         get_index(NULL, ftp2fn(efNDX, NFILE, fnm), 2, isize, index, grpname);
 226     }
 227
 228     if (bCom)
 229     {
 230         nrefgrp = 1;
 231         nrefat  = isize[0];
 232     }
 233     else
 234     {
 235         nrefgrp = isize[0];
 236         nrefat  = 1;
 237     }
 238
 239     if (isize[1] % 3)
 240     {
 241         gmx_fatal(FARGS, "The number of solvent atoms (%d) is not a multiple of 3",
 242                   isize[1]);
 243     }
 244
 245     /* initialize reading trajectory:                         */
 246     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 247
 248     rmin2 = gmx::square(rmin);
 249     rmax2 = gmx::square(rmax);
 250     rcut  = 0.99*std::sqrt(max_cutoff2(guess_ePBC(box), box));
 251     if (rcut == 0)
 252     {
 253         rcut = 10*rmax;
 254     }
 255     rcut2 = gmx::square(rcut);
 256
 257     invbw = 1/binwidth;
 258     nbin1 = 1+static_cast<int>(2*invbw + 0.5);
 259     nbin2 = 1+static_cast<int>(invbw + 0.5);
 260
 261     invrbw = 1/rbinw;
 262
 263     snew(hist1, nbin1);
 264     snew(hist2, nbin2);
 265     nrbin = 1+static_cast<int>(rcut/rbinw);
 266     if (nrbin == 0)
 267     {
 268         nrbin = 1;
 269     }
 270     snew(histi1, nrbin);
 271     snew(histi2, nrbin);
 272     snew(histn, nrbin);
 273
 274     ntot = 0;
 275     nf   = 0;
 276     sum1 = 0;
 277     sum2 = 0;
 278
 279     if (bTPS)
 280     {
 281         /* make molecules whole again */
 282         gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 283     }
 284     /* start analysis of trajectory */
 285     do
 286     {
 287         if (bTPS)
 288         {
 289             /* make molecules whole again */
 290             gmx_rmpbc(gpbc, natoms, box, x);
 291         }
 292
 293         set_pbc(&pbc, ePBC, box);
 294         n    = 0;
 295         inp  = 0;
 296         for (p = 0; (p < nrefgrp); p++)
 297         {
 298             if (bCom)
 299             {
 300                 calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC);
 301             }
 302             else
 303             {
 304                 copy_rvec(x[index[0][p]], xref);
 305             }
 306
 307             for (m = 0; m < isize[1]; m += 3)
 308             {
 309                 sa0 = index[1][m];
 310                 sa1 = index[1][m+1];
 311                 sa2 = index[1][m+2];
 312                 range_check(sa0, 0, natoms);
 313                 range_check(sa1, 0, natoms);
 314                 range_check(sa2, 0, natoms);
 315                 pbc_dx(&pbc, x[sa0], xref, dx);
 316                 r2  = norm2(dx);
 317                 if (r2 < rcut2)
 318                 {
 319                     r = std::sqrt(r2);
 320                     if (!bVec23)
 321                     {
 322                         /* Determine the normal to the plain */
 323                         rvec_sub(x[sa1], x[sa0], dxh1);
 324                         rvec_sub(x[sa2], x[sa0], dxh2);
 325                         rvec_inc(dxh1, dxh2);
 326                         svmul(1/r, dx, dx);
 327                         unitv(dxh1, dxh1);
 328                         inp = iprod(dx, dxh1);
 329                         cprod(dxh1, dxh2, outer);
 330                         unitv(outer, outer);
 331                         outp = iprod(dx, outer);
 332                     }
 333                     else
 334                     {
 335                         /* Use the vector between the 2nd and 3rd atom */
 336                         rvec_sub(x[sa2], x[sa1], dxh2);
 337                         unitv(dxh2, dxh2);
 338                         outp = iprod(dx, dxh2)/r;
 339                     }
 340                     {
 341                         int ii = static_cast<int>(invrbw*r);
 342                         range_check(ii, 0, nrbin);
 343                         histi1[ii] += inp;
 344                         histi2[ii] += 3*gmx::square(outp) - 1;
 345                         histn[ii]++;
 346                     }
 347                     if ((r2 >= rmin2) && (r2 < rmax2))
 348                     {
 349                         int ii1 = static_cast<int>(invbw*(inp + 1));
 350                         int ii2 = static_cast<int>(invbw*std::abs(outp));
 351
 352                         range_check(ii1, 0, nbin1);
 353                         range_check(ii2, 0, nbin2);
 354                         hist1[ii1]++;
 355                         hist2[ii2]++;
 356                         sum1 += inp;
 357                         sum2 += outp;
 358                         n++;
 359                     }
 360                 }
 361             }
 362         }
 363         ntot += n;
 364         nf++;
 365
 366     }
 367     while (read_next_x(oenv, status, &t, x, box));
 368
 369     /* clean up */
 370     sfree(x);
 371     close_trj(status);
 372     gmx_rmpbc_done(gpbc);
 373
 374     /* Add the bin for the exact maximum to the previous bin */
 375     hist1[nbin1-1] += hist1[nbin1];
 376     hist2[nbin2-1] += hist2[nbin2];
 377
 378     nav     = static_cast<real>(ntot)/(nrefgrp*nf);
 379     normfac = invbw/ntot;
 380
 381     fprintf(stderr,  "Average nr of molecules between %g and %g nm: %.1f\n",
 382             rmin, rmax, nav);
 383     if (ntot > 0)
 384     {
 385         sum1 /= ntot;
 386         sum2 /= ntot;
 387         fprintf(stderr, "Average cos(theta1)     between %g and %g nm: %6.3f\n",
 388                 rmin, rmax, sum1);
 389         fprintf(stderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
 390                 rmin, rmax, sum2);
 391     }
 392
 393     sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax);
 394     fp = xvgropen(opt2fn("-o", NFILE, fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv);
 395     if (output_env_get_print_xvgr_codes(oenv))
 396     {
 397         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 398     }
 399     for (i = 0; i < nbin1; i++)
 400     {
 401         fprintf(fp, "%g %g\n", (i+0.5)*binwidth-1, 2*normfac*hist1[i]);
 402     }
 403     xvgrclose(fp);
 404
 405     sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax);
 406     fp = xvgropen(opt2fn("-no", NFILE, fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv);
 407     if (output_env_get_print_xvgr_codes(oenv))
 408     {
 409         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 410     }
 411     for (i = 0; i < nbin2; i++)
 412     {
 413         fprintf(fp, "%g %g\n", (i+0.5)*binwidth, normfac*hist2[i]);
 414     }
 415     xvgrclose(fp);
 416
 417
 418     sprintf(str, "Solvent orientation");
 419     fp = xvgropen(opt2fn("-ro", NFILE, fnm), str, "r (nm)", "", oenv);
 420     if (output_env_get_print_xvgr_codes(oenv))
 421     {
 422         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 423     }
 424     xvgr_legend(fp, 2, legr, oenv);
 425     for (i = 0; i < nrbin; i++)
 426     {
 427         fprintf(fp, "%g %g %g\n", (i+0.5)*rbinw,
 428                 histn[i] ? histi1[i]/histn[i] : 0,
 429                 histn[i] ? histi2[i]/histn[i] : 0);
 430     }
 431     xvgrclose(fp);
 432
 433     sprintf(str, "Cumulative solvent orientation");
 434     fp = xvgropen(opt2fn("-co", NFILE, fnm), str, "r (nm)", "", oenv);
 435     if (output_env_get_print_xvgr_codes(oenv))
 436     {
 437         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 438     }
 439     xvgr_legend(fp, 2, legc, oenv);
 440     normfac = 1.0/(nrefgrp*nf);
 441     c1      = 0;
 442     c2      = 0;
 443     fprintf(fp, "%g %g %g\n", 0.0, c1, c2);
 444     for (i = 0; i < nrbin; i++)
 445     {
 446         c1 += histi1[i]*normfac;
 447         c2 += histi2[i]*normfac;
 448         fprintf(fp, "%g %g %g\n", (i+1)*rbinw, c1, c2);
 449     }
 450     xvgrclose(fp);
 451
 452     sprintf(str, "Solvent distribution");
 453     fp = xvgropen(opt2fn("-rc", NFILE, fnm), str, "r (nm)", "molecules/nm", oenv);
 454     if (output_env_get_print_xvgr_codes(oenv))
 455     {
 456         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 457     }
 458     normfac = 1.0/(rbinw*nf);
 459     for (i = 0; i < nrbin; i++)
 460     {
 461         fprintf(fp, "%g %g\n", (i+0.5)*rbinw, histn[i]*normfac);
 462     }
 463     xvgrclose(fp);
 464
 465     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 466     do_view(oenv, opt2fn("-no", NFILE, fnm), NULL);
 467     do_view(oenv, opt2fn("-ro", NFILE, fnm), "-nxy");
 468     do_view(oenv, opt2fn("-co", NFILE, fnm), "-nxy");
 469
 470     return 0;
 471 }