src/gromacs/gmxpreprocess/readrot.c

   1 /*
   2  *                This source code is part of
   3  *
   4  *                 G   R   O   M   A   C   S
   5  *
   6  *          GROningen MAchine for Chemical Simulations
   7  *
   8  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
   9  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  10  * Copyright (c) 2001-2004, The GROMACS development team,
  11  * check out http://www.gromacs.org for more information.
  12
  13  * This program is free software; you can redistribute it and/or
  14  * modify it under the terms of the GNU General Public License
  15  * as published by the Free Software Foundation; either version 2
  16  * of the License, or (at your option) any later version.
  17  *
  18  * If you want to redistribute modifications, please consider that
  19  * scientific software is very special. Version control is crucial -
  20  * bugs must be traceable. We will be happy to consider code for
  21  * inclusion in the official distribution, but derived work must not
  22  * be called official GROMACS. Details are found in the README & COPYING
  23  * files - if they are missing, get the official version at www.gromacs.org.
  24  *
  25  * To help us fund GROMACS development, we humbly ask that you cite
  26  * the papers on the package - you can find them in the top README file.
  27  *
  28  * For more info, check our website at http://www.gromacs.org
  29  *
  30  * And Hey:
  31  * GROwing Monsters And Cloning Shrimps
  32  */
  33 #ifdef HAVE_CONFIG_H
  34 #include <config.h>
  35 #endif
  36
  37 #include "vec.h"
  38 #include "smalloc.h"
  39 #include "readir.h"
  40 #include "names.h"
  41 #include "futil.h"
  42 #include "trnio.h"
  43 #include "txtdump.h"
  44
  45 static char *RotStr = {"Enforced rotation:"};
  46
  47
  48 static char s_vec[STRLEN];
  49
  50
  51 static void string2dvec(char buf[], dvec nums)
  52 {
  53     if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3)
  54     {
  55         gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
  56     }
  57 }
  58
  59
  60 extern char **read_rotparams(int *ninp_p, t_inpfile **inp_p, t_rot *rot,
  61                              warninp_t wi)
  62 {
  63     int         ninp, g, m;
  64     t_inpfile  *inp;
  65     const char *tmp;
  66     char      **grpbuf;
  67     char        buf[STRLEN];
  68     char        warn_buf[STRLEN];
  69     dvec        vec;
  70     t_rotgrp   *rotg;
  71
  72     ninp   = *ninp_p;
  73     inp    = *inp_p;
  74
  75     /* read rotation parameters */
  76     CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
  77     ITYPE("rot_nstrout",     rot->nstrout, 100);
  78     CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
  79     ITYPE("rot_nstsout",     rot->nstsout, 1000);
  80     CTYPE("Number of rotation groups");
  81     ITYPE("rot_ngroups",     rot->ngrp, 1);
  82
  83     if (rot->ngrp < 1)
  84     {
  85         gmx_fatal(FARGS, "rot_ngroups should be >= 1");
  86     }
  87
  88     snew(rot->grp, rot->ngrp);
  89
  90     /* Read the rotation groups */
  91     snew(grpbuf, rot->ngrp);
  92     for (g = 0; g < rot->ngrp; g++)
  93     {
  94         rotg = &rot->grp[g];
  95         snew(grpbuf[g], STRLEN);
  96         CTYPE("Rotation group name");
  97         sprintf(buf, "rot_group%d", g);
  98         STYPE(buf, grpbuf[g], "");
  99
 100         CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
 101         sprintf(buf, "rot_type%d", g);
 102         ETYPE(buf, rotg->eType, erotg_names);
 103
 104         CTYPE("Use mass-weighting of the rotation group positions");
 105         sprintf(buf, "rot_massw%d", g);
 106         ETYPE(buf, rotg->bMassW, yesno_names);
 107
 108         CTYPE("Rotation vector, will get normalized");
 109         sprintf(buf, "rot_vec%d", g);
 110         STYPE(buf, s_vec, "1.0 0.0 0.0");
 111         string2dvec(s_vec, vec);
 112         /* Normalize the rotation vector */
 113         if (dnorm(vec) != 0)
 114         {
 115             dsvmul(1.0/dnorm(vec), vec, vec);
 116         }
 117         else
 118         {
 119             sprintf(warn_buf, "rot_vec%d = 0", g);
 120             warning_error(wi, warn_buf);
 121         }
 122         fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
 123                 RotStr, g, erotg_names[rotg->eType], vec[0], vec[1], vec[2]);
 124         for (m = 0; m < DIM; m++)
 125         {
 126             rotg->vec[m] = vec[m];
 127         }
 128
 129         CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
 130         sprintf(buf, "rot_pivot%d", g);
 131         STYPE(buf, s_vec, "0.0 0.0 0.0");
 132         clear_dvec(vec);
 133         if ( (rotg->eType == erotgISO) || (rotg->eType == erotgPM) || (rotg->eType == erotgRM) || (rotg->eType == erotgRM2) )
 134         {
 135             string2dvec(s_vec, vec);
 136         }
 137         for (m = 0; m < DIM; m++)
 138         {
 139             rotg->pivot[m] = vec[m];
 140         }
 141
 142         CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
 143         sprintf(buf, "rot_rate%d", g);
 144         RTYPE(buf, rotg->rate, 0.0);
 145
 146         sprintf(buf, "rot_k%d", g);
 147         RTYPE(buf, rotg->k, 0.0);
 148         if (rotg->k <= 0.0)
 149         {
 150             sprintf(warn_buf, "rot_k%d <= 0", g);
 151             warning_note(wi, warn_buf);
 152         }
 153
 154         CTYPE("Slab distance for flexible axis rotation (nm)");
 155         sprintf(buf, "rot_slab_dist%d", g);
 156         RTYPE(buf, rotg->slab_dist, 1.5);
 157         if (rotg->slab_dist <= 0.0)
 158         {
 159             sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
 160             warning_error(wi, warn_buf);
 161         }
 162
 163         CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
 164         sprintf(buf, "rot_min_gauss%d", g);
 165         RTYPE(buf, rotg->min_gaussian, 1e-3);
 166         if (rotg->min_gaussian <= 0.0)
 167         {
 168             sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
 169             warning_error(wi, warn_buf);
 170         }
 171
 172         CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
 173         sprintf(buf, "rot_eps%d", g);
 174         RTYPE(buf, rotg->eps, 1e-4);
 175         if ( (rotg->eps <= 0.0) && (rotg->eType == erotgRM2 || rotg->eType == erotgFLEX2) )
 176         {
 177             sprintf(warn_buf, "rot_eps%d <= 0", g);
 178             warning_error(wi, warn_buf);
 179         }
 180
 181         CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
 182         sprintf(buf, "rot_fit_method%d", g);
 183         ETYPE(buf, rotg->eFittype, erotg_fitnames);
 184         CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
 185         sprintf(buf, "rot_potfit_nsteps%d", g);
 186         ITYPE(buf, rotg->PotAngle_nstep, 21);
 187         if ( (rotg->eFittype == erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
 188         {
 189             sprintf(warn_buf, "rot_potfit_nsteps%d < 1", g);
 190             warning_error(wi, warn_buf);
 191         }
 192         CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
 193         sprintf(buf, "rot_potfit_step%d", g);
 194         RTYPE(buf, rotg->PotAngle_step, 0.25);
 195     }
 196
 197     *ninp_p   = ninp;
 198     *inp_p    = inp;
 199
 200     return grpbuf;
 201 }
 202
 203
 204 /* Check whether the box is unchanged */
 205 static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
 206 {
 207     int      i, ii;
 208     gmx_bool bSame = TRUE;
 209     char     warn_buf[STRLEN];
 210
 211
 212     for (i = 0; i < DIM; i++)
 213     {
 214         for (ii = 0; ii < DIM; ii++)
 215         {
 216             if (f_box[i][ii] != box[i][ii])
 217             {
 218                 bSame = FALSE;
 219             }
 220         }
 221     }
 222     if (!bSame)
 223     {
 224         sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
 225                 RotStr, fn);
 226         warning(wi, warn_buf);
 227         pr_rvecs(stderr, 0, "Your box is:", box, 3);
 228         pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
 229     }
 230 }
 231
 232
 233 /* Extract the reference positions for the rotation group(s) */
 234 extern void set_reference_positions(
 235         t_rot *rot, gmx_mtop_t *mtop, rvec *x, matrix box,
 236         const char *fn, gmx_bool bSet, warninp_t wi)
 237 {
 238     int         g, i, ii;
 239     t_rotgrp   *rotg;
 240     t_trnheader header;    /* Header information of reference file */
 241     char        base[STRLEN], extension[STRLEN], reffile[STRLEN];
 242     char       *extpos;
 243     rvec        f_box[3];  /* Box from reference file */
 244
 245
 246     /* Base name and extension of the reference file: */
 247     strncpy(base, fn, STRLEN - 1);
 248     base[STRLEN-1] = '\0';
 249     extpos         = strrchr(base, '.');
 250     strcpy(extension, extpos+1);
 251     *extpos = '\0';
 252
 253
 254     for (g = 0; g < rot->ngrp; g++)
 255     {
 256         rotg = &rot->grp[g];
 257         fprintf(stderr, "%s group %d has %d reference positions.\n", RotStr, g, rotg->nat);
 258         snew(rotg->x_ref, rotg->nat);
 259
 260         /* Construct the name for the file containing the reference positions for this group: */
 261         sprintf(reffile, "%s.%d.%s", base, g, extension);
 262
 263         /* If the base filename for the reference position files was explicitly set by
 264          * the user, we issue a fatal error if the group file can not be found */
 265         if (bSet && !gmx_fexist(reffile))
 266         {
 267             gmx_fatal(FARGS, "%s The file containing the reference positions was not found.\n"
 268                       "Expected the file '%s' for group %d.\n",
 269                       RotStr, reffile, g);
 270         }
 271
 272         if (gmx_fexist(reffile))
 273         {
 274             fprintf(stderr, "  Reading them from %s.\n", reffile);
 275             read_trnheader(reffile, &header);
 276             if (rotg->nat != header.natoms)
 277             {
 278                 gmx_fatal(FARGS, "Number of atoms in file %s (%d) does not match the number of atoms in rotation group (%d)!\n",
 279                           reffile, header.natoms, rotg->nat);
 280             }
 281             read_trn(reffile, &header.step, &header.t, &header.lambda, f_box, &header.natoms, rotg->x_ref, NULL, NULL);
 282
 283             /* Check whether the box is unchanged and output a warning if not: */
 284             check_box_unchanged(f_box, box, reffile, wi);
 285         }
 286         else
 287         {
 288             fprintf(stderr, " Saving them to %s.\n", reffile);
 289             for (i = 0; i < rotg->nat; i++)
 290             {
 291                 ii = rotg->ind[i];
 292                 copy_rvec(x[ii], rotg->x_ref[i]);
 293             }
 294             write_trn(reffile, g, 0.0, 0.0, box, rotg->nat, rotg->x_ref, NULL, NULL);
 295         }
 296     }
 297 }
 298
 299
 300 extern void make_rotation_groups(t_rot *rot, char **rotgnames, t_blocka *grps, char **gnames)
 301 {
 302     int       g, ig = -1, i;
 303     t_rotgrp *rotg;
 304
 305
 306     for (g = 0; g < rot->ngrp; g++)
 307     {
 308         rotg      = &rot->grp[g];
 309         ig        = search_string(rotgnames[g], grps->nr, gnames);
 310         rotg->nat = grps->index[ig+1] - grps->index[ig];
 311
 312         if (rotg->nat > 0)
 313         {
 314             fprintf(stderr, "Rotation group %d '%s' has %d atoms\n", g, rotgnames[g], rotg->nat);
 315             snew(rotg->ind, rotg->nat);
 316             for (i = 0; i < rotg->nat; i++)
 317             {
 318                 rotg->ind[i] = grps->a[grps->index[ig]+i];
 319             }
 320         }
 321         else
 322         {
 323             gmx_fatal(FARGS, "Rotation group %d '%s' is empty", g, rotgnames[g]);
 324         }
 325     }
 326 }