src/gromacs/gmxana/gmx_h2order.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 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/commandline/viewit.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxana/princ.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/arraysize.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 /****************************************************************************/
  59 /* This program calculates the ordering of water molecules across a box, as */
  60 /* function of the z-coordinate. This implies averaging over slices and over*/
  61 /* time. Output is the average cos of the angle of the dipole with the      */
  62 /* normal, per slice.                                                       */
  63 /* In addition, it calculates the average dipole moment itself in three     */
  64 /* directions.                                                              */
  65 /****************************************************************************/
  66
  67 static void calc_h2order(const char*             fn,
  68                          const int               index[],
  69                          int                     ngx,
  70                          rvec**                  slDipole,
  71                          real**                  slOrder,
  72                          real*                   slWidth,
  73                          int*                    nslices,
  74                          const t_topology*       top,
  75                          PbcType                 pbcType,
  76                          int                     axis,
  77                          gmx_bool                bMicel,
  78                          int                     micel[],
  79                          int                     nmic,
  80                          const gmx_output_env_t* oenv)
  81 {
  82     rvec *x0,            /* coordinates with pbc */
  83             dipole,      /* dipole moment due to one molecules */
  84             normal, com; /* center of mass of micel, with bMicel */
  85     rvec*        dip;    /* sum of dipoles, unnormalized */
  86     matrix       box;    /* box (3x3) */
  87     t_trxstatus* status;
  88     real         t,                      /* time from trajectory */
  89             *sum,                        /* sum of all cosines of dipoles, per slice */
  90             *frame;                      /* order over one frame */
  91     int natoms,                          /* nr. atoms in trj */
  92             i, j, teller = 0, slice = 0, /* current slice number */
  93             *count;                      /* nr. of atoms in one slice */
  94     gmx_rmpbc_t gpbc = nullptr;
  95
  96     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
  97     {
  98         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
  99     }
 100
 101     if (!*nslices)
 102     {
 103         *nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
 104     }
 105     switch (axis)
 106     {
 107         case 0:
 108             normal[0] = 1;
 109             normal[1] = 0;
 110             normal[2] = 0;
 111             break;
 112         case 1:
 113             normal[0] = 0;
 114             normal[1] = 1;
 115             normal[2] = 0;
 116             break;
 117         case 2:
 118             normal[0] = 0;
 119             normal[1] = 0;
 120             normal[2] = 1;
 121             break;
 122         default: gmx_fatal(FARGS, "No valid value for -axis-. Exiting.\n");
 123     }
 124
 125     clear_rvec(dipole);
 126     snew(count, *nslices);
 127     snew(sum, *nslices);
 128     snew(dip, *nslices);
 129     snew(frame, *nslices);
 130
 131     *slWidth = box[axis][axis] / (*nslices);
 132     fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", *nslices, *slWidth);
 133
 134     teller = 0;
 135
 136     gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
 137     /*********** Start processing trajectory ***********/
 138     do
 139     {
 140         *slWidth = box[axis][axis] / (*nslices);
 141         teller++;
 142
 143         gmx_rmpbc(gpbc, natoms, box, x0);
 144
 145         if (bMicel)
 146         {
 147             calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE);
 148         }
 149
 150         for (i = 0; i < ngx / 3; i++)
 151         {
 152             /* put all waters in box */
 153             for (j = 0; j < DIM; j++)
 154             {
 155                 if (x0[index[3 * i]][j] < 0)
 156                 {
 157                     x0[index[3 * i]][j] += box[j][j];
 158                     x0[index[3 * i + 1]][j] += box[j][j];
 159                     x0[index[3 * i + 2]][j] += box[j][j];
 160                 }
 161                 if (x0[index[3 * i]][j] > box[j][j])
 162                 {
 163                     x0[index[3 * i]][j] -= box[j][j];
 164                     x0[index[3 * i + 1]][j] -= box[j][j];
 165                     x0[index[3 * i + 2]][j] -= box[j][j];
 166                 }
 167             }
 168
 169             for (j = 0; j < DIM; j++)
 170             {
 171                 dipole[j] = x0[index[3 * i]][j] * top->atoms.atom[index[3 * i]].q
 172                             + x0[index[3 * i + 1]][j] * top->atoms.atom[index[3 * i + 1]].q
 173                             + x0[index[3 * i + 2]][j] * top->atoms.atom[index[3 * i + 2]].q;
 174             }
 175
 176             /* now we have a dipole vector. Might as well safe it. Then the
 177                rest depends on whether we're dealing with
 178                a flat or a spherical interface.
 179              */
 180
 181             if (bMicel)
 182             {                                            /* this is for spherical interfaces */
 183                 rvec_sub(com, x0[index[3 * i]], normal); /* vector from Oxygen to COM */
 184                 slice = static_cast<int>(norm(normal) / (*slWidth)); /* spherical slice           */
 185
 186                 sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 187                 frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 188                 count[slice]++;
 189             }
 190             else
 191             {
 192                 /* this is for flat interfaces      */
 193
 194                 /* determine which slice atom is in */
 195                 slice = static_cast<int>(x0[index[3 * i]][axis] / (*slWidth));
 196                 if (slice < 0 || slice >= *nslices)
 197                 {
 198                     fprintf(stderr, "Coordinate: %f ", x0[index[3 * i]][axis]);
 199                     fprintf(stderr, "HELP PANIC! slice = %d, OUT OF RANGE!\n", slice);
 200                 }
 201                 else
 202                 {
 203                     rvec_add(dipole, dip[slice], dip[slice]);
 204                     /* Add dipole to total. mag[slice] is total dipole in axis direction */
 205                     sum[slice] += iprod(dipole, normal) / norm(dipole);
 206                     frame[slice] += iprod(dipole, normal) / norm(dipole);
 207                     /* increase count for that slice */
 208                     count[slice]++;
 209                 }
 210             }
 211         }
 212
 213     } while (read_next_x(oenv, status, &t, x0, box));
 214     /*********** done with status file **********/
 215
 216     fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
 217     gmx_rmpbc_done(gpbc);
 218
 219     for (i = 0; i < *nslices; i++) /* average over frames */
 220     {
 221         fprintf(stderr, "%d waters in slice %d\n", count[i], i);
 222         if (count[i] > 0) /* divide by number of molecules in each slice */
 223         {
 224             sum[i]     = sum[i] / count[i];
 225             dip[i][XX] = dip[i][XX] / count[i];
 226             dip[i][YY] = dip[i][YY] / count[i];
 227             dip[i][ZZ] = dip[i][ZZ] / count[i];
 228         }
 229         else
 230         {
 231             fprintf(stderr, "No water in slice %d\n", i);
 232         }
 233     }
 234
 235     *slOrder  = sum; /* copy a pointer, I hope */
 236     *slDipole = dip;
 237     sfree(x0); /* free memory used by coordinate arrays */
 238 }
 239
 240 static void h2order_plot(rvec dipole[], real order[], const char* afile, int nslices, real slWidth, const gmx_output_env_t* oenv)
 241 {
 242     FILE* ord;      /* xvgr files with order parameters  */
 243     int   slice;    /* loop index     */
 244     char  buf[256]; /* for xvgr title */
 245     real  factor;   /* conversion to Debye from electron*nm */
 246
 247     /*  factor = 1e-9*1.60217733e-19/3.336e-30 */
 248     factor = 1.60217733 / 3.336e-2;
 249     fprintf(stderr, "%d slices\n", nslices);
 250     sprintf(buf, "Water orientation with respect to normal");
 251     ord = xvgropen(afile, buf, "box (nm)", "mu_x, mu_y, mu_z (D), cosine with normal", oenv);
 252
 253     for (slice = 0; slice < nslices; slice++)
 254     {
 255         fprintf(ord, "%8.3f %8.3f %8.3f %8.3f %e\n", slWidth * slice, factor * dipole[slice][XX],
 256                 factor * dipole[slice][YY], factor * dipole[slice][ZZ], order[slice]);
 257     }
 258
 259     xvgrclose(ord);
 260 }
 261
 262 enum
 263 {
 264     axisSEL,
 265     axisZ,
 266     axisY,
 267     axisX,
 268     axisNR
 269 };
 270
 271 int gmx_h2order(int argc, char* argv[])
 272 {
 273     const char* desc[] = {
 274         "[THISMODULE] computes the orientation of water molecules with respect to the normal",
 275         "of the box. The program determines the average cosine of the angle",
 276         "between the dipole moment of water and an axis of the box. The box is",
 277         "divided in slices and the average orientation per slice is printed.",
 278         "Each water molecule is assigned to a slice, per time frame, based on the",
 279         "position of the oxygen. When [TT]-nm[tt] is used, the angle between the water",
 280         "dipole and the axis from the center of mass to the oxygen is calculated",
 281         "instead of the angle between the dipole and a box axis."
 282     };
 283     static const char* axisOption[axisNR + 1] = { nullptr, "Z", "Y", "X", nullptr };
 284     static int         nslices                = 0; /* nr of slices defined       */
 285     // The struct that will hold the parsed user input
 286     t_pargs     pa[]   = { { "-d",
 287                        FALSE,
 288                        etENUM,
 289                        { axisOption },
 290                        "Take the normal on the membrane in direction X, Y or Z." },
 291                      { "-sl",
 292                        FALSE,
 293                        etINT,
 294                        { &nslices },
 295                        "Calculate order parameter as function of boxlength, dividing the box"
 296                        " in this number of slices." } };
 297     const char* bugs[] = {
 298         "The program assigns whole water molecules to a slice, based on the first "
 299         "atom of three in the index file group. It assumes an order O,H,H. "
 300         "Name is not important, but the order is. If this demand is not met, "
 301         "assigning molecules to slices is different."
 302     };
 303
 304     gmx_output_env_t* oenv;
 305     real *            slOrder, /* av. cosine, per slice      */
 306             slWidth = 0.0;     /* width of a slice           */
 307     rvec* slDipole;
 308     char *grpname, /* groupnames                 */
 309             *micname;
 310     int ngx,          /* nr. of atomsin sol group   */
 311             nmic = 0; /* nr. of atoms in micelle    */
 312     t_topology* top;  /* topology           */
 313     PbcType     pbcType;
 314     int *       index, /* indices for solvent group  */
 315             *micelle = nullptr;
 316     gmx_bool bMicel  = FALSE; /* think we're a micel        */
 317     t_filenm fnm[]   = {
 318         /* files for g_order      */
 319         { efTRX, "-f", nullptr, ffREAD },    /* trajectory file            */
 320         { efNDX, nullptr, nullptr, ffREAD }, /* index file         */
 321         { efNDX, "-nm", nullptr, ffOPTRD },  /* index with micelle atoms   */
 322         { efTPR, nullptr, nullptr, ffREAD }, /* topology file              */
 323         { efXVG, "-o", "order", ffWRITE },   /* xvgr output file       */
 324     };
 325
 326 #define NFILE asize(fnm)
 327
 328     // Parse the user input in argv into pa
 329     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
 330                            asize(desc), desc, asize(bugs), bugs, &oenv))
 331     {
 332         return 0;
 333     }
 334
 335     // Process the axis option chosen by the user to set the
 336     // axis used for the computation. The useful choice is an
 337     // axis normal to the membrane. Default is z-axis.
 338     int axis = ZZ;
 339     switch (nenum(axisOption))
 340     {
 341         case axisX: axis = XX; break;
 342         case axisY: axis = YY; break;
 343         case axisZ: axis = ZZ; break;
 344         default: axis = ZZ;
 345     }
 346
 347     bMicel = opt2bSet("-nm", NFILE, fnm);
 348
 349     top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
 350
 351     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &ngx, &index, &grpname);
 352
 353     if (bMicel)
 354     {
 355         rd_index(opt2fn("-nm", NFILE, fnm), 1, &nmic, &micelle, &micname);
 356     }
 357
 358     calc_h2order(ftp2fn(efTRX, NFILE, fnm), index, ngx, &slDipole, &slOrder, &slWidth, &nslices,
 359                  top, pbcType, axis, bMicel, micelle, nmic, oenv);
 360
 361     h2order_plot(slDipole, slOrder, opt2fn("-o", NFILE, fnm), nslices, slWidth, oenv);
 362
 363     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */
 364
 365     return 0;
 366 }