src/gromacs/gmxpreprocess/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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "protonate.h"
  40
  41 #include <math.h>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/gmxpreprocess/genhydro.h"
  48 #include "gromacs/gmxpreprocess/hackblock.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/legacyheaders/typedefs.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 int gmx_protonate(int argc, char *argv[])
  58 {
  59     const char     *desc[] = {
  60         "[THISMODULE] reads (a) conformation(s) and adds all missing",
  61         "hydrogens as defined in [TT]oplsaa.ff/aminoacids.hdb[tt]. If only [TT]-s[tt] is",
  62         "specified, this conformation will be protonated, if also [TT]-f[tt]",
  63         "is specified, the conformation(s) will be read from this file, ",
  64         "which can be either a single conformation or a trajectory.",
  65         "[PAR]",
  66         "If a [TT].pdb[tt] file is supplied, residue names might not correspond to",
  67         "to the GROMACS naming conventions, in which case these residues will",
  68         "probably not be properly protonated.",
  69         "[PAR]",
  70         "If an index file is specified, please note that the atom numbers",
  71         "should correspond to the [BB]protonated[bb] state."
  72     };
  73
  74     char            title[STRLEN+1];
  75     const char     *infile;
  76     char           *grpnm;
  77     t_topology      top;
  78     int             ePBC;
  79     t_atoms        *atoms, *iatoms;
  80     t_protonate     protdata;
  81     atom_id        *index;
  82     t_trxstatus    *status;
  83     t_trxstatus    *out;
  84     t_trxframe      fr, frout;
  85     rvec           *x, *ix;
  86     int             nidx, natoms, natoms_out;
  87     matrix          box;
  88     int             i, frame, resind;
  89     gmx_bool        bReadMultiple;
  90     output_env_t    oenv;
  91
  92     const char     *bugs[] = {
  93         "For the moment, only .pdb files are accepted to the -s flag"
  94     };
  95
  96     t_filenm        fnm[] = {
  97         { efTPS, NULL, NULL,         ffREAD  },
  98         { efTRX, "-f", NULL,         ffOPTRD },
  99         { efNDX, NULL, NULL,         ffOPTRD },
 100         { efTRO, "-o", "protonated", ffWRITE }
 101     };
 102 #define NFILE asize(fnm)
 103
 104     if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
 105                            NFILE, fnm, 0, NULL, asize(desc), desc, asize(bugs), bugs, &oenv))
 106     {
 107         return 0;
 108     }
 109
 110     infile = opt2fn("-s", NFILE, fnm);
 111     read_tps_conf(infile, title, &top, &ePBC, &x, NULL, box, FALSE);
 112     atoms = &(top.atoms);
 113     printf("Select group to process:\n");
 114     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
 115     bReadMultiple = opt2bSet("-f", NFILE, fnm);
 116     if (bReadMultiple)
 117     {
 118         infile = opt2fn("-f", NFILE, fnm);
 119         if (!read_first_frame(oenv, &status, infile, &fr, TRX_NEED_X ) )
 120         {
 121             gmx_fatal(FARGS, "cannot read coordinate file %s", infile);
 122         }
 123         natoms = fr.natoms;
 124     }
 125     else
 126     {
 127         clear_trxframe(&fr, TRUE);
 128         fr.natoms = atoms->nr;
 129         fr.bTitle = TRUE;
 130         fr.title  = title;
 131         fr.bX     = TRUE;
 132         fr.x      = x;
 133         fr.bBox   = TRUE;
 134         copy_mat(box, fr.box);
 135         natoms = fr.natoms;
 136     }
 137
 138     /* check input */
 139     if (natoms == 0)
 140     {
 141         gmx_fatal(FARGS, "no atoms in coordinate file %s", infile);
 142     }
 143
 144     if (natoms > atoms->nr)
 145     {
 146         gmx_fatal(FARGS, "topology with %d atoms does not match "
 147                   "coordinates with %d atoms", atoms->nr, natoms);
 148     }
 149
 150     for (i = 0; i < nidx; i++)
 151     {
 152         if (index[i] > natoms)
 153         {
 154             gmx_fatal(FARGS, "An atom number in group %s is larger than the number of "
 155                       "atoms (%d) in the coordinate file %s", grpnm, natoms, infile);
 156         }
 157     }
 158
 159     /* get indexed copy of atoms */
 160     snew(iatoms, 1);
 161     init_t_atoms(iatoms, nidx, FALSE);
 162     snew(iatoms->atom, iatoms->nr);
 163     resind = 0;
 164     for (i = 0; i < nidx; i++)
 165     {
 166         iatoms->atom[i]     = atoms->atom[index[i]];
 167         iatoms->atomname[i] = atoms->atomname[index[i]];
 168         if (i > 0 && (atoms->atom[index[i]].resind != atoms->atom[index[i-1]].resind) )
 169         {
 170             resind++;
 171         }
 172         iatoms->atom[i].resind  = resind;
 173         iatoms->resinfo[resind] = atoms->resinfo[atoms->atom[index[i]].resind];
 174         /* allocate some space for the rtp name and copy from name */
 175         snew(iatoms->resinfo[resind].rtp, 1);
 176         *iatoms->resinfo[resind].rtp = gmx_strdup(*atoms->resinfo[resind].name);
 177
 178         iatoms->nres = max(iatoms->nres, iatoms->atom[i].resind+1);
 179     }
 180
 181     init_t_protonate(&protdata);
 182
 183     out = open_trx(opt2fn("-o", NFILE, fnm), "w");
 184     snew(ix, nidx);
 185     frame = 0;
 186     do
 187     {
 188         if (debug)
 189         {
 190             fprintf(debug, "FRAME %d (%d %g)\n", frame, fr.step, fr.time);
 191         }
 192         /* get indexed copy of x */
 193         for (i = 0; i < nidx; i++)
 194         {
 195             copy_rvec(fr.x[index[i]], ix[i]);
 196         }
 197         /* protonate */
 198         natoms_out = protonate(&iatoms, &ix, &protdata);
 199
 200         /* setup output frame */
 201         frout        = fr;
 202         frout.natoms = natoms_out;
 203         frout.bAtoms = TRUE;
 204         frout.atoms  = iatoms;
 205         frout.bV     = FALSE;
 206         frout.bF     = FALSE;
 207         frout.x      = ix;
 208
 209         /* write output */
 210         write_trxframe(out, &frout, NULL);
 211         frame++;
 212     }
 213     while (bReadMultiple && read_next_frame(oenv, status, &fr) );
 214
 215     sfree(ix);
 216     sfree(iatoms);
 217
 218     return 0;
 219 }