src/kernel/g_protonate.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 "string2.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "copyrite.h"
  47 #include "smalloc.h"
  48 #include "statutil.h"
  49 #include "confio.h"
  50 #include "genhydro.h"
  51 #include "tpxio.h"
  52 #include "index.h"
  53 #include "vec.h"
  54 #include "hackblock.h"
  55
  56 int cmain (int argc,char *argv[])
  57 {
  58   const char *desc[] = {
  59     "[TT]g_protonate[tt] reads (a) conformation(s) and adds all missing",
  60     "hydrogens as defined in [TT]gmx2.ff/aminoacids.hdb[tt]. If only [TT]-s[tt] is",
  61     "specified, this conformation will be protonated, if also [TT]-f[tt]",
  62     "is specified, the conformation(s) will be read from this file, ",
  63     "which can be either a single conformation or a trajectory.",
  64     "[PAR]",
  65     "If a [TT].pdb[tt] file is supplied, residue names might not correspond to",
  66     "to the GROMACS naming conventions, in which case these residues will",
  67     "probably not be properly protonated.",
  68     "[PAR]",
  69     "If an index file is specified, please note that the atom numbers",
  70     "should correspond to the [BB]protonated[bb] state."
  71   };
  72
  73   char        title[STRLEN+1];
  74   const char  *infile;
  75   char        *grpnm;
  76   t_topology  top;
  77   int         ePBC;
  78   t_atoms     *atoms,*iatoms;
  79   t_protonate protdata;
  80   atom_id     *index;
  81   t_trxstatus *status;
  82   t_trxstatus *out;
  83   t_trxframe  fr,frout;
  84   rvec        *x,*ix;
  85   int         nidx,natoms,natoms_out;
  86   matrix      box;
  87   int         i,frame,resind;
  88   gmx_bool        bReadMultiple;
  89   output_env_t oenv;
  90
  91   const char *bugs[] = {
  92     "For the moment, only .pdb files are accepted to the -s flag"
  93   };
  94
  95   t_filenm fnm[] = {
  96     { efTPS, NULL, NULL,         ffREAD  },
  97     { efTRX, "-f", NULL,         ffOPTRD },
  98     { efNDX, NULL, NULL,         ffOPTRD },
  99     { efTRO, "-o", "protonated", ffWRITE }
 100   };
 101 #define NFILE asize(fnm)
 102
 103   CopyRight(stderr,argv[0]);
 104   parse_common_args(&argc,argv,PCA_CAN_TIME,
 105                     NFILE,fnm,0,NULL,asize(desc),desc,asize(bugs),bugs,&oenv);
 106
 107   infile=opt2fn("-s",NFILE,fnm);
 108   read_tps_conf(infile,title,&top,&ePBC,&x,NULL,box,FALSE);
 109   atoms=&(top.atoms);
 110   printf("Select group to process:\n");
 111   get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nidx,&index,&grpnm);
 112   bReadMultiple = opt2bSet("-f",NFILE,fnm);
 113   if (bReadMultiple) {
 114     infile = opt2fn("-f",NFILE,fnm);
 115     if ( !read_first_frame(oenv,&status, infile, &fr, TRX_NEED_X ) ) {
 116       gmx_fatal(FARGS,"cannot read coordinate file %s",infile);
 117     }
 118     natoms = fr.natoms;
 119   } else {
 120     clear_trxframe(&fr,TRUE);
 121     fr.natoms = atoms->nr;
 122     fr.bTitle = TRUE;
 123     fr.title  = title;
 124     fr.bX     = TRUE;
 125     fr.x      = x;
 126     fr.bBox   = TRUE;
 127     copy_mat(box, fr.box);
 128     natoms = fr.natoms;
 129   }
 130
 131   /* check input */
 132   if ( natoms == 0 ) {
 133     gmx_fatal(FARGS,"no atoms in coordinate file %s",infile);
 134   }
 135
 136   if ( natoms > atoms->nr ) {
 137     gmx_fatal(FARGS,"topology with %d atoms does not match "
 138                 "coordinates with %d atoms",atoms->nr,natoms);
 139   }
 140
 141   for(i=0; i<nidx; i++) {
 142     if (index[i] > natoms) {
 143       gmx_fatal(FARGS,"An atom number in group %s is larger than the number of "
 144                   "atoms (%d) in the coordinate file %s",grpnm,natoms,infile);
 145     }
 146   }
 147
 148   /* get indexed copy of atoms */
 149   snew(iatoms,1);
 150   init_t_atoms(iatoms,nidx,FALSE);
 151   snew(iatoms->atom, iatoms->nr);
 152   resind = 0;
 153   for(i=0; i<nidx; i++) {
 154     iatoms->atom[i] = atoms->atom[index[i]];
 155     iatoms->atomname[i] = atoms->atomname[index[i]];
 156     if ( i>0 && (atoms->atom[index[i]].resind!=atoms->atom[index[i-1]].resind) ) {
 157       resind++;
 158     }
 159     iatoms->atom[i].resind = resind;
 160     iatoms->resinfo[resind] = atoms->resinfo[atoms->atom[index[i]].resind];
 161     /* allocate some space for the rtp name and copy from name */
 162     snew(iatoms->resinfo[resind].rtp,1);
 163     *iatoms->resinfo[resind].rtp = gmx_strdup(*atoms->resinfo[resind].name);
 164
 165     iatoms->nres = max(iatoms->nres, iatoms->atom[i].resind+1);
 166   }
 167
 168   init_t_protonate(&protdata);
 169
 170   out = open_trx(opt2fn("-o",NFILE,fnm),"w");
 171   snew(ix, nidx);
 172   frame=0;
 173   do {
 174     if (debug) {
 175         fprintf(debug,"FRAME %d (%d %g)\n",frame,fr.step,fr.time);
 176     }
 177     /* get indexed copy of x */
 178     for(i=0; i<nidx; i++) {
 179       copy_rvec(fr.x[index[i]], ix[i]);
 180     }
 181     /* protonate */
 182     natoms_out = protonate(&iatoms, &ix, &protdata);
 183
 184     /* setup output frame */
 185     frout = fr;
 186     frout.natoms = natoms_out;
 187     frout.bAtoms = TRUE;
 188     frout.atoms  = iatoms;
 189     frout.bV     = FALSE;
 190     frout.bF     = FALSE;
 191     frout.x      = ix;
 192
 193     /* write output */
 194     write_trxframe(out,&frout,NULL);
 195     frame++;
 196   } while ( bReadMultiple && read_next_frame(oenv,status, &fr) );
 197
 198   sfree(ix);
 199   sfree(iatoms);
 200
 201   thanx(stderr);
 202
 203   return 0;
 204 }
 205