src/programs/gmx/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 "smalloc.h"
  44 #include "statutil.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "genhydro.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "index.h"
  50 #include "vec.h"
  51 #include "hackblock.h"
  52
  53 int gmx_protonate(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     if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
 101                            NFILE, fnm, 0, NULL, asize(desc), desc, asize(bugs), bugs, &oenv))
 102     {
 103         return 0;
 104     }
 105
 106     infile = opt2fn("-s", NFILE, fnm);
 107     read_tps_conf(infile, title, &top, &ePBC, &x, NULL, box, FALSE);
 108     atoms = &(top.atoms);
 109     printf("Select group to process:\n");
 110     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
 111     bReadMultiple = opt2bSet("-f", NFILE, fnm);
 112     if (bReadMultiple)
 113     {
 114         infile = opt2fn("-f", NFILE, fnm);
 115         if (!read_first_frame(oenv, &status, infile, &fr, TRX_NEED_X ) )
 116         {
 117             gmx_fatal(FARGS, "cannot read coordinate file %s", infile);
 118         }
 119         natoms = fr.natoms;
 120     }
 121     else
 122     {
 123         clear_trxframe(&fr, TRUE);
 124         fr.natoms = atoms->nr;
 125         fr.bTitle = TRUE;
 126         fr.title  = title;
 127         fr.bX     = TRUE;
 128         fr.x      = x;
 129         fr.bBox   = TRUE;
 130         copy_mat(box, fr.box);
 131         natoms = fr.natoms;
 132     }
 133
 134     /* check input */
 135     if (natoms == 0)
 136     {
 137         gmx_fatal(FARGS, "no atoms in coordinate file %s", infile);
 138     }
 139
 140     if (natoms > atoms->nr)
 141     {
 142         gmx_fatal(FARGS, "topology with %d atoms does not match "
 143                   "coordinates with %d atoms", atoms->nr, natoms);
 144     }
 145
 146     for (i = 0; i < nidx; i++)
 147     {
 148         if (index[i] > natoms)
 149         {
 150             gmx_fatal(FARGS, "An atom number in group %s is larger than the number of "
 151                       "atoms (%d) in the coordinate file %s", grpnm, natoms, infile);
 152         }
 153     }
 154
 155     /* get indexed copy of atoms */
 156     snew(iatoms, 1);
 157     init_t_atoms(iatoms, nidx, FALSE);
 158     snew(iatoms->atom, iatoms->nr);
 159     resind = 0;
 160     for (i = 0; i < nidx; i++)
 161     {
 162         iatoms->atom[i]     = atoms->atom[index[i]];
 163         iatoms->atomname[i] = atoms->atomname[index[i]];
 164         if (i > 0 && (atoms->atom[index[i]].resind != atoms->atom[index[i-1]].resind) )
 165         {
 166             resind++;
 167         }
 168         iatoms->atom[i].resind  = resind;
 169         iatoms->resinfo[resind] = atoms->resinfo[atoms->atom[index[i]].resind];
 170         /* allocate some space for the rtp name and copy from name */
 171         snew(iatoms->resinfo[resind].rtp, 1);
 172         *iatoms->resinfo[resind].rtp = gmx_strdup(*atoms->resinfo[resind].name);
 173
 174         iatoms->nres = max(iatoms->nres, iatoms->atom[i].resind+1);
 175     }
 176
 177     init_t_protonate(&protdata);
 178
 179     out = open_trx(opt2fn("-o", NFILE, fnm), "w");
 180     snew(ix, nidx);
 181     frame = 0;
 182     do
 183     {
 184         if (debug)
 185         {
 186             fprintf(debug, "FRAME %d (%d %g)\n", frame, fr.step, fr.time);
 187         }
 188         /* get indexed copy of x */
 189         for (i = 0; i < nidx; i++)
 190         {
 191             copy_rvec(fr.x[index[i]], ix[i]);
 192         }
 193         /* protonate */
 194         natoms_out = protonate(&iatoms, &ix, &protdata);
 195
 196         /* setup output frame */
 197         frout        = fr;
 198         frout.natoms = natoms_out;
 199         frout.bAtoms = TRUE;
 200         frout.atoms  = iatoms;
 201         frout.bV     = FALSE;
 202         frout.bF     = FALSE;
 203         frout.x      = ix;
 204
 205         /* write output */
 206         write_trxframe(out, &frout, NULL);
 207         frame++;
 208     }
 209     while (bReadMultiple && read_next_frame(oenv, status, &fr) );
 210
 211     sfree(ix);
 212     sfree(iatoms);
 213
 214     return 0;
 215 }