src/gromacs/gmxpreprocess/readrot.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,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include <string>
  41
  42 #include "gromacs/fileio/readinp.h"
  43 #include "gromacs/fileio/trrio.h"
  44 #include "gromacs/fileio/warninp.h"
  45 #include "gromacs/gmxpreprocess/readir.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/math/vecdump.h"
  48 #include "gromacs/mdtypes/inputrec.h"
  49 #include "gromacs/mdtypes/md_enums.h"
  50 #include "gromacs/topology/block.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/utility/path.h"
  56 #include "gromacs/utility/smalloc.h"
  57 #include "gromacs/utility/stringutil.h"
  58
  59 static const std::string RotStr = "Enforced rotation:";
  60
  61 static void string2dvec(char buf[], dvec nums)
  62 {
  63     if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3)
  64     {
  65         gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
  66     }
  67 }
  68
  69
  70 extern std::vector<std::string> read_rotparams(std::vector<t_inpfile>* inp, t_rot* rot, warninp_t wi)
  71 {
  72     int       g, m;
  73     char      buf[STRLEN];
  74     char      warn_buf[STRLEN];
  75     dvec      vec;
  76     t_rotgrp* rotg;
  77
  78     /* read rotation parameters */
  79     printStringNoNewline(
  80             inp,
  81             "Output frequency for angle, torque and rotation potential energy for the whole group");
  82     rot->nstrout = get_eint(inp, "rot-nstrout", 100, wi);
  83     printStringNoNewline(inp,
  84                          "Output frequency for per-slab data (angles, torques and slab centers)");
  85     rot->nstsout = get_eint(inp, "rot-nstsout", 1000, wi);
  86     printStringNoNewline(inp, "Number of rotation groups");
  87     rot->ngrp = get_eint(inp, "rot-ngroups", 1, wi);
  88
  89     if (rot->ngrp < 1)
  90     {
  91         gmx_fatal(FARGS, "rot-ngroups should be >= 1");
  92     }
  93
  94     snew(rot->grp, rot->ngrp);
  95
  96     /* Read the rotation groups */
  97     std::vector<std::string> rotateGroups(rot->ngrp);
  98     char                     readBuffer[STRLEN];
  99     char                     s_vec[STRLEN];
 100     for (g = 0; g < rot->ngrp; g++)
 101     {
 102         rotg = &rot->grp[g];
 103         printStringNoNewline(inp, "Rotation group name");
 104         sprintf(buf, "rot-group%d", g);
 105         setStringEntry(inp, buf, readBuffer, "");
 106         rotateGroups[g] = readBuffer;
 107         printStringNoNewline(inp,
 108                              "Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, "
 109                              "rm2-pf, flex, flex-t, flex2, flex2-t");
 110         sprintf(buf, "rot-type%d", g);
 111         rotg->eType = getEnum<EnforcedRotationGroupType>(inp, buf, wi);
 112
 113         printStringNoNewline(inp, "Use mass-weighting of the rotation group positions");
 114         sprintf(buf, "rot-massw%d", g);
 115         rotg->bMassW = getEnum<Boolean>(inp, buf, wi) != Boolean::No;
 116
 117         printStringNoNewline(inp, "Rotation vector, will get normalized");
 118         sprintf(buf, "rot-vec%d", g);
 119         setStringEntry(inp, buf, s_vec, "1.0 0.0 0.0");
 120         string2dvec(s_vec, vec);
 121         /* Normalize the rotation vector */
 122         if (dnorm(vec) != 0)
 123         {
 124             dsvmul(1.0 / dnorm(vec), vec, vec);
 125         }
 126         else
 127         {
 128             sprintf(warn_buf, "rot-vec%d = 0", g);
 129             warning_error(wi, warn_buf);
 130         }
 131         fprintf(stderr,
 132                 "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
 133                 RotStr.c_str(),
 134                 g,
 135                 enumValueToString(rotg->eType),
 136                 vec[0],
 137                 vec[1],
 138                 vec[2]);
 139         for (m = 0; m < DIM; m++)
 140         {
 141             rotg->inputVec[m] = vec[m];
 142         }
 143
 144         printStringNoNewline(inp, "Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
 145         sprintf(buf, "rot-pivot%d", g);
 146         setStringEntry(inp, buf, s_vec, "0.0 0.0 0.0");
 147         clear_dvec(vec);
 148         if ((rotg->eType == EnforcedRotationGroupType::Iso)
 149             || (rotg->eType == EnforcedRotationGroupType::Pm)
 150             || (rotg->eType == EnforcedRotationGroupType::Rm)
 151             || (rotg->eType == EnforcedRotationGroupType::Rm2))
 152         {
 153             string2dvec(s_vec, vec);
 154         }
 155         for (m = 0; m < DIM; m++)
 156         {
 157             rotg->pivot[m] = vec[m];
 158         }
 159
 160         printStringNoNewline(inp, "Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
 161         sprintf(buf, "rot-rate%d", g);
 162         rotg->rate = get_ereal(inp, buf, 0.0, wi);
 163
 164         sprintf(buf, "rot-k%d", g);
 165         rotg->k = get_ereal(inp, buf, 0.0, wi);
 166         if (rotg->k <= 0.0)
 167         {
 168             sprintf(warn_buf, "rot-k%d <= 0", g);
 169             warning_note(wi, warn_buf);
 170         }
 171
 172         printStringNoNewline(inp, "Slab distance for flexible axis rotation (nm)");
 173         sprintf(buf, "rot-slab-dist%d", g);
 174         rotg->slab_dist = get_ereal(inp, buf, 1.5, wi);
 175         if (rotg->slab_dist <= 0.0)
 176         {
 177             sprintf(warn_buf, "rot-slab-dist%d <= 0", g);
 178             warning_error(wi, warn_buf);
 179         }
 180
 181         printStringNoNewline(inp,
 182                              "Minimum value of Gaussian function for the force to be evaluated "
 183                              "(for flex* potentials)");
 184         sprintf(buf, "rot-min-gauss%d", g);
 185         rotg->min_gaussian = get_ereal(inp, buf, 1e-3, wi);
 186         if (rotg->min_gaussian <= 0.0)
 187         {
 188             sprintf(warn_buf, "rot-min-gauss%d <= 0", g);
 189             warning_error(wi, warn_buf);
 190         }
 191
 192         printStringNoNewline(
 193                 inp, "Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
 194         sprintf(buf, "rot-eps%d", g);
 195         rotg->eps = get_ereal(inp, buf, 1e-4, wi);
 196         if ((rotg->eps <= 0.0)
 197             && (rotg->eType == EnforcedRotationGroupType::Rm2
 198                 || rotg->eType == EnforcedRotationGroupType::Flex2))
 199         {
 200             sprintf(warn_buf, "rot-eps%d <= 0", g);
 201             warning_error(wi, warn_buf);
 202         }
 203
 204         printStringNoNewline(
 205                 inp,
 206                 "Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
 207         sprintf(buf, "rot-fit-method%d", g);
 208         rotg->eFittype = getEnum<RotationGroupFitting>(inp, buf, wi);
 209         printStringNoNewline(inp,
 210                              "For fit type 'potential', nr. of angles around the reference for "
 211                              "which the pot. is evaluated");
 212         sprintf(buf, "rot-potfit-nsteps%d", g);
 213         rotg->PotAngle_nstep = get_eint(inp, buf, 21, wi);
 214         if ((rotg->eFittype == RotationGroupFitting::Pot) && (rotg->PotAngle_nstep < 1))
 215         {
 216             sprintf(warn_buf, "rot-potfit-nsteps%d < 1", g);
 217             warning_error(wi, warn_buf);
 218         }
 219         printStringNoNewline(
 220                 inp, "For fit type 'potential', distance in degrees between two consecutive angles");
 221         sprintf(buf, "rot-potfit-step%d", g);
 222         rotg->PotAngle_step = get_ereal(inp, buf, 0.25, wi);
 223     }
 224
 225     return rotateGroups;
 226 }
 227
 228
 229 /* Check whether the box is unchanged */
 230 static void check_box_unchanged(matrix f_box, matrix box, const char fn[], warninp_t wi)
 231 {
 232     int  i, ii;
 233     bool bSame = TRUE;
 234     char warn_buf[STRLEN];
 235
 236
 237     for (i = 0; i < DIM; i++)
 238     {
 239         for (ii = 0; ii < DIM; ii++)
 240         {
 241             if (f_box[i][ii] != box[i][ii])
 242             {
 243                 bSame = FALSE;
 244             }
 245         }
 246     }
 247     if (!bSame)
 248     {
 249         sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!", RotStr.c_str(), fn);
 250         warning(wi, warn_buf);
 251         pr_rvecs(stderr, 0, "Your box is:", box, 3);
 252         pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
 253     }
 254 }
 255
 256
 257 /* Extract the reference positions for the rotation group(s) */
 258 extern void set_reference_positions(t_rot* rot, rvec* x, matrix box, const char* fn, bool bSet, warninp_t wi)
 259 {
 260     int              g, i, ii;
 261     t_rotgrp*        rotg;
 262     gmx_trr_header_t header;   /* Header information of reference file */
 263     rvec             f_box[3]; /* Box from reference file */
 264
 265     for (g = 0; g < rot->ngrp; g++)
 266     {
 267         rotg = &rot->grp[g];
 268         fprintf(stderr, "%s group %d has %d reference positions.\n", RotStr.c_str(), g, rotg->nat);
 269         snew(rotg->x_ref, rotg->nat);
 270
 271         /* Construct the name for the file containing the reference positions for this group: */
 272         std::string reffileString =
 273                 gmx::Path::concatenateBeforeExtension(fn, gmx::formatString(".%d", g));
 274         const char* reffile = reffileString.c_str();
 275
 276         /* If the base filename for the reference position files was explicitly set by
 277          * the user, we issue a fatal error if the group file can not be found */
 278         if (bSet && !gmx_fexist(reffile))
 279         {
 280             gmx_fatal(FARGS,
 281                       "%s The file containing the reference positions was not found.\n"
 282                       "Expected the file '%s' for group %d.\n",
 283                       RotStr.c_str(),
 284                       reffile,
 285                       g);
 286         }
 287
 288         if (gmx_fexist(reffile))
 289         {
 290             fprintf(stderr, "  Reading them from %s.\n", reffile);
 291             gmx_trr_read_single_header(reffile, &header);
 292             if (rotg->nat != header.natoms)
 293             {
 294                 gmx_fatal(FARGS,
 295                           "Number of atoms in file %s (%d) does not match the number of atoms in "
 296                           "rotation group (%d)!\n",
 297                           reffile,
 298                           header.natoms,
 299                           rotg->nat);
 300             }
 301             gmx_trr_read_single_frame(
 302                     reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, nullptr, nullptr);
 303
 304             /* Check whether the box is unchanged and output a warning if not: */
 305             check_box_unchanged(f_box, box, reffile, wi);
 306         }
 307         else
 308         {
 309             fprintf(stderr, " Saving them to %s.\n", reffile);
 310             for (i = 0; i < rotg->nat; i++)
 311             {
 312                 ii = rotg->ind[i];
 313                 copy_rvec(x[ii], rotg->x_ref[i]);
 314             }
 315             gmx_trr_write_single_frame(reffile, g, 0.0, 0.0, box, rotg->nat, rotg->x_ref, nullptr, nullptr);
 316         }
 317     }
 318 }
 319
 320
 321 extern void make_rotation_groups(t_rot*                           rot,
 322                                  gmx::ArrayRef<const std::string> rotateGroupNames,
 323                                  t_blocka*                        grps,
 324                                  char**                           gnames)
 325 {
 326     int       g, ig = -1, i;
 327     t_rotgrp* rotg;
 328
 329
 330     for (g = 0; g < rot->ngrp; g++)
 331     {
 332         rotg      = &rot->grp[g];
 333         ig        = search_string(rotateGroupNames[g].c_str(), grps->nr, gnames);
 334         rotg->nat = grps->index[ig + 1] - grps->index[ig];
 335
 336         if (rotg->nat > 0)
 337         {
 338             fprintf(stderr, "Rotation group %d '%s' has %d atoms\n", g, rotateGroupNames[g].c_str(), rotg->nat);
 339             snew(rotg->ind, rotg->nat);
 340             for (i = 0; i < rotg->nat; i++)
 341             {
 342                 rotg->ind[i] = grps->a[grps->index[ig] + i];
 343             }
 344         }
 345         else
 346         {
 347             gmx_fatal(FARGS, "Rotation group %d '%s' is empty", g, rotateGroupNames[g].c_str());
 348         }
 349     }
 350 }