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