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,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/commandline/viewit.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxana/gstat.h"
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/pbcutil/rmpbc.h"
  53 #include "gromacs/topology/index.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 static void calc_com_pbc(int nrefat, t_topology* top, rvec x[], t_pbc* pbc, const 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         } while (bChanged);
 112     }
 113 }
 114
 115 int gmx_sorient(int argc, char* argv[])
 116 {
 117     t_topology   top;
 118     PbcType      pbcType = PbcType::Unset;
 119     t_trxstatus* status;
 120     int          natoms;
 121     real         t;
 122     rvec *       xtop, *x;
 123     matrix       box;
 124
 125     FILE*       fp;
 126     int         i, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
 127     real *      histi1, *histi2, invbw, invrbw;
 128     double      sum1, sum2;
 129     int *       isize, nrefgrp, nrefat;
 130     int**       index;
 131     char**      grpname;
 132     real        inp, outp, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r;
 133     real        c1, c2;
 134     char        str[STRLEN];
 135     gmx_bool    bTPS;
 136     rvec        xref, dx, dxh1, dxh2, outer;
 137     gmx_rmpbc_t gpbc = nullptr;
 138     t_pbc       pbc;
 139     const char* legr[] = { "<cos(\\8q\\4\\s1\\N)>", "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>" };
 140     const char* legc[] = { "cos(\\8q\\4\\s1\\N)", "3cos\\S2\\N(\\8q\\4\\s2\\N)-1" };
 141
 142     const char* desc[] = {
 143         "[THISMODULE] analyzes solvent orientation around solutes.",
 144         "It calculates two angles between the vector from one or more",
 145         "reference positions to the first atom of each solvent molecule:",
 146         "",
 147         " * [GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of "
 148         "the solvent",
 149         "   molecule to the midpoint between atoms 2 and 3.",
 150         " * [GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, "
 151         "defined by the",
 152         "   same three atoms, or, when the option [TT]-v23[tt] is set, ",
 153         "   the angle with the vector between atoms 2 and 3.",
 154         "",
 155         "The reference can be a set of atoms or",
 156         "the center of mass of a set of atoms. The group of solvent atoms should",
 157         "consist of 3 atoms per solvent molecule.",
 158         "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
 159         "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
 160         "[TT]-o[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for "
 161         "rmin<=r<=rmax.[PAR]",
 162         "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for "
 163         "rmin<=r<=rmax.[PAR]",
 164         "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] "
 165         "and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a "
 166         "function of the",
 167         "distance.[PAR]",
 168         "[TT]-co[tt]: the sum over all solvent molecules within distance r",
 169         "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and "
 170         "[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 }, "Use the center of mass as the reference position" },
 179         { "-v23", FALSE, etBOOL, { &bVec23 }, "Use the vector between atoms 2 and 3" },
 180         { "-rmin", FALSE, etREAL, { &rmin }, "Minimum distance (nm)" },
 181         { "-rmax", FALSE, etREAL, { &rmax }, "Maximum distance (nm)" },
 182         { "-cbin", FALSE, etREAL, { &binwidth }, "Binwidth for the cosine" },
 183         { "-rbin", FALSE, etREAL, { &rbinw }, "Binwidth for r (nm)" },
 184         { "-pbc",
 185           FALSE,
 186           etBOOL,
 187           { &bPBC },
 188           "Check PBC for the center of mass calculation. Only necessary when your reference group "
 189           "consists of several molecules." }
 190     };
 191
 192     t_filenm fnm[] = { { efTRX, nullptr, nullptr, ffREAD },  { efTPS, nullptr, nullptr, ffREAD },
 193                        { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, nullptr, "sori", ffWRITE },
 194                        { efXVG, "-no", "snor", ffWRITE },    { efXVG, "-ro", "sord", ffWRITE },
 195                        { efXVG, "-co", "scum", ffWRITE },    { efXVG, "-rc", "scount", ffWRITE } };
 196 #define NFILE asize(fnm)
 197
 198     if (!parse_common_args(
 199                 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 200     {
 201         return 0;
 202     }
 203
 204     bTPS = (opt2bSet("-s", NFILE, fnm) || !opt2bSet("-n", NFILE, fnm) || bCom);
 205     if (bTPS)
 206     {
 207         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, box, bCom);
 208     }
 209
 210     /* get index groups */
 211     printf("Select a group of reference particles and a solvent group:\n");
 212     snew(grpname, 2);
 213     snew(index, 2);
 214     snew(isize, 2);
 215     if (bTPS)
 216     {
 217         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 218     }
 219     else
 220     {
 221         get_index(nullptr, ftp2fn(efNDX, NFILE, fnm), 2, isize, index, grpname);
 222     }
 223
 224     if (bCom)
 225     {
 226         nrefgrp = 1;
 227         nrefat  = isize[0];
 228     }
 229     else
 230     {
 231         nrefgrp = isize[0];
 232         nrefat  = 1;
 233     }
 234
 235     if (isize[1] % 3)
 236     {
 237         gmx_fatal(FARGS, "The number of solvent atoms (%d) is not a multiple of 3", isize[1]);
 238     }
 239
 240     /* initialize reading trajectory:                         */
 241     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 242
 243     rmin2 = gmx::square(rmin);
 244     rmax2 = gmx::square(rmax);
 245     rcut  = 0.99 * std::sqrt(max_cutoff2(guessPbcType(box), box));
 246     if (rcut == 0)
 247     {
 248         rcut = 10 * rmax;
 249     }
 250     rcut2 = gmx::square(rcut);
 251
 252     invbw = 1 / binwidth;
 253     nbin1 = 1 + gmx::roundToInt(2 * invbw);
 254     nbin2 = 1 + gmx::roundToInt(invbw);
 255
 256     invrbw = 1 / rbinw;
 257
 258     snew(hist1, nbin1);
 259     snew(hist2, nbin2);
 260     nrbin = 1 + static_cast<int>(rcut / rbinw);
 261     if (nrbin == 0)
 262     {
 263         nrbin = 1;
 264     }
 265     snew(histi1, nrbin);
 266     snew(histi2, nrbin);
 267     snew(histn, nrbin);
 268
 269     ntot = 0;
 270     nf   = 0;
 271     sum1 = 0;
 272     sum2 = 0;
 273
 274     if (bTPS)
 275     {
 276         /* make molecules whole again */
 277         gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
 278     }
 279     /* start analysis of trajectory */
 280     do
 281     {
 282         if (bTPS)
 283         {
 284             /* make molecules whole again */
 285             gmx_rmpbc(gpbc, natoms, box, x);
 286         }
 287
 288         set_pbc(&pbc, pbcType, box);
 289         n   = 0;
 290         inp = 0;
 291         for (p = 0; (p < nrefgrp); p++)
 292         {
 293             if (bCom)
 294             {
 295                 calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC);
 296             }
 297             else
 298             {
 299                 copy_rvec(x[index[0][p]], xref);
 300             }
 301
 302             for (m = 0; m < isize[1]; m += 3)
 303             {
 304                 sa0 = index[1][m];
 305                 sa1 = index[1][m + 1];
 306                 sa2 = index[1][m + 2];
 307                 range_check(sa0, 0, natoms);
 308                 range_check(sa1, 0, natoms);
 309                 range_check(sa2, 0, natoms);
 310                 pbc_dx(&pbc, x[sa0], xref, dx);
 311                 r2 = norm2(dx);
 312                 if (r2 < rcut2)
 313                 {
 314                     r = std::sqrt(r2);
 315                     if (!bVec23)
 316                     {
 317                         /* Determine the normal to the plain */
 318                         rvec_sub(x[sa1], x[sa0], dxh1);
 319                         rvec_sub(x[sa2], x[sa0], dxh2);
 320                         rvec_inc(dxh1, dxh2);
 321                         svmul(1 / r, dx, dx);
 322                         unitv(dxh1, dxh1);
 323                         inp = iprod(dx, dxh1);
 324                         cprod(dxh1, dxh2, outer);
 325                         unitv(outer, outer);
 326                         outp = iprod(dx, outer);
 327                     }
 328                     else
 329                     {
 330                         /* Use the vector between the 2nd and 3rd atom */
 331                         rvec_sub(x[sa2], x[sa1], dxh2);
 332                         unitv(dxh2, dxh2);
 333                         outp = iprod(dx, dxh2) / r;
 334                     }
 335                     {
 336                         int ii = static_cast<int>(invrbw * r);
 337                         range_check(ii, 0, nrbin);
 338                         histi1[ii] += inp;
 339                         histi2[ii] += 3 * gmx::square(outp) - 1;
 340                         histn[ii]++;
 341                     }
 342                     if ((r2 >= rmin2) && (r2 < rmax2))
 343                     {
 344                         int ii1 = static_cast<int>(invbw * (inp + 1));
 345                         int ii2 = static_cast<int>(invbw * std::abs(outp));
 346
 347                         range_check(ii1, 0, nbin1);
 348                         range_check(ii2, 0, nbin2);
 349                         hist1[ii1]++;
 350                         hist2[ii2]++;
 351                         sum1 += inp;
 352                         sum2 += outp;
 353                         n++;
 354                     }
 355                 }
 356             }
 357         }
 358         ntot += n;
 359         nf++;
 360
 361     } while (read_next_x(oenv, status, &t, x, box));
 362
 363     /* clean up */
 364     sfree(x);
 365     close_trx(status);
 366     gmx_rmpbc_done(gpbc);
 367
 368     /* Add the bin for the exact maximum to the previous bin */
 369     hist1[nbin1 - 1] += hist1[nbin1];
 370     hist2[nbin2 - 1] += hist2[nbin2];
 371
 372     nav     = static_cast<real>(ntot) / (nrefgrp * nf);
 373     normfac = invbw / ntot;
 374
 375     fprintf(stderr, "Average nr of molecules between %g and %g nm: %.1f\n", rmin, rmax, nav);
 376     if (ntot > 0)
 377     {
 378         sum1 /= ntot;
 379         sum2 /= ntot;
 380         fprintf(stderr, "Average cos(theta1)     between %g and %g nm: %6.3f\n", rmin, rmax, sum1);
 381         fprintf(stderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n", rmin, rmax, sum2);
 382     }
 383
 384     sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax);
 385     fp = xvgropen(opt2fn("-o", NFILE, fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv);
 386     if (output_env_get_print_xvgr_codes(oenv))
 387     {
 388         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 389     }
 390     for (i = 0; i < nbin1; i++)
 391     {
 392         fprintf(fp, "%g %g\n", (i + 0.5) * binwidth - 1, 2 * normfac * hist1[i]);
 393     }
 394     xvgrclose(fp);
 395
 396     sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax);
 397     fp = xvgropen(opt2fn("-no", NFILE, fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv);
 398     if (output_env_get_print_xvgr_codes(oenv))
 399     {
 400         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 401     }
 402     for (i = 0; i < nbin2; i++)
 403     {
 404         fprintf(fp, "%g %g\n", (i + 0.5) * binwidth, normfac * hist2[i]);
 405     }
 406     xvgrclose(fp);
 407
 408
 409     sprintf(str, "Solvent orientation");
 410     fp = xvgropen(opt2fn("-ro", NFILE, fnm), str, "r (nm)", "", oenv);
 411     if (output_env_get_print_xvgr_codes(oenv))
 412     {
 413         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 414     }
 415     xvgr_legend(fp, 2, legr, oenv);
 416     for (i = 0; i < nrbin; i++)
 417     {
 418         fprintf(fp,
 419                 "%g %g %g\n",
 420                 (i + 0.5) * rbinw,
 421                 histn[i] ? histi1[i] / histn[i] : 0,
 422                 histn[i] ? histi2[i] / histn[i] : 0);
 423     }
 424     xvgrclose(fp);
 425
 426     sprintf(str, "Cumulative solvent orientation");
 427     fp = xvgropen(opt2fn("-co", NFILE, fnm), str, "r (nm)", "", oenv);
 428     if (output_env_get_print_xvgr_codes(oenv))
 429     {
 430         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 431     }
 432     xvgr_legend(fp, 2, legc, oenv);
 433     normfac = 1.0 / (nrefgrp * nf);
 434     c1      = 0;
 435     c2      = 0;
 436     fprintf(fp, "%g %g %g\n", 0.0, c1, c2);
 437     for (i = 0; i < nrbin; i++)
 438     {
 439         c1 += histi1[i] * normfac;
 440         c2 += histi2[i] * normfac;
 441         fprintf(fp, "%g %g %g\n", (i + 1) * rbinw, c1, c2);
 442     }
 443     xvgrclose(fp);
 444
 445     sprintf(str, "Solvent distribution");
 446     fp = xvgropen(opt2fn("-rc", NFILE, fnm), str, "r (nm)", "molecules/nm", oenv);
 447     if (output_env_get_print_xvgr_codes(oenv))
 448     {
 449         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 450     }
 451     normfac = 1.0 / (rbinw * nf);
 452     for (i = 0; i < nrbin; i++)
 453     {
 454         fprintf(fp, "%g %g\n", (i + 0.5) * rbinw, histn[i] * normfac);
 455     }
 456     xvgrclose(fp);
 457
 458     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 459     do_view(oenv, opt2fn("-no", NFILE, fnm), nullptr);
 460     do_view(oenv, opt2fn("-ro", NFILE, fnm), "-nxy");
 461     do_view(oenv, opt2fn("-co", NFILE, fnm), "-nxy");
 462
 463     return 0;
 464 }