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