src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44 #include "statutil.h"
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "pbc.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "index.h"
  55 #include "mshift.h"
  56 #include "xvgr.h"
  57 #include "princ.h"
  58 #include "rmpbc.h"
  59 #include "txtdump.h"
  60 #include "tpxio.h"
  61 #include "gmx_ana.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     return -1;
  79   else
  80     return 1;
  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   CopyRight(stderr,argv[0]);
 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     isize_ref = isize[0];
 168     isize_sol = isize[1];
 169     ind_ref   = index[0];
 170     ind_sol   = index[1];
 171   } else {
 172     isize_sol = isize[0];
 173     ind_sol   = index[0];
 174   }
 175
 176   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 177   if (natoms > top.atoms.nr)
 178     gmx_fatal(FARGS,"Number of atoms in the run input file is larger than in the trjactory");
 179   for(i=0; (i<2); i++)
 180     for(j=0; (j<isize[i]); j++)
 181       if (index[i][j] > natoms)
 182         gmx_fatal(FARGS,"An atom number in group %s is larger than the number of atoms in the trajectory");
 183
 184   if ((isize_sol % na) != 0)
 185     gmx_fatal(FARGS,"Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 186                 isize[1],na);
 187
 188   nwat = isize_sol/na;
 189   if (ref_a > na)
 190     gmx_fatal(FARGS,"The reference atom can not be larger than the number of atoms in a molecule");
 191   ref_a--;
 192   snew(xsol,nwat);
 193   snew(order,nwat);
 194   snew(swi,natoms);
 195   for(i=0; (i<natoms); i++)
 196     swi[i] = i;
 197
 198   out     = NULL;
 199   fp      = NULL;
 200   bNShell = ((opt2bSet("-nshell",NFILE,fnm)) ||
 201              (opt2parg_bSet("-r",asize(pa),pa)));
 202   bPDBout = FALSE;
 203   if (bNShell) {
 204     rcut2   = rcut*rcut;
 205     fp = xvgropen(opt2fn("-nshell",NFILE,fnm),"Number of molecules",
 206                   "Time (ps)","N",oenv);
 207     printf("Will compute the number of molecules within a radius of %g\n",
 208            rcut);
 209   }
 210   if (!bNShell || opt2bSet("-o",NFILE,fnm)) {
 211     bPDBout = (fn2ftp(opt2fn("-o",NFILE,fnm)) == efPDB);
 212     if (bPDBout && !top.atoms.pdbinfo) {
 213       fprintf(stderr,"Creating pdbfino records\n");
 214       snew(top.atoms.pdbinfo,top.atoms.nr);
 215     }
 216     out = open_trx(opt2fn("-o",NFILE,fnm),"w");
 217   }
 218   gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 219   do {
 220     gmx_rmpbc(gpbc,natoms,box,x);
 221     set_pbc(&pbc,ePBC,box);
 222
 223     if (ref_a == -1) {
 224       /* Calculate the COM of all solvent molecules */
 225       for(i=0; i<nwat; i++) {
 226         totmass = 0;
 227         clear_rvec(xsol[i]);
 228         for(j=0; j<na; j++) {
 229           sa = ind_sol[i*na+j];
 230           mass = top.atoms.atom[sa].m;
 231           totmass += mass;
 232           for(d=0; d<DIM; d++) {
 233             xsol[i][d] += mass*x[sa][d];
 234           }
 235         }
 236         svmul(1/totmass,xsol[i],xsol[i]);
 237       }
 238     } else {
 239       /* Copy the reference atom of all solvent molecules */
 240       for(i=0; i<nwat; i++) {
 241         copy_rvec(x[ind_sol[i*na+ref_a]],xsol[i]);
 242       }
 243     }
 244
 245     if (bZ) {
 246       for(i=0; (i<nwat); i++) {
 247         sa = ind_sol[na*i];
 248         order[i].i   = sa;
 249         order[i].d2  = xsol[i][ZZ];
 250       }
 251     } else if (bCOM) {
 252       totmass = 0;
 253       clear_rvec(xcom);
 254       for(i=0; i<isize_ref; i++) {
 255         mass = top.atoms.atom[ind_ref[i]].m;
 256         totmass += mass;
 257         for(j=0; j<DIM; j++)
 258           xcom[j] += mass*x[ind_ref[i]][j];
 259       }
 260       svmul(1/totmass,xcom,xcom);
 261       for(i=0; (i<nwat); i++) {
 262         sa = ind_sol[na*i];
 263         pbc_dx(&pbc,xcom,xsol[i],dx);
 264         order[i].i   = sa;
 265         order[i].d2  = norm2(dx);
 266       }
 267     } else {
 268       /* Set distance to first atom */
 269       for(i=0; (i<nwat); i++) {
 270         sa = ind_sol[na*i];
 271         pbc_dx(&pbc,x[ind_ref[0]],xsol[i],dx);
 272         order[i].i   = sa;
 273         order[i].d2  = norm2(dx);
 274       }
 275       for(j=1; (j<isize_ref); j++) {
 276         sr = ind_ref[j];
 277         for(i=0; (i<nwat); i++) {
 278           sa = ind_sol[na*i];
 279           pbc_dx(&pbc,x[sr],xsol[i],dx);
 280           n2 = norm2(dx);
 281           if (n2 < order[i].d2)
 282             order[i].d2  = n2;
 283         }
 284       }
 285     }
 286
 287     if (bNShell) {
 288       ncut = 0;
 289       for(i=0; (i<nwat); i++)
 290         if (order[i].d2 <= rcut2)
 291           ncut++;
 292       fprintf(fp,"%10.3f  %8d\n",t,ncut);
 293     }
 294     if (out) {
 295       qsort(order,nwat,sizeof(*order),ocomp);
 296       for(i=0; (i<nwat); i++)
 297         for(j=0; (j<na); j++)
 298           swi[ind_sol[na*i]+j] = order[i].i+j;
 299
 300       /* Store the distance as the B-factor */
 301       if (bPDBout) {
 302         for(i=0; (i<nwat); i++) {
 303           for(j=0; (j<na); j++) {
 304             top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
 305           }
 306         }
 307       }
 308       write_trx(out,natoms,swi,&top.atoms,0,t,box,x,NULL,NULL);
 309     }
 310   } while(read_next_x(oenv,status,&t,natoms,x,box));
 311   close_trj(status);
 312   if (out)
 313     close_trx(out);
 314   if (fp)
 315     ffclose(fp);
 316   gmx_rmpbc_done(gpbc);
 317
 318   thanx(stderr);
 319
 320   return 0;
 321 }