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 "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 "tpxio.h"
  57 #include "gmx_ana.h"
  58
  59 typedef struct {
  60     atom_id i;
  61     real    d2;
  62 } t_order;
  63
  64 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         "[TT]trjorder[tt] 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         "[TT]trjorder[tt] 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 [TT].pdb[tt] 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, title[256];
 138     int             i, j, d, *isize, isize_ref = 0, isize_sol;
 139     atom_id         sa, sr, *swi, **index, *ind_ref = NULL, *ind_sol;
 140     output_env_t    oenv;
 141     t_filenm        fnm[] = {
 142         { efTRX, "-f", NULL, ffREAD  },
 143         { efTPS, NULL, NULL, ffREAD  },
 144         { efNDX, NULL, NULL, ffOPTRD },
 145         { efTRO, "-o", "ordered", ffOPTWR },
 146         { efXVG, "-nshell", "nshell", ffOPTWR }
 147     };
 148 #define NFILE asize(fnm)
 149
 150     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
 151                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
 152
 153     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box, TRUE);
 154     sfree(x);
 155
 156     /* get index groups */
 157     printf("Select %sa group of molecules to be ordered:\n",
 158            bZ ? "" : "a group of reference atoms and ");
 159     snew(grpname, 2);
 160     snew(index, 2);
 161     snew(isize, 2);
 162     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bZ ? 1 : 2,
 163               isize, index, grpname);
 164
 165     if (!bZ)
 166     {
 167         isize_ref = isize[0];
 168         isize_sol = isize[1];
 169         ind_ref   = index[0];
 170         ind_sol   = index[1];
 171     }
 172     else
 173     {
 174         isize_sol = isize[0];
 175         ind_sol   = index[0];
 176     }
 177
 178     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 179     if (natoms > top.atoms.nr)
 180     {
 181         gmx_fatal(FARGS, "Number of atoms in the run input file is larger than in the trjactory");
 182     }
 183     for (i = 0; (i < 2); i++)
 184     {
 185         for (j = 0; (j < isize[i]); j++)
 186         {
 187             if (index[i][j] > natoms)
 188             {
 189                 gmx_fatal(FARGS, "An atom number in group %s is larger than the number of atoms in the trajectory");
 190             }
 191         }
 192     }
 193
 194     if ((isize_sol % na) != 0)
 195     {
 196         gmx_fatal(FARGS, "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 197                   isize[1], na);
 198     }
 199
 200     nwat = isize_sol/na;
 201     if (ref_a > na)
 202     {
 203         gmx_fatal(FARGS, "The reference atom can not be larger than the number of atoms in a molecule");
 204     }
 205     ref_a--;
 206     snew(xsol, nwat);
 207     snew(order, nwat);
 208     snew(swi, natoms);
 209     for (i = 0; (i < natoms); i++)
 210     {
 211         swi[i] = i;
 212     }
 213
 214     out     = NULL;
 215     fp      = NULL;
 216     bNShell = ((opt2bSet("-nshell", NFILE, fnm)) ||
 217                (opt2parg_bSet("-r", asize(pa), pa)));
 218     bPDBout = FALSE;
 219     if (bNShell)
 220     {
 221         rcut2   = rcut*rcut;
 222         fp      = xvgropen(opt2fn("-nshell", NFILE, fnm), "Number of molecules",
 223                            "Time (ps)", "N", oenv);
 224         printf("Will compute the number of molecules within a radius of %g\n",
 225                rcut);
 226     }
 227     if (!bNShell || opt2bSet("-o", NFILE, fnm))
 228     {
 229         bPDBout = (fn2ftp(opt2fn("-o", NFILE, fnm)) == efPDB);
 230         if (bPDBout && !top.atoms.pdbinfo)
 231         {
 232             fprintf(stderr, "Creating pdbfino records\n");
 233             snew(top.atoms.pdbinfo, top.atoms.nr);
 234         }
 235         out = open_trx(opt2fn("-o", NFILE, fnm), "w");
 236     }
 237     gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 238     do
 239     {
 240         gmx_rmpbc(gpbc, natoms, box, x);
 241         set_pbc(&pbc, ePBC, box);
 242
 243         if (ref_a == -1)
 244         {
 245             /* Calculate the COM of all solvent molecules */
 246             for (i = 0; i < nwat; i++)
 247             {
 248                 totmass = 0;
 249                 clear_rvec(xsol[i]);
 250                 for (j = 0; j < na; j++)
 251                 {
 252                     sa       = ind_sol[i*na+j];
 253                     mass     = top.atoms.atom[sa].m;
 254                     totmass += mass;
 255                     for (d = 0; d < DIM; d++)
 256                     {
 257                         xsol[i][d] += mass*x[sa][d];
 258                     }
 259                 }
 260                 svmul(1/totmass, xsol[i], xsol[i]);
 261             }
 262         }
 263         else
 264         {
 265             /* Copy the reference atom of all solvent molecules */
 266             for (i = 0; i < nwat; i++)
 267             {
 268                 copy_rvec(x[ind_sol[i*na+ref_a]], xsol[i]);
 269             }
 270         }
 271
 272         if (bZ)
 273         {
 274             for (i = 0; (i < nwat); i++)
 275             {
 276                 sa           = ind_sol[na*i];
 277                 order[i].i   = sa;
 278                 order[i].d2  = xsol[i][ZZ];
 279             }
 280         }
 281         else if (bCOM)
 282         {
 283             totmass = 0;
 284             clear_rvec(xcom);
 285             for (i = 0; i < isize_ref; i++)
 286             {
 287                 mass     = top.atoms.atom[ind_ref[i]].m;
 288                 totmass += mass;
 289                 for (j = 0; j < DIM; j++)
 290                 {
 291                     xcom[j] += mass*x[ind_ref[i]][j];
 292                 }
 293             }
 294             svmul(1/totmass, xcom, xcom);
 295             for (i = 0; (i < nwat); i++)
 296             {
 297                 sa = ind_sol[na*i];
 298                 pbc_dx(&pbc, xcom, xsol[i], dx);
 299                 order[i].i   = sa;
 300                 order[i].d2  = norm2(dx);
 301             }
 302         }
 303         else
 304         {
 305             /* Set distance to first atom */
 306             for (i = 0; (i < nwat); i++)
 307             {
 308                 sa = ind_sol[na*i];
 309                 pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
 310                 order[i].i   = sa;
 311                 order[i].d2  = norm2(dx);
 312             }
 313             for (j = 1; (j < isize_ref); j++)
 314             {
 315                 sr = ind_ref[j];
 316                 for (i = 0; (i < nwat); i++)
 317                 {
 318                     sa = ind_sol[na*i];
 319                     pbc_dx(&pbc, x[sr], xsol[i], dx);
 320                     n2 = norm2(dx);
 321                     if (n2 < order[i].d2)
 322                     {
 323                         order[i].d2  = n2;
 324                     }
 325                 }
 326             }
 327         }
 328
 329         if (bNShell)
 330         {
 331             ncut = 0;
 332             for (i = 0; (i < nwat); i++)
 333             {
 334                 if (order[i].d2 <= rcut2)
 335                 {
 336                     ncut++;
 337                 }
 338             }
 339             fprintf(fp, "%10.3f  %8d\n", t, ncut);
 340         }
 341         if (out)
 342         {
 343             qsort(order, nwat, sizeof(*order), ocomp);
 344             for (i = 0; (i < nwat); i++)
 345             {
 346                 for (j = 0; (j < na); j++)
 347                 {
 348                     swi[ind_sol[na*i]+j] = order[i].i+j;
 349                 }
 350             }
 351
 352             /* Store the distance as the B-factor */
 353             if (bPDBout)
 354             {
 355                 for (i = 0; (i < nwat); i++)
 356                 {
 357                     for (j = 0; (j < na); j++)
 358                     {
 359                         top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
 360                     }
 361                 }
 362             }
 363             write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, NULL, NULL);
 364         }
 365     }
 366     while (read_next_x(oenv, status, &t, x, box));
 367     close_trj(status);
 368     if (out)
 369     {
 370         close_trx(out);
 371     }
 372     if (fp)
 373     {
 374         ffclose(fp);
 375     }
 376     gmx_rmpbc_done(gpbc);
 377
 378     return 0;
 379 }