src/kernel/protonate.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * GROningen Mixture of Alchemy and Childrens' Stories
  31  */
  32 static char *SRCID_protonate_c = "$Id$";
  33 #include <math.h>
  34 #include "typedefs.h"
  35 #include "macros.h"
  36 #include "copyrite.h"
  37 #include "smalloc.h"
  38 #include "statutil.h"
  39 #include "confio.h"
  40 #include "genhydro.h"
  41 #include "tpxio.h"
  42 #include "rdgroup.h"
  43 #include "vec.h"
  44 #include "hackblock.h"
  45
  46 int main (int argc,char *argv[])
  47 {
  48   static char *desc[] = {
  49     "[TT]protonate[tt] reads (a) conformation(s) and adds all missing",
  50     "hydrogens as defined in [TT]ffgmx2.hdb[tt]. If only [TT]-s[tt] is",
  51     "specified, this conformation will be protonated, if also [TT]-f[tt]",
  52     "is specified, the conformation(s) will be read from this file",
  53     "which can be either a single conformation or a trajectory.",
  54     "[PAR]",
  55     "If a pdb file is supplied, residue names might not correspond to",
  56     "to the GROMACS naming conventions, in which case these residues will",
  57     "probably not be properly protonated.",
  58     "[PAR]",
  59     "If an index file is specified, please note that the atom numbers",
  60     "should correspond to the [BB]protonated[bb] state."
  61   };
  62
  63   char        title[STRLEN+1];
  64   char        *infile;
  65   char        *grpnm;
  66   t_topology  top;
  67   t_atoms     *atoms,*iatoms;
  68   t_protonate protdata;
  69   atom_id     *index;
  70   int         status,out;
  71   t_trxframe  fr,frout;
  72   rvec        *x,*ix;
  73   int         nidx,natoms,natoms_out;
  74   matrix      box;
  75   int         i,frame,resnr;
  76   bool        bReadMultiple;
  77
  78   t_filenm fnm[] = {
  79     { efTPS, NULL, NULL,         ffREAD  },
  80     { efTRX, "-f", NULL,         ffOPTRD },
  81     { efNDX, NULL, NULL,         ffOPTRD },
  82     { efTRX, "-o", "protonated", ffWRITE }
  83   };
  84 #define NFILE asize(fnm)
  85
  86   CopyRight(stderr,argv[0]);
  87   parse_common_args(&argc,argv,PCA_CAN_TIME,
  88                     NFILE,fnm,0,NULL,asize(desc),desc,0,NULL);
  89
  90   infile=opt2fn("-s",NFILE,fnm);
  91   read_tps_conf(infile,title,&top,&x,NULL,box,FALSE);
  92   atoms=&(top.atoms);
  93   printf("Select group to process:\n");
  94   get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nidx,&index,&grpnm);
  95   bReadMultiple = opt2bSet("-f",NFILE,fnm);
  96   if (bReadMultiple) {
  97     infile = opt2fn("-f",NFILE,fnm);
  98     if ( !read_first_frame(&status, infile, &fr, TRX_NEED_X ) )
  99       fatal_error(0,"cannot read coordinate file %s",infile);
 100     natoms = fr.natoms;
 101   } else {
 102     clear_trxframe(&fr,TRUE);
 103     fr.natoms = atoms->nr;
 104     fr.bTitle = TRUE;
 105     fr.title  = title;
 106     fr.bX     = TRUE;
 107     fr.x      = x;
 108     fr.bBox   = TRUE;
 109     copy_mat(box, fr.box);
 110     natoms = fr.natoms;
 111   }
 112
 113   /* check input */
 114   if ( natoms == 0 )
 115     fatal_error(0,"no atoms in coordinate file %s",infile);
 116   if ( natoms > atoms->nr )
 117     fatal_error(0,"topology with %d atoms does not match "
 118                 "coordinates with %d atoms",atoms->nr,natoms);
 119   for(i=0; i<nidx; i++)
 120     if (index[i] > natoms)
 121       fatal_error(0,"An atom number in group %s is larger than the number of "
 122                   "atoms (%d) in the coordinate file %s",grpnm,natoms,infile);
 123
 124   /* get indexed copy of atoms */
 125   snew(iatoms,1);
 126   init_t_atoms(iatoms,nidx,FALSE);
 127   snew(iatoms->atom, iatoms->nr);
 128   resnr = 0;
 129   for(i=0; i<nidx; i++) {
 130     iatoms->atom[i] = atoms->atom[index[i]];
 131     iatoms->atomname[i] = atoms->atomname[index[i]];
 132     if ( i>0 && (atoms->atom[index[i]].resnr!=atoms->atom[index[i-1]].resnr) )
 133       resnr++;
 134     iatoms->atom[i].resnr = resnr;
 135     iatoms->resname[resnr] = atoms->resname[atoms->atom[index[i]].resnr];
 136     iatoms->nres = max(iatoms->nres, iatoms->atom[i].resnr+1);
 137   }
 138
 139   init_t_protonate(&protdata);
 140
 141   out = open_trx(opt2fn("-o",NFILE,fnm),"w");
 142   snew(ix, nidx);
 143   frame=0;
 144   do {
 145     if (debug) fprintf(debug,"FRAME %d (%d %g)\n",frame,fr.step,fr.time);
 146     /* get indexed copy of x */
 147     for(i=0; i<nidx; i++)
 148       copy_rvec(fr.x[index[i]], ix[i]);
 149     /* protonate */
 150     natoms_out = protonate(&iatoms, &ix, &protdata);
 151
 152     /* setup output frame */
 153     frout = fr;
 154     frout.natoms = natoms_out;
 155     frout.bAtoms = TRUE;
 156     frout.atoms  = iatoms;
 157     frout.x      = ix;
 158
 159     /* write output */
 160     write_trxframe(out,&frout);
 161     frame++;
 162   } while ( bReadMultiple && read_next_frame(status, &fr) );
 163
 164   thanx(stderr);
 165
 166   return 0;
 167 }
 168