src/gromacs/gmxana/gmx_trjorder.c

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