src/gromacs/gmxpreprocess/genrestr.cpp

   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,2015,2017,2018,2019, 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 "genrestr.h"
  40
  41 #include <cmath>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/topology/index.h"
  50 #include "gromacs/topology/mtop_util.h"
  51 #include "gromacs/topology/topology.h"
  52 #include "gromacs/utility/arraysize.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 int gmx_genrestr(int argc, char* argv[])
  59 {
  60     const char* desc[] = {
  61         "[THISMODULE] produces an #include file for a topology containing",
  62         "a list of atom numbers and three force constants for the",
  63         "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-direction based on",
  64         "the contents of the [TT]-f[tt] file. A single isotropic force constant may",
  65         "be given on the command line instead of three components.[PAR]",
  66         "WARNING: Position restraints are interactions within molecules, therefore",
  67         "they must be included within the correct [TT][ moleculetype ][tt]",
  68         "block in the topology. The atom indices within the",
  69         "[TT][ position_restraints ][tt] block must be within the range of the",
  70         "atom indices for that molecule type. Since the atom numbers in every",
  71         "moleculetype in the topology start at 1 and the numbers in the input file",
  72         "for [THISMODULE] number consecutively from 1, [THISMODULE] will only",
  73         "produce a useful file for the first molecule. You may wish to",
  74         "edit the resulting index file to remove the lines for later atoms,",
  75         "or construct a suitable index group to provide",
  76         "as input to [THISMODULE].[PAR]",
  77         "The [TT]-of[tt] option produces an index file that can be used for",
  78         "freezing atoms. In this case, the input file must be a [REF].pdb[ref] file.[PAR]",
  79         "With the [TT]-disre[tt] option, half a matrix of distance restraints",
  80         "is generated instead of position restraints. With this matrix, that",
  81         "one typically would apply to C[GRK]alpha[grk] atoms in a protein, one can",
  82         "maintain the overall conformation of a protein without tieing it to",
  83         "a specific position (as with position restraints)."
  84     };
  85     static rvec     fc           = { 1000.0, 1000.0, 1000.0 };
  86     static real     freeze_level = 0.0;
  87     static real     disre_dist   = 0.1;
  88     static real     disre_frac   = 0.0;
  89     static real     disre_up2    = 1.0;
  90     static gmx_bool bDisre       = FALSE;
  91     static gmx_bool bConstr      = FALSE;
  92     static real     cutoff       = -1.0;
  93
  94     t_pargs pa[] = {
  95         { "-fc", FALSE, etRVEC, { fc }, "Force constants (kJ/mol nm^2)" },
  96         { "-freeze",
  97           FALSE,
  98           etREAL,
  99           { &freeze_level },
 100           "If the [TT]-of[tt] option or this one is given an index file will be written containing "
 101           "atom numbers of all atoms that have a B-factor less than the level given here" },
 102         { "-disre",
 103           FALSE,
 104           etBOOL,
 105           { &bDisre },
 106           "Generate a distance restraint matrix for all the atoms in index" },
 107         { "-disre_dist",
 108           FALSE,
 109           etREAL,
 110           { &disre_dist },
 111           "Distance range around the actual distance for generating distance restraints" },
 112         { "-disre_frac",
 113           FALSE,
 114           etREAL,
 115           { &disre_frac },
 116           "Fraction of distance to be used as interval rather than a fixed distance. If the "
 117           "fraction of the distance that you specify here is less than the distance given in the "
 118           "previous option, that one is used instead." },
 119         { "-disre_up2",
 120           FALSE,
 121           etREAL,
 122           { &disre_up2 },
 123           "Distance between upper bound for distance restraints, and the distance at which the "
 124           "force becomes constant (see manual)" },
 125         { "-cutoff",
 126           FALSE,
 127           etREAL,
 128           { &cutoff },
 129           "Only generate distance restraints for atoms pairs within cutoff (nm)" },
 130         { "-constr",
 131           FALSE,
 132           etBOOL,
 133           { &bConstr },
 134           "Generate a constraint matrix rather than distance restraints. Constraints of type 2 "
 135           "will be generated that do generate exclusions." }
 136     };
 137 #define npargs asize(pa)
 138
 139     gmx_output_env_t* oenv;
 140     t_atoms           atoms;
 141     int               i, j, k;
 142     FILE*             out;
 143     int               igrp;
 144     real              d, dd, lo, hi;
 145     int*              ind_grp;
 146     const char *      xfn, *nfn;
 147     char*             gn_grp;
 148     matrix            box;
 149     gmx_bool          bFreeze;
 150     rvec              dx, *x = nullptr, *v = nullptr;
 151
 152     t_filenm fnm[] = { { efSTX, "-f", nullptr, ffREAD },
 153                        { efNDX, "-n", nullptr, ffOPTRD },
 154                        { efITP, "-o", "posre", ffWRITE },
 155                        { efNDX, "-of", "freeze", ffOPTWR } };
 156 #define NFILE asize(fnm)
 157
 158     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, npargs, pa, asize(desc), desc, 0, nullptr, &oenv))
 159     {
 160         return 0;
 161     }
 162
 163     bFreeze = opt2bSet("-of", NFILE, fnm) || opt2parg_bSet("-freeze", asize(pa), pa);
 164     bDisre  = bDisre || opt2parg_bSet("-disre_dist", npargs, pa);
 165     xfn     = opt2fn_null("-f", NFILE, fnm);
 166     nfn     = opt2fn_null("-n", NFILE, fnm);
 167
 168     if ((nfn == nullptr) && (xfn == nullptr))
 169     {
 170         gmx_fatal(FARGS, "no index file and no structure file supplied");
 171     }
 172
 173     if ((disre_frac < 0) || (disre_frac >= 1))
 174     {
 175         gmx_fatal(FARGS, "disre_frac should be between 0 and 1");
 176     }
 177     if (disre_dist < 0)
 178     {
 179         gmx_fatal(FARGS, "disre_dist should be >= 0");
 180     }
 181
 182     const char* title        = "";
 183     bool        haveTopology = false;
 184     if (xfn != nullptr)
 185     {
 186         fprintf(stderr, "\nReading structure file\n");
 187         gmx_mtop_t mtop;
 188         readConfAndTopology(xfn, &haveTopology, &mtop, nullptr, &x, &v, box);
 189         title = *mtop.name;
 190         atoms = gmx_mtop_global_atoms(&mtop);
 191         if (atoms.pdbinfo == nullptr)
 192         {
 193             snew(atoms.pdbinfo, atoms.nr);
 194         }
 195         haveTopology = true;
 196     }
 197
 198     if (bFreeze)
 199     {
 200         if (!haveTopology || !atoms.pdbinfo)
 201         {
 202             gmx_fatal(FARGS, "No B-factors in input file %s, use a pdb file next time.", xfn);
 203         }
 204
 205         out = opt2FILE("-of", NFILE, fnm, "w");
 206         fprintf(out, "[ freeze ]\n");
 207         for (i = 0; (i < atoms.nr); i++)
 208         {
 209             if (atoms.pdbinfo[i].bfac <= freeze_level)
 210             {
 211                 fprintf(out, "%d\n", i + 1);
 212             }
 213         }
 214         gmx_ffclose(out);
 215     }
 216     else if ((bDisre || bConstr) && x)
 217     {
 218         printf("Select group to generate %s matrix from\n", bConstr ? "constraint" : "distance restraint");
 219         get_index(&atoms, nfn, 1, &igrp, &ind_grp, &gn_grp);
 220
 221         out = ftp2FILE(efITP, NFILE, fnm, "w");
 222         if (bConstr)
 223         {
 224             fprintf(out, "; constraints for %s of %s\n\n", gn_grp, title);
 225             fprintf(out, "[ constraints ]\n");
 226             fprintf(out, ";%4s %5s %1s %10s\n", "i", "j", "tp", "dist");
 227         }
 228         else
 229         {
 230             fprintf(out, "; distance restraints for %s of %s\n\n", gn_grp, title);
 231             fprintf(out, "[ distance_restraints ]\n");
 232             fprintf(out, ";%4s %5s %1s %5s %10s %10s %10s %10s %10s\n", "i", "j", "?", "label",
 233                     "funct", "lo", "up1", "up2", "weight");
 234         }
 235         for (i = k = 0; i < igrp; i++)
 236         {
 237             for (j = i + 1; j < igrp; j++, k++)
 238             {
 239                 rvec_sub(x[ind_grp[i]], x[ind_grp[j]], dx);
 240                 d = norm(dx);
 241                 if (bConstr)
 242                 {
 243                     fprintf(out, "%5d %5d %1d %10g\n", ind_grp[i] + 1, ind_grp[j] + 1, 2, d);
 244                 }
 245                 else
 246                 {
 247                     if (cutoff < 0 || d < cutoff)
 248                     {
 249                         if (disre_frac > 0)
 250                         {
 251                             dd = std::min(disre_dist, disre_frac * d);
 252                         }
 253                         else
 254                         {
 255                             dd = disre_dist;
 256                         }
 257                         lo = std::max(0.0_real, d - dd);
 258                         hi = d + dd;
 259                         fprintf(out, "%5d %5d %1d %5d %10d %10g %10g %10g %10g\n", ind_grp[i] + 1,
 260                                 ind_grp[j] + 1, 1, k, 1, lo, hi, hi + disre_up2, 1.0);
 261                     }
 262                 }
 263             }
 264         }
 265         gmx_ffclose(out);
 266     }
 267     else
 268     {
 269         printf("Select group to position restrain\n");
 270         get_index(&atoms, nfn, 1, &igrp, &ind_grp, &gn_grp);
 271
 272         out = ftp2FILE(efITP, NFILE, fnm, "w");
 273         fprintf(out, "; position restraints for %s of %s\n\n", gn_grp, title);
 274         fprintf(out, "[ position_restraints ]\n");
 275         fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
 276         for (i = 0; i < igrp; i++)
 277         {
 278             fprintf(out, "%4d %4d %10g %10g %10g\n", ind_grp[i] + 1, 1, fc[XX], fc[YY], fc[ZZ]);
 279         }
 280         gmx_ffclose(out);
 281     }
 282     if (xfn)
 283     {
 284         sfree(x);
 285         sfree(v);
 286         done_atom(&atoms);
 287     }
 288
 289     return 0;
 290 }