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