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",
  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     "trjorder 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 pdb 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 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   int        status,out;
 124   bool       bNShell,bPDBout;
 125   t_topology top;
 126   int        ePBC;
 127   rvec       *x,*xsol,xcom,dx;
 128   matrix     box;
 129   t_pbc      pbc;
 130   real       t,totmass,mass,rcut2=0,n2;
 131   int        natoms,nwat,ncut;
 132   char       **grpname,title[256];
 133   int        i,j,d,*isize,isize_ref=0,isize_sol;
 134   atom_id    sa,sr,*swi,**index,*ind_ref=NULL,*ind_sol;
 135   output_env_t oenv;
 136   t_filenm fnm[] = {
 137     { efTRX, "-f", NULL, ffREAD  },
 138     { efTPS, NULL, NULL, ffREAD  },
 139     { efNDX, NULL, NULL, ffOPTRD },
 140     { efTRO, "-o", "ordered", ffOPTWR },
 141     { efXVG, "-nshell", "nshell", ffOPTWR }
 142   };
 143 #define NFILE asize(fnm)
 144
 145   CopyRight(stderr,argv[0]);
 146   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 147                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 148
 149   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
 150   sfree(x);
 151
 152   /* get index groups */
 153   printf("Select %sa group of molecules to be ordered:\n",
 154          bZ ? "" : "a group of reference atoms and ");
 155   snew(grpname,2);
 156   snew(index,2);
 157   snew(isize,2);
 158   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),bZ ? 1 : 2,
 159             isize,index,grpname);
 160
 161   if (!bZ) {
 162     isize_ref = isize[0];
 163     isize_sol = isize[1];
 164     ind_ref   = index[0];
 165     ind_sol   = index[1];
 166   } else {
 167     isize_sol = isize[0];
 168     ind_sol   = index[0];
 169   }
 170
 171   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 172   if (natoms > top.atoms.nr)
 173     gmx_fatal(FARGS,"Number of atoms in the run input file is larger than in the trjactory");
 174   for(i=0; (i<2); i++)
 175     for(j=0; (j<isize[i]); j++)
 176       if (index[i][j] > natoms)
 177         gmx_fatal(FARGS,"An atom number in group %s is larger than the number of atoms in the trajectory");
 178
 179   if ((isize_sol % na) != 0)
 180     gmx_fatal(FARGS,"Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
 181                 isize[1],na);
 182
 183   nwat = isize_sol/na;
 184   if (ref_a > na)
 185     gmx_fatal(FARGS,"The reference atom can not be larger than the number of atoms in a molecule");
 186   ref_a--;
 187   snew(xsol,nwat);
 188   snew(order,nwat);
 189   snew(swi,natoms);
 190   for(i=0; (i<natoms); i++)
 191     swi[i] = i;
 192
 193   out     = -1;
 194   fp      = NULL;
 195   bNShell = ((opt2bSet("-nshell",NFILE,fnm)) ||
 196              (opt2parg_bSet("-r",asize(pa),pa)));
 197   bPDBout = FALSE;
 198   if (bNShell) {
 199     rcut2   = rcut*rcut;
 200     fp = xvgropen(opt2fn("-nshell",NFILE,fnm),"Number of molecules",
 201                   "Time (ps)","N",oenv);
 202     printf("Will compute the number of molecules within a radius of %g\n",
 203            rcut);
 204   }
 205   if (!bNShell || opt2bSet("-o",NFILE,fnm)) {
 206     bPDBout = (fn2ftp(opt2fn("-o",NFILE,fnm)) == efPDB);
 207     if (bPDBout && !top.atoms.pdbinfo) {
 208       fprintf(stderr,"Creating pdbfino records\n");
 209       snew(top.atoms.pdbinfo,top.atoms.nr);
 210     }
 211     out = open_trx(opt2fn("-o",NFILE,fnm),"w");
 212   }
 213   do {
 214     rm_pbc(&top.idef,ePBC,natoms,box,x,x);
 215     set_pbc(&pbc,ePBC,box);
 216
 217     if (ref_a == -1) {
 218       /* Calculate the COM of all solvent molecules */
 219       for(i=0; i<nwat; i++) {
 220         totmass = 0;
 221         clear_rvec(xsol[i]);
 222         for(j=0; j<na; j++) {
 223           sa = ind_sol[i*na+j];
 224           mass = top.atoms.atom[sa].m;
 225           totmass += mass;
 226           for(d=0; d<DIM; d++) {
 227             xsol[i][d] += mass*x[sa][d];
 228           }
 229         }
 230         svmul(1/totmass,xsol[i],xsol[i]);
 231       }
 232     } else {
 233       /* Copy the reference atom of all solvent molecules */
 234       for(i=0; i<nwat; i++) {
 235         copy_rvec(x[ind_sol[i*na+ref_a]],xsol[i]);
 236       }
 237     }
 238
 239     if (bZ) {
 240       for(i=0; (i<nwat); i++) {
 241         sa = ind_sol[na*i];
 242         order[i].i   = sa;
 243         order[i].d2  = xsol[i][ZZ];
 244       }
 245     } else if (bCOM) {
 246       totmass = 0;
 247       clear_rvec(xcom);
 248       for(i=0; i<isize_ref; i++) {
 249         mass = top.atoms.atom[ind_ref[i]].m;
 250         totmass += mass;
 251         for(j=0; j<DIM; j++)
 252           xcom[j] += mass*x[ind_ref[i]][j];
 253       }
 254       svmul(1/totmass,xcom,xcom);
 255       for(i=0; (i<nwat); i++) {
 256         sa = ind_sol[na*i];
 257         pbc_dx(&pbc,xcom,xsol[i],dx);
 258         order[i].i   = sa;
 259         order[i].d2  = norm2(dx);
 260       }
 261     } else {
 262       /* Set distance to first atom */
 263       for(i=0; (i<nwat); i++) {
 264         sa = ind_sol[na*i];
 265         pbc_dx(&pbc,x[ind_ref[0]],xsol[i],dx);
 266         order[i].i   = sa;
 267         order[i].d2  = norm2(dx);
 268       }
 269       for(j=1; (j<isize_ref); j++) {
 270         sr = ind_ref[j];
 271         for(i=0; (i<nwat); i++) {
 272           sa = ind_sol[na*i];
 273           pbc_dx(&pbc,x[sr],xsol[i],dx);
 274           n2 = norm2(dx);
 275           if (n2 < order[i].d2)
 276             order[i].d2  = n2;
 277         }
 278       }
 279     }
 280
 281     if (bNShell) {
 282       ncut = 0;
 283       for(i=0; (i<nwat); i++)
 284         if (order[i].d2 <= rcut2)
 285           ncut++;
 286       fprintf(fp,"%10.3f  %8d\n",t,ncut);
 287     }
 288     if (out != -1) {
 289       qsort(order,nwat,sizeof(*order),ocomp);
 290       for(i=0; (i<nwat); i++)
 291         for(j=0; (j<na); j++)
 292           swi[ind_sol[na*i]+j] = order[i].i+j;
 293
 294       /* Store the distance as the B-factor */
 295       if (bPDBout) {
 296         for(i=0; (i<nwat); i++) {
 297           for(j=0; (j<na); j++) {
 298             top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
 299           }
 300         }
 301       }
 302       write_trx(out,natoms,swi,&top.atoms,0,t,box,x,NULL,NULL);
 303     }
 304   } while(read_next_x(oenv,status,&t,natoms,x,box));
 305   close_trj(status);
 306   if (out != -1)
 307     close_trx(out);
 308   if (fp)
 309     ffclose(fp);
 310
 311   thanx(stderr);
 312
 313   return 0;
 314 }