src/gromacs/gmxana/gmx_trjorder.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,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020,2021, 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 <cstdlib>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxana/gmx_ana.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/trajectory/trajectoryframe.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 typedef struct
  61 {
  62     int  i;
  63     real d2;
  64 } t_order;
  65
  66 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  67 static t_order* order;
  68
  69 static int ocomp(const void* a, const void* b)
  70 {
  71     const t_order *oa, *ob;
  72
  73     oa = reinterpret_cast<const t_order*>(a);
  74     ob = reinterpret_cast<const t_order*>(b);
  75
  76     if (oa->d2 < ob->d2)
  77     {
  78         return -1;
  79     }
  80     else
  81     {
  82         return 1;
  83     }
  84 }
  85
  86 int gmx_trjorder(int argc, char* argv[])
  87 {
  88     const char* desc[] = {
  89         "[THISMODULE] orders molecules according to the smallest distance",
  90         "to atoms in a reference group",
  91         "or on z-coordinate (with option [TT]-z[tt]).",
  92         "With distance ordering, it will ask for a group of reference",
  93         "atoms and a group of molecules. For each frame of the trajectory",
  94         "the selected molecules will be reordered according to the shortest",
  95         "distance between atom number [TT]-da[tt] in the molecule and all the",
  96         "atoms in the reference group. The center of mass of the molecules can",
  97         "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
  98         "All atoms in the trajectory are written",
  99         "to the output trajectory.[PAR]",
 100         "[THISMODULE] can be useful for e.g. analyzing the n waters closest to a",
 101         "protein.",
 102         "In that case the reference group would be the protein and the group",
 103         "of molecules would consist of all the water atoms. When an index group",
 104         "of the first n waters is made, the ordered trajectory can be used",
 105         "with any GROMACS program to analyze the n closest waters.",
 106         "[PAR]",
 107         "If the output file is a [REF].pdb[ref] file, the distance to the reference target",
 108         "will be stored in the B-factor field in order to color with e.g. Rasmol.",
 109         "[PAR]",
 110         "With option [TT]-nshell[tt] the number of molecules within a shell",
 111         "of radius [TT]-r[tt] around the reference group are printed."
 112     };
 113     static int      na = 3, ref_a = 1;
 114     static real     rcut = 0;
 115     static gmx_bool bCOM = FALSE, bZ = FALSE;
 116     t_pargs         pa[] = {
 117         { "-na", FALSE, etINT, { &na }, "Number of atoms in a molecule" },
 118         { "-da", FALSE, etINT, { &ref_a }, "Atom used for the distance calculation, 0 is COM" },
 119         { "-com",
 120           FALSE,
 121           etBOOL,
 122           { &bCOM },
 123           "Use the distance to the center of mass of the reference group" },
 124         { "-r",
 125           FALSE,
 126           etREAL,
 127           { &rcut },
 128           "Cutoff used for the distance calculation when computing the number of molecules in a "
 129           "shell around e.g. a protein" },
 130         { "-z", FALSE, etBOOL, { &bZ }, "Order molecules on z-coordinate" }
 131     };
 132     FILE*             fp;
 133     t_trxstatus*      out;
 134     t_trxstatus*      status;
 135     gmx_bool          bNShell, bPDBout;
 136     t_topology        top;
 137     PbcType           pbcType;
 138     rvec *            x, *xsol, xcom, dx;
 139     matrix            box;
 140     t_pbc             pbc;
 141     gmx_rmpbc_t       gpbc;
 142     real              t, totmass, mass, rcut2 = 0, n2;
 143     int               natoms, nwat, ncut;
 144     char**            grpname;
 145     int               i, j, d, *isize, isize_ref      = 0, isize_sol;
 146     int               sa, sr, *swi, **index, *ind_ref = nullptr, *ind_sol;
 147     gmx_output_env_t* oenv;
 148     t_filenm          fnm[] = { { efTRX, "-f", nullptr, ffREAD },
 149                        { efTPS, nullptr, nullptr, ffREAD },
 150                        { efNDX, nullptr, nullptr, ffOPTRD },
 151                        { efTRO, "-o", "ordered", ffOPTWR },
 152                        { efXVG, "-nshell", "nshell", ffOPTWR } };
 153 #define NFILE asize(fnm)
 154
 155     if (!parse_common_args(
 156                 &argc, argv, PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 157     {
 158         return 0;
 159     }
 160
 161     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &x, nullptr, box, TRUE);
 162     sfree(x);
 163
 164     /* get index groups */
 165     printf("Select %sa group of molecules to be ordered:\n", bZ ? "" : "a group of reference atoms and ");
 166     snew(grpname, 2);
 167     snew(index, 2);
 168     snew(isize, 2);
 169     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bZ ? 1 : 2, isize, index, grpname);
 170
 171     if (!bZ)
 172     {
 173         isize_ref = isize[0];
 174         isize_sol = isize[1];
 175         ind_ref   = index[0];
 176         ind_sol   = index[1];
 177     }
 178     else
 179     {
 180         isize_sol = isize[0];
 181         ind_sol   = index[0];
 182     }
 183
 184     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 185     if (natoms > top.atoms.nr)
 186     {
 187         gmx_fatal(FARGS, "Number of atoms in the run input file is larger than in the trjactory");
 188     }
 189     for (i = 0; (i < 2); i++)
 190     {
 191         for (j = 0; (j < isize[i]); j++)
 192         {
 193             if (index[i][j] > natoms)
 194             {
 195                 gmx_fatal(FARGS,
 196                           "An atom number in group %s is larger than the number of atoms in the "
 197                           "trajectory",
 198                           grpname[i]);
 199             }
 200         }
 201     }
 202
 203     if ((isize_sol % na) != 0)
 204     {
 205         gmx_fatal(FARGS,
 206                   "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 207                   isize[1],
 208                   na);
 209     }
 210
 211     nwat = isize_sol / na;
 212     if (ref_a > na)
 213     {
 214         gmx_fatal(FARGS,
 215                   "The reference atom can not be larger than the number of atoms in a molecule");
 216     }
 217     ref_a--;
 218     snew(xsol, nwat);
 219     snew(order, nwat);
 220     snew(swi, natoms);
 221     for (i = 0; (i < natoms); i++)
 222     {
 223         swi[i] = i;
 224     }
 225
 226     out     = nullptr;
 227     fp      = nullptr;
 228     bNShell = ((opt2bSet("-nshell", NFILE, fnm)) || (opt2parg_bSet("-r", asize(pa), pa)));
 229     bPDBout = FALSE;
 230     if (bNShell)
 231     {
 232         rcut2 = rcut * rcut;
 233         fp = xvgropen(opt2fn("-nshell", NFILE, fnm), "Number of molecules", "Time (ps)", "N", oenv);
 234         printf("Will compute the number of molecules within a radius of %g\n", rcut);
 235     }
 236     if (!bNShell || opt2bSet("-o", NFILE, fnm))
 237     {
 238         bPDBout = (fn2ftp(opt2fn("-o", NFILE, fnm)) == efPDB);
 239         if (bPDBout && !top.atoms.pdbinfo)
 240         {
 241             fprintf(stderr, "Creating pdbfino records\n");
 242             snew(top.atoms.pdbinfo, top.atoms.nr);
 243         }
 244         out = open_trx(opt2fn("-o", NFILE, fnm), "w");
 245     }
 246     gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
 247     do
 248     {
 249         gmx_rmpbc(gpbc, natoms, box, x);
 250         set_pbc(&pbc, pbcType, box);
 251
 252         if (ref_a == -1)
 253         {
 254             /* Calculate the COM of all solvent molecules */
 255             for (i = 0; i < nwat; i++)
 256             {
 257                 totmass = 0;
 258                 clear_rvec(xsol[i]);
 259                 for (j = 0; j < na; j++)
 260                 {
 261                     sa   = ind_sol[i * na + j];
 262                     mass = top.atoms.atom[sa].m;
 263                     totmass += mass;
 264                     for (d = 0; d < DIM; d++)
 265                     {
 266                         xsol[i][d] += mass * x[sa][d];
 267                     }
 268                 }
 269                 svmul(1.0 / totmass, xsol[i], xsol[i]);
 270             }
 271         }
 272         else
 273         {
 274             /* Copy the reference atom of all solvent molecules */
 275             for (i = 0; i < nwat; i++)
 276             {
 277                 copy_rvec(x[ind_sol[i * na + ref_a]], xsol[i]);
 278             }
 279         }
 280
 281         if (bZ)
 282         {
 283             for (i = 0; (i < nwat); i++)
 284             {
 285                 sa          = ind_sol[na * i];
 286                 order[i].i  = sa;
 287                 order[i].d2 = xsol[i][ZZ];
 288             }
 289         }
 290         else if (bCOM)
 291         {
 292             totmass = 0;
 293             clear_rvec(xcom);
 294             for (i = 0; i < isize_ref; i++)
 295             {
 296                 mass = top.atoms.atom[ind_ref[i]].m;
 297                 totmass += mass;
 298                 for (j = 0; j < DIM; j++)
 299                 {
 300                     xcom[j] += mass * x[ind_ref[i]][j];
 301                 }
 302             }
 303             svmul(1 / totmass, xcom, xcom);
 304             for (i = 0; (i < nwat); i++)
 305             {
 306                 sa = ind_sol[na * i];
 307                 pbc_dx(&pbc, xcom, xsol[i], dx);
 308                 order[i].i  = sa;
 309                 order[i].d2 = norm2(dx);
 310             }
 311         }
 312         else
 313         {
 314             /* Set distance to first atom */
 315             for (i = 0; (i < nwat); i++)
 316             {
 317                 sa = ind_sol[na * i];
 318                 pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
 319                 order[i].i  = sa;
 320                 order[i].d2 = norm2(dx);
 321             }
 322             for (j = 1; (j < isize_ref); j++)
 323             {
 324                 sr = ind_ref[j];
 325                 for (i = 0; (i < nwat); i++)
 326                 {
 327                     pbc_dx(&pbc, x[sr], xsol[i], dx);
 328                     n2 = norm2(dx);
 329                     if (n2 < order[i].d2)
 330                     {
 331                         order[i].d2 = n2;
 332                     }
 333                 }
 334             }
 335         }
 336
 337         if (bNShell)
 338         {
 339             ncut = 0;
 340             for (i = 0; (i < nwat); i++)
 341             {
 342                 if (order[i].d2 <= rcut2)
 343                 {
 344                     ncut++;
 345                 }
 346             }
 347             fprintf(fp, "%10.3f  %8d\n", t, ncut);
 348         }
 349         if (out)
 350         {
 351             qsort(order, nwat, sizeof(*order), ocomp);
 352             for (i = 0; (i < nwat); i++)
 353             {
 354                 for (j = 0; (j < na); j++)
 355                 {
 356                     swi[ind_sol[na * i] + j] = order[i].i + j;
 357                 }
 358             }
 359
 360             /* Store the distance as the B-factor */
 361             if (bPDBout)
 362             {
 363                 for (i = 0; (i < nwat); i++)
 364                 {
 365                     for (j = 0; (j < na); j++)
 366                     {
 367                         top.atoms.pdbinfo[order[i].i + j].bfac = std::sqrt(order[i].d2);
 368                     }
 369                 }
 370             }
 371             write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, nullptr, nullptr);
 372         }
 373     } while (read_next_x(oenv, status, &t, x, box));
 374     close_trx(status);
 375     if (out)
 376     {
 377         close_trx(out);
 378     }
 379     if (fp)
 380     {
 381         xvgrclose(fp);
 382     }
 383     gmx_rmpbc_done(gpbc);
 384
 385     return 0;
 386 }