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