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