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