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