src/gromacs/gmxana/gmx_trjorder.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41 #include "statutil.h"
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "pbc.h"
  48 #include "gromacs/fileio/futil.h"
  49 #include "statutil.h"
  50 #include "index.h"
  51 #include "mshift.h"
  52 #include "xvgr.h"
  53 #include "princ.h"
  54 #include "rmpbc.h"
  55 #include "txtdump.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "gromacs/fileio/trxio.h"
  58 #include "gmx_ana.h"
  59
  60 typedef struct {
  61     atom_id i;
  62     real    d2;
  63 } t_order;
  64
  65 t_order *order;
  66
  67 static int ocomp(const void *a, const void *b)
  68 {
  69     t_order *oa, *ob;
  70
  71     oa = (t_order *)a;
  72     ob = (t_order *)b;
  73
  74     if (oa->d2 < ob->d2)
  75     {
  76         return -1;
  77     }
  78     else
  79     {
  80         return 1;
  81     }
  82 }
  83
  84 int gmx_trjorder(int argc, char *argv[])
  85 {
  86     const char     *desc[] = {
  87         "[TT]trjorder[tt] orders molecules according to the smallest distance",
  88         "to atoms in a reference group",
  89         "or on z-coordinate (with option [TT]-z[tt]).",
  90         "With distance ordering, it will ask for a group of reference",
  91         "atoms and a group of molecules. For each frame of the trajectory",
  92         "the selected molecules will be reordered according to the shortest",
  93         "distance between atom number [TT]-da[tt] in the molecule and all the",
  94         "atoms in the reference group. The center of mass of the molecules can",
  95         "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
  96         "All atoms in the trajectory are written",
  97         "to the output trajectory.[PAR]",
  98         "[TT]trjorder[tt] can be useful for e.g. analyzing the n waters closest to a",
  99         "protein.",
 100         "In that case the reference group would be the protein and the group",
 101         "of molecules would consist of all the water atoms. When an index group",
 102         "of the first n waters is made, the ordered trajectory can be used",
 103         "with any Gromacs program to analyze the n closest waters.",
 104         "[PAR]",
 105         "If the output file is a [TT].pdb[tt] file, the distance to the reference target",
 106         "will be stored in the B-factor field in order to color with e.g. Rasmol.",
 107         "[PAR]",
 108         "With option [TT]-nshell[tt] the number of molecules within a shell",
 109         "of radius [TT]-r[tt] around the reference group are printed."
 110     };
 111     static int      na   = 3, ref_a = 1;
 112     static real     rcut = 0;
 113     static gmx_bool bCOM = FALSE, bZ = FALSE;
 114     t_pargs         pa[] = {
 115         { "-na", FALSE, etINT,  {&na},
 116           "Number of atoms in a molecule" },
 117         { "-da", FALSE, etINT,  {&ref_a},
 118           "Atom used for the distance calculation, 0 is COM" },
 119         { "-com", FALSE, etBOOL, {&bCOM},
 120           "Use the distance to the center of mass of the reference group" },
 121         { "-r",  FALSE, etREAL, {&rcut},
 122           "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
 123         { "-z", FALSE, etBOOL, {&bZ},
 124           "Order molecules on z-coordinate" }
 125     };
 126     FILE           *fp;
 127     t_trxstatus    *out;
 128     t_trxstatus    *status;
 129     gmx_bool        bNShell, bPDBout;
 130     t_topology      top;
 131     int             ePBC;
 132     rvec           *x, *xsol, xcom, dx;
 133     matrix          box;
 134     t_pbc           pbc;
 135     gmx_rmpbc_t     gpbc;
 136     real            t, totmass, mass, rcut2 = 0, n2;
 137     int             natoms, nwat, ncut;
 138     char          **grpname, title[256];
 139     int             i, j, d, *isize, isize_ref = 0, isize_sol;
 140     atom_id         sa, sr, *swi, **index, *ind_ref = NULL, *ind_sol;
 141     output_env_t    oenv;
 142     t_filenm        fnm[] = {
 143         { efTRX, "-f", NULL, ffREAD  },
 144         { efTPS, NULL, NULL, ffREAD  },
 145         { efNDX, NULL, NULL, ffOPTRD },
 146         { efTRO, "-o", "ordered", ffOPTWR },
 147         { efXVG, "-nshell", "nshell", ffOPTWR }
 148     };
 149 #define NFILE asize(fnm)
 150
 151     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
 152                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 153     {
 154         return 0;
 155     }
 156
 157     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box, TRUE);
 158     sfree(x);
 159
 160     /* get index groups */
 161     printf("Select %sa group of molecules to be ordered:\n",
 162            bZ ? "" : "a group of reference atoms and ");
 163     snew(grpname, 2);
 164     snew(index, 2);
 165     snew(isize, 2);
 166     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bZ ? 1 : 2,
 167               isize, index, grpname);
 168
 169     if (!bZ)
 170     {
 171         isize_ref = isize[0];
 172         isize_sol = isize[1];
 173         ind_ref   = index[0];
 174         ind_sol   = index[1];
 175     }
 176     else
 177     {
 178         isize_sol = isize[0];
 179         ind_sol   = index[0];
 180     }
 181
 182     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 183     if (natoms > top.atoms.nr)
 184     {
 185         gmx_fatal(FARGS, "Number of atoms in the run input file is larger than in the trjactory");
 186     }
 187     for (i = 0; (i < 2); i++)
 188     {
 189         for (j = 0; (j < isize[i]); j++)
 190         {
 191             if (index[i][j] > natoms)
 192             {
 193                 gmx_fatal(FARGS, "An atom number in group %s is larger than the number of atoms in the trajectory");
 194             }
 195         }
 196     }
 197
 198     if ((isize_sol % na) != 0)
 199     {
 200         gmx_fatal(FARGS, "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 201                   isize[1], na);
 202     }
 203
 204     nwat = isize_sol/na;
 205     if (ref_a > na)
 206     {
 207         gmx_fatal(FARGS, "The reference atom can not be larger than the number of atoms in a molecule");
 208     }
 209     ref_a--;
 210     snew(xsol, nwat);
 211     snew(order, nwat);
 212     snew(swi, natoms);
 213     for (i = 0; (i < natoms); i++)
 214     {
 215         swi[i] = i;
 216     }
 217
 218     out     = NULL;
 219     fp      = NULL;
 220     bNShell = ((opt2bSet("-nshell", NFILE, fnm)) ||
 221                (opt2parg_bSet("-r", asize(pa), pa)));
 222     bPDBout = FALSE;
 223     if (bNShell)
 224     {
 225         rcut2   = rcut*rcut;
 226         fp      = xvgropen(opt2fn("-nshell", NFILE, fnm), "Number of molecules",
 227                            "Time (ps)", "N", oenv);
 228         printf("Will compute the number of molecules within a radius of %g\n",
 229                rcut);
 230     }
 231     if (!bNShell || opt2bSet("-o", NFILE, fnm))
 232     {
 233         bPDBout = (fn2ftp(opt2fn("-o", NFILE, fnm)) == efPDB);
 234         if (bPDBout && !top.atoms.pdbinfo)
 235         {
 236             fprintf(stderr, "Creating pdbfino records\n");
 237             snew(top.atoms.pdbinfo, top.atoms.nr);
 238         }
 239         out = open_trx(opt2fn("-o", NFILE, fnm), "w");
 240     }
 241     gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 242     do
 243     {
 244         gmx_rmpbc(gpbc, natoms, box, x);
 245         set_pbc(&pbc, ePBC, box);
 246
 247         if (ref_a == -1)
 248         {
 249             /* Calculate the COM of all solvent molecules */
 250             for (i = 0; i < nwat; i++)
 251             {
 252                 totmass = 0;
 253                 clear_rvec(xsol[i]);
 254                 for (j = 0; j < na; j++)
 255                 {
 256                     sa       = ind_sol[i*na+j];
 257                     mass     = top.atoms.atom[sa].m;
 258                     totmass += mass;
 259                     for (d = 0; d < DIM; d++)
 260                     {
 261                         xsol[i][d] += mass*x[sa][d];
 262                     }
 263                 }
 264                 svmul(1/totmass, xsol[i], xsol[i]);
 265             }
 266         }
 267         else
 268         {
 269             /* Copy the reference atom of all solvent molecules */
 270             for (i = 0; i < nwat; i++)
 271             {
 272                 copy_rvec(x[ind_sol[i*na+ref_a]], xsol[i]);
 273             }
 274         }
 275
 276         if (bZ)
 277         {
 278             for (i = 0; (i < nwat); i++)
 279             {
 280                 sa           = ind_sol[na*i];
 281                 order[i].i   = sa;
 282                 order[i].d2  = xsol[i][ZZ];
 283             }
 284         }
 285         else if (bCOM)
 286         {
 287             totmass = 0;
 288             clear_rvec(xcom);
 289             for (i = 0; i < isize_ref; i++)
 290             {
 291                 mass     = top.atoms.atom[ind_ref[i]].m;
 292                 totmass += mass;
 293                 for (j = 0; j < DIM; j++)
 294                 {
 295                     xcom[j] += mass*x[ind_ref[i]][j];
 296                 }
 297             }
 298             svmul(1/totmass, xcom, xcom);
 299             for (i = 0; (i < nwat); i++)
 300             {
 301                 sa = ind_sol[na*i];
 302                 pbc_dx(&pbc, xcom, xsol[i], dx);
 303                 order[i].i   = sa;
 304                 order[i].d2  = norm2(dx);
 305             }
 306         }
 307         else
 308         {
 309             /* Set distance to first atom */
 310             for (i = 0; (i < nwat); i++)
 311             {
 312                 sa = ind_sol[na*i];
 313                 pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
 314                 order[i].i   = sa;
 315                 order[i].d2  = norm2(dx);
 316             }
 317             for (j = 1; (j < isize_ref); j++)
 318             {
 319                 sr = ind_ref[j];
 320                 for (i = 0; (i < nwat); i++)
 321                 {
 322                     sa = ind_sol[na*i];
 323                     pbc_dx(&pbc, x[sr], xsol[i], dx);
 324                     n2 = norm2(dx);
 325                     if (n2 < order[i].d2)
 326                     {
 327                         order[i].d2  = n2;
 328                     }
 329                 }
 330             }
 331         }
 332
 333         if (bNShell)
 334         {
 335             ncut = 0;
 336             for (i = 0; (i < nwat); i++)
 337             {
 338                 if (order[i].d2 <= rcut2)
 339                 {
 340                     ncut++;
 341                 }
 342             }
 343             fprintf(fp, "%10.3f  %8d\n", t, ncut);
 344         }
 345         if (out)
 346         {
 347             qsort(order, nwat, sizeof(*order), ocomp);
 348             for (i = 0; (i < nwat); i++)
 349             {
 350                 for (j = 0; (j < na); j++)
 351                 {
 352                     swi[ind_sol[na*i]+j] = order[i].i+j;
 353                 }
 354             }
 355
 356             /* Store the distance as the B-factor */
 357             if (bPDBout)
 358             {
 359                 for (i = 0; (i < nwat); i++)
 360                 {
 361                     for (j = 0; (j < na); j++)
 362                     {
 363                         top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
 364                     }
 365                 }
 366             }
 367             write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, NULL, NULL);
 368         }
 369     }
 370     while (read_next_x(oenv, status, &t, x, box));
 371     close_trj(status);
 372     if (out)
 373     {
 374         close_trx(out);
 375     }
 376     if (fp)
 377     {
 378         ffclose(fp);
 379     }
 380     gmx_rmpbc_done(gpbc);
 381
 382     return 0;
 383 }