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