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