src/tools/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     return -1;
  76   else
  77     return 1;
  78 }
  79
  80 int gmx_trjorder(int argc,char *argv[])
  81 {
  82   const char *desc[] = {
  83     "[TT]trjorder[tt] orders molecules according to the smallest distance",
  84     "to atoms in a reference group",
  85     "or on z-coordinate (with option [TT]-z[tt]).",
  86     "With distance ordering, it will ask for a group of reference",
  87     "atoms and a group of molecules. For each frame of the trajectory",
  88     "the selected molecules will be reordered according to the shortest",
  89     "distance between atom number [TT]-da[tt] in the molecule and all the",
  90     "atoms in the reference group. The center of mass of the molecules can",
  91     "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
  92     "All atoms in the trajectory are written",
  93     "to the output trajectory.[PAR]",
  94     "[TT]trjorder[tt] can be useful for e.g. analyzing the n waters closest to a",
  95     "protein.",
  96     "In that case the reference group would be the protein and the group",
  97     "of molecules would consist of all the water atoms. When an index group",
  98     "of the first n waters is made, the ordered trajectory can be used",
  99     "with any Gromacs program to analyze the n closest waters.",
 100     "[PAR]",
 101     "If the output file is a [TT].pdb[tt] file, the distance to the reference target",
 102     "will be stored in the B-factor field in order to color with e.g. rasmol.",
 103     "[PAR]",
 104     "With option [TT]-nshell[tt] the number of molecules within a shell",
 105     "of radius [TT]-r[tt] around the reference group are printed."
 106   };
 107   static int na=3,ref_a=1;
 108   static real rcut=0;
 109   static gmx_bool bCOM=FALSE,bZ=FALSE;
 110   t_pargs pa[] = {
 111     { "-na", FALSE, etINT,  {&na},
 112       "Number of atoms in a molecule" },
 113     { "-da", FALSE, etINT,  {&ref_a},
 114       "Atom used for the distance calculation, 0 is COM" },
 115     { "-com", FALSE, etBOOL, {&bCOM},
 116       "Use the distance to the center of mass of the reference group" },
 117     { "-r",  FALSE, etREAL, {&rcut},
 118       "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
 119     { "-z", FALSE, etBOOL, {&bZ},
 120       "Order molecules on z-coordinate" }
 121   };
 122   FILE       *fp;
 123   t_trxstatus *out;
 124   t_trxstatus *status;
 125   gmx_bool       bNShell,bPDBout;
 126   t_topology top;
 127   int        ePBC;
 128   rvec       *x,*xsol,xcom,dx;
 129   matrix     box;
 130   t_pbc      pbc;
 131   gmx_rmpbc_t gpbc;
 132   real       t,totmass,mass,rcut2=0,n2;
 133   int        natoms,nwat,ncut;
 134   char       **grpname,title[256];
 135   int        i,j,d,*isize,isize_ref=0,isize_sol;
 136   atom_id    sa,sr,*swi,**index,*ind_ref=NULL,*ind_sol;
 137   output_env_t oenv;
 138   t_filenm fnm[] = {
 139     { efTRX, "-f", NULL, ffREAD  },
 140     { efTPS, NULL, NULL, ffREAD  },
 141     { efNDX, NULL, NULL, ffOPTRD },
 142     { efTRO, "-o", "ordered", ffOPTWR },
 143     { efXVG, "-nshell", "nshell", ffOPTWR }
 144   };
 145 #define NFILE asize(fnm)
 146
 147   CopyRight(stderr,argv[0]);
 148   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 149                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 150
 151   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
 152   sfree(x);
 153
 154   /* get index groups */
 155   printf("Select %sa group of molecules to be ordered:\n",
 156          bZ ? "" : "a group of reference atoms and ");
 157   snew(grpname,2);
 158   snew(index,2);
 159   snew(isize,2);
 160   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),bZ ? 1 : 2,
 161             isize,index,grpname);
 162
 163   if (!bZ) {
 164     isize_ref = isize[0];
 165     isize_sol = isize[1];
 166     ind_ref   = index[0];
 167     ind_sol   = index[1];
 168   } else {
 169     isize_sol = isize[0];
 170     ind_sol   = index[0];
 171   }
 172
 173   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 174   if (natoms > top.atoms.nr)
 175     gmx_fatal(FARGS,"Number of atoms in the run input file is larger than in the trjactory");
 176   for(i=0; (i<2); i++)
 177     for(j=0; (j<isize[i]); j++)
 178       if (index[i][j] > natoms)
 179         gmx_fatal(FARGS,"An atom number in group %s is larger than the number of atoms in the trajectory");
 180
 181   if ((isize_sol % na) != 0)
 182     gmx_fatal(FARGS,"Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 183                 isize[1],na);
 184
 185   nwat = isize_sol/na;
 186   if (ref_a > na)
 187     gmx_fatal(FARGS,"The reference atom can not be larger than the number of atoms in a molecule");
 188   ref_a--;
 189   snew(xsol,nwat);
 190   snew(order,nwat);
 191   snew(swi,natoms);
 192   for(i=0; (i<natoms); i++)
 193     swi[i] = i;
 194
 195   out     = NULL;
 196   fp      = NULL;
 197   bNShell = ((opt2bSet("-nshell",NFILE,fnm)) ||
 198              (opt2parg_bSet("-r",asize(pa),pa)));
 199   bPDBout = FALSE;
 200   if (bNShell) {
 201     rcut2   = rcut*rcut;
 202     fp = xvgropen(opt2fn("-nshell",NFILE,fnm),"Number of molecules",
 203                   "Time (ps)","N",oenv);
 204     printf("Will compute the number of molecules within a radius of %g\n",
 205            rcut);
 206   }
 207   if (!bNShell || opt2bSet("-o",NFILE,fnm)) {
 208     bPDBout = (fn2ftp(opt2fn("-o",NFILE,fnm)) == efPDB);
 209     if (bPDBout && !top.atoms.pdbinfo) {
 210       fprintf(stderr,"Creating pdbfino records\n");
 211       snew(top.atoms.pdbinfo,top.atoms.nr);
 212     }
 213     out = open_trx(opt2fn("-o",NFILE,fnm),"w");
 214   }
 215   gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 216   do {
 217     gmx_rmpbc(gpbc,natoms,box,x);
 218     set_pbc(&pbc,ePBC,box);
 219
 220     if (ref_a == -1) {
 221       /* Calculate the COM of all solvent molecules */
 222       for(i=0; i<nwat; i++) {
 223         totmass = 0;
 224         clear_rvec(xsol[i]);
 225         for(j=0; j<na; j++) {
 226           sa = ind_sol[i*na+j];
 227           mass = top.atoms.atom[sa].m;
 228           totmass += mass;
 229           for(d=0; d<DIM; d++) {
 230             xsol[i][d] += mass*x[sa][d];
 231           }
 232         }
 233         svmul(1/totmass,xsol[i],xsol[i]);
 234       }
 235     } else {
 236       /* Copy the reference atom of all solvent molecules */
 237       for(i=0; i<nwat; i++) {
 238         copy_rvec(x[ind_sol[i*na+ref_a]],xsol[i]);
 239       }
 240     }
 241
 242     if (bZ) {
 243       for(i=0; (i<nwat); i++) {
 244         sa = ind_sol[na*i];
 245         order[i].i   = sa;
 246         order[i].d2  = xsol[i][ZZ];
 247       }
 248     } else if (bCOM) {
 249       totmass = 0;
 250       clear_rvec(xcom);
 251       for(i=0; i<isize_ref; i++) {
 252         mass = top.atoms.atom[ind_ref[i]].m;
 253         totmass += mass;
 254         for(j=0; j<DIM; j++)
 255           xcom[j] += mass*x[ind_ref[i]][j];
 256       }
 257       svmul(1/totmass,xcom,xcom);
 258       for(i=0; (i<nwat); i++) {
 259         sa = ind_sol[na*i];
 260         pbc_dx(&pbc,xcom,xsol[i],dx);
 261         order[i].i   = sa;
 262         order[i].d2  = norm2(dx);
 263       }
 264     } else {
 265       /* Set distance to first atom */
 266       for(i=0; (i<nwat); i++) {
 267         sa = ind_sol[na*i];
 268         pbc_dx(&pbc,x[ind_ref[0]],xsol[i],dx);
 269         order[i].i   = sa;
 270         order[i].d2  = norm2(dx);
 271       }
 272       for(j=1; (j<isize_ref); j++) {
 273         sr = ind_ref[j];
 274         for(i=0; (i<nwat); i++) {
 275           sa = ind_sol[na*i];
 276           pbc_dx(&pbc,x[sr],xsol[i],dx);
 277           n2 = norm2(dx);
 278           if (n2 < order[i].d2)
 279             order[i].d2  = n2;
 280         }
 281       }
 282     }
 283
 284     if (bNShell) {
 285       ncut = 0;
 286       for(i=0; (i<nwat); i++)
 287         if (order[i].d2 <= rcut2)
 288           ncut++;
 289       fprintf(fp,"%10.3f  %8d\n",t,ncut);
 290     }
 291     if (out) {
 292       qsort(order,nwat,sizeof(*order),ocomp);
 293       for(i=0; (i<nwat); i++)
 294         for(j=0; (j<na); j++)
 295           swi[ind_sol[na*i]+j] = order[i].i+j;
 296
 297       /* Store the distance as the B-factor */
 298       if (bPDBout) {
 299         for(i=0; (i<nwat); i++) {
 300           for(j=0; (j<na); j++) {
 301             top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
 302           }
 303         }
 304       }
 305       write_trx(out,natoms,swi,&top.atoms,0,t,box,x,NULL,NULL);
 306     }
 307   } while(read_next_x(oenv,status,&t,natoms,x,box));
 308   close_trj(status);
 309   if (out)
 310     close_trx(out);
 311   if (fp)
 312     ffclose(fp);
 313   gmx_rmpbc_done(gpbc);
 314
 315   thanx(stderr);
 316
 317   return 0;
 318 }