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