src/gromacs/gmxpreprocess/readpull.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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <string.h>
  42 #include <stdlib.h>
  43
  44 #include "gromacs/utility/cstringutil.h"
  45 #include "gromacs/utility/futil.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "typedefs.h"
  49 #include "names.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "macros.h"
  52 #include "index.h"
  53 #include "readinp.h"
  54 #include "readir.h"
  55 #include "mdatoms.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/pulling/pull.h"
  58
  59
  60 static char pulldim[STRLEN];
  61
  62 static void string2dvec(const char buf[], dvec nums)
  63 {
  64     double dum;
  65
  66     if (sscanf(buf, "%lf%lf%lf%lf", &nums[0], &nums[1], &nums[2], &dum) != 3)
  67     {
  68         gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
  69     }
  70 }
  71
  72 static void init_pull_group(t_pull_group *pg,
  73                             const char   *wbuf)
  74 {
  75     double d;
  76     int    n, m;
  77
  78     pg->nweight = 0;
  79     while (sscanf(wbuf, "%lf %n", &d, &n) == 1)
  80     {
  81         if (pg->nweight % 100 == 0)
  82         {
  83             srenew(pg->weight, pg->nweight+100);
  84         }
  85         pg->weight[pg->nweight++] = d;
  86         wbuf += n;
  87     }
  88 }
  89
  90 static void init_pull_coord(t_pull_coord *pcrd, int eGeom,
  91                             const char *origin_buf, const char *vec_buf)
  92 {
  93     int    m;
  94     dvec   origin, vec;
  95
  96     string2dvec(origin_buf, origin);
  97     if (pcrd->group[0] != 0 && dnorm(origin) > 0)
  98     {
  99         gmx_fatal(FARGS, "The pull origin can only be set with an absolute reference");
 100     }
 101
 102     if (eGeom == epullgDIST)
 103     {
 104         clear_dvec(vec);
 105     }
 106     else
 107     {
 108         string2dvec(vec_buf, vec);
 109         if (eGeom == epullgDIR || eGeom == epullgCYL)
 110         {
 111             /* Normalize the direction vector */
 112             dsvmul(1/dnorm(vec), vec, vec);
 113         }
 114     }
 115     for (m = 0; m < DIM; m++)
 116     {
 117         pcrd->origin[m] = origin[m];
 118         pcrd->vec[m]    = vec[m];
 119     }
 120 }
 121
 122 char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
 123                        t_pull *pull, gmx_bool *bStart,
 124                        warninp_t wi)
 125 {
 126     int           ninp, nerror = 0, i, nchar, nscan, m, idum;
 127     t_inpfile    *inp;
 128     const char   *tmp;
 129     char        **grpbuf;
 130     char          dummy[STRLEN], buf[STRLEN], groups[STRLEN], init[STRLEN];
 131     const char   *init_def1 = "0.0", *init_def3 = "0.0 0.0 0.0";
 132     char          wbuf[STRLEN], origin_buf[STRLEN], vec_buf[STRLEN];
 133
 134     t_pull_group *pgrp;
 135     t_pull_coord *pcrd;
 136
 137     ninp   = *ninp_p;
 138     inp    = *inp_p;
 139
 140     /* read pull parameters */
 141     CTYPE("Pull geometry: distance, direction, direction-periodic or cylinder");
 142     EETYPE("pull-geometry",   pull->eGeom, epullg_names);
 143     CTYPE("Select components for the pull vector. default: Y Y Y");
 144     STYPE("pull-dim",         pulldim, "Y Y Y");
 145     CTYPE("Cylinder radius for dynamic reaction force groups (nm)");
 146     RTYPE("pull-r1",          pull->cyl_r1, 1.0);
 147     CTYPE("Switch from r1 to r0 in case of dynamic reaction force");
 148     RTYPE("pull-r0",          pull->cyl_r0, 1.5);
 149     RTYPE("pull_constr_tol",  pull->constr_tol, 1E-6);
 150     EETYPE("pull-start",      *bStart, yesno_names);
 151     EETYPE("pull-print-reference", pull->bPrintRef, yesno_names);
 152     ITYPE("pull-nstxout",     pull->nstxout, 10);
 153     ITYPE("pull-nstfout",     pull->nstfout,  1);
 154     CTYPE("Number of pull groups");
 155     ITYPE("pull-ngroups",     pull->ngroup, 1);
 156     CTYPE("Number of pull coordinates");
 157     ITYPE("pull-ncoords",     pull->ncoord, 1);
 158
 159     if (pull->cyl_r1 > pull->cyl_r0)
 160     {
 161         warning_error(wi, "pull-r1 > pull_r0");
 162     }
 163
 164     if (pull->ngroup < 1)
 165     {
 166         gmx_fatal(FARGS, "pull-ngroups should be >= 1");
 167     }
 168     /* We always add an absolute reference group (index 0), even if not used */
 169     pull->ngroup += 1;
 170
 171     if (pull->ncoord < 1)
 172     {
 173         gmx_fatal(FARGS, "pull-ncoords should be >= 1");
 174     }
 175
 176     snew(pull->group, pull->ngroup);
 177
 178     snew(pull->coord, pull->ncoord);
 179
 180     /* pull group options */
 181     CTYPE("Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2)");
 182
 183     /* Read the pull groups */
 184     snew(grpbuf, pull->ngroup);
 185     /* Group 0 is the absolute reference, we don't read anything for 0 */
 186     for (i = 1; i < pull->ngroup; i++)
 187     {
 188         pgrp = &pull->group[i];
 189         snew(grpbuf[i], STRLEN);
 190         sprintf(buf, "pull-group%d-name", i);
 191         STYPE(buf,              grpbuf[i], "");
 192         sprintf(buf, "pull-group%d-weights", i);
 193         STYPE(buf,              wbuf, "");
 194         sprintf(buf, "pull-group%d-pbcatom", i);
 195         ITYPE(buf,              pgrp->pbcatom, 0);
 196
 197         /* Initialize the pull group */
 198         init_pull_group(pgrp, wbuf);
 199     }
 200
 201     /* Read the pull coordinates */
 202     for (i = 1; i < pull->ncoord + 1; i++)
 203     {
 204         pcrd = &pull->coord[i-1];
 205         sprintf(buf, "pull-coord%d-groups", i);
 206         STYPE(buf,              groups, "");
 207         nscan = sscanf(groups, "%d %d %d", &pcrd->group[0], &pcrd->group[1], &idum);
 208         if (nscan != 2)
 209         {
 210             fprintf(stderr, "ERROR: %s should have %d components\n", buf, 2);
 211             nerror++;
 212         }
 213         sprintf(buf, "pull-coord%d-origin", i);
 214         STYPE(buf,              origin_buf, "0.0 0.0 0.0");
 215         sprintf(buf, "pull-coord%d-vec", i);
 216         STYPE(buf,              vec_buf,    "0.0 0.0 0.0");
 217         sprintf(buf, "pull-coord%d-init", i);
 218         RTYPE(buf,              pcrd->init, 0.0);
 219         sprintf(buf, "pull-coord%d-rate", i);
 220         RTYPE(buf,              pcrd->rate, 0.0);
 221         sprintf(buf, "pull-coord%d-k", i);
 222         RTYPE(buf,              pcrd->k, 0.0);
 223         sprintf(buf, "pull-coord%d-kB", i);
 224         RTYPE(buf,              pcrd->kB, pcrd->k);
 225
 226         /* Initialize the pull coordinate */
 227         init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf);
 228     }
 229
 230     *ninp_p   = ninp;
 231     *inp_p    = inp;
 232
 233     return grpbuf;
 234 }
 235
 236 void make_pull_groups(t_pull *pull,
 237                       char **pgnames,
 238                       const t_blocka *grps, char **gnames)
 239 {
 240     int           g, ig = -1, i;
 241     t_pull_group *pgrp;
 242
 243     /* Absolute reference group (might not be used) is special */
 244     pgrp          = &pull->group[0];
 245     pgrp->nat     = 0;
 246     pgrp->pbcatom = -1;
 247
 248     for (g = 1; g < pull->ngroup; g++)
 249     {
 250         pgrp = &pull->group[g];
 251
 252         ig        = search_string(pgnames[g], grps->nr, gnames);
 253         pgrp->nat = grps->index[ig+1] - grps->index[ig];
 254
 255         fprintf(stderr, "Pull group %d '%s' has %d atoms\n",
 256                 g, pgnames[g], pgrp->nat);
 257
 258         if (pgrp->nat == 0)
 259         {
 260             gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pgnames[g]);
 261         }
 262
 263         snew(pgrp->ind, pgrp->nat);
 264         for (i = 0; i < pgrp->nat; i++)
 265         {
 266             pgrp->ind[i] = grps->a[grps->index[ig]+i];
 267         }
 268
 269         if (pull->eGeom == epullgCYL && g == 1 && pgrp->nweight > 0)
 270         {
 271             gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling");
 272         }
 273         if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat)
 274         {
 275             gmx_fatal(FARGS, "Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)",
 276                       pgrp->nweight, g, pgnames[g], pgrp->nat);
 277         }
 278
 279         if (pgrp->nat == 1)
 280         {
 281             /* No pbc is required for this group */
 282             pgrp->pbcatom = -1;
 283         }
 284         else
 285         {
 286             if (pgrp->pbcatom > 0)
 287             {
 288                 pgrp->pbcatom -= 1;
 289             }
 290             else if (pgrp->pbcatom == 0)
 291             {
 292                 pgrp->pbcatom = pgrp->ind[(pgrp->nat-1)/2];
 293             }
 294             else
 295             {
 296                 /* Use cosine weighting */
 297                 pgrp->pbcatom = -1;
 298             }
 299         }
 300     }
 301 }
 302
 303 void make_pull_coords(t_pull *pull)
 304 {
 305     int           ndim, d, nchar, c;
 306     char         *ptr, pulldim1[STRLEN];
 307     t_pull_coord *pcrd;
 308
 309     ptr  = pulldim;
 310     ndim = 0;
 311     for (d = 0; d < DIM; d++)
 312     {
 313         if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1)
 314         {
 315             gmx_fatal(FARGS, "Less than 3 pull dimensions given in pull_dim: '%s'",
 316                       pulldim);
 317         }
 318
 319         if (gmx_strncasecmp(pulldim1, "N", 1) == 0)
 320         {
 321             pull->dim[d] = 0;
 322         }
 323         else if (gmx_strncasecmp(pulldim1, "Y", 1) == 0)
 324         {
 325             pull->dim[d] = 1;
 326             ndim++;
 327         }
 328         else
 329         {
 330             gmx_fatal(FARGS, "Please use Y(ES) or N(O) for pull_dim only (not %s)",
 331                       pulldim1);
 332         }
 333         ptr += nchar;
 334     }
 335     if (ndim == 0)
 336     {
 337         gmx_fatal(FARGS, "All entries in pull_dim are N");
 338     }
 339
 340     for (c = 0; c < pull->ncoord; c++)
 341     {
 342         pcrd = &pull->coord[c];
 343
 344         if (pcrd->group[0] < 0 || pcrd->group[0] >= pull->ngroup ||
 345             pcrd->group[1] < 0 || pcrd->group[1] >= pull->ngroup)
 346         {
 347             gmx_fatal(FARGS, "Pull group index in pull-coord%d-groups out of range, should be between %d and %d", c+1, 0, pull->ngroup+1);
 348         }
 349
 350         if (pcrd->group[0] == pcrd->group[1])
 351         {
 352             gmx_fatal(FARGS, "Identical pull group indices in pull-coord%d-groups", c+1);
 353         }
 354
 355         if (pull->eGeom == epullgCYL && pcrd->group[0] != 1)
 356         {
 357             gmx_fatal(FARGS, "With pull geometry %s, the first pull group should always be 1", EPULLGEOM(pull->eGeom));
 358         }
 359
 360         if (pull->eGeom != epullgDIST)
 361         {
 362             for (d = 0; d < DIM; d++)
 363             {
 364                 if (pcrd->vec[d] != 0 && pull->dim[d] == 0)
 365                 {
 366                     gmx_fatal(FARGS, "ERROR: pull-group%d-vec has non-zero %c-component while pull_dim for the %c-dimension is N\n", c+1, 'x'+d, 'x'+d);
 367                 }
 368             }
 369         }
 370
 371         if ((pull->eGeom == epullgDIR || pull->eGeom == epullgCYL) &&
 372             norm2(pcrd->vec) == 0)
 373         {
 374             gmx_fatal(FARGS, "pull-group%d-vec can not be zero with geometry %s",
 375                       c+1, EPULLGEOM(pull->eGeom));
 376         }
 377     }
 378 }
 379
 380 void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
 381                    const output_env_t oenv, gmx_bool bStart)
 382 {
 383     t_mdatoms    *md;
 384     t_pull       *pull;
 385     t_pull_coord *pcrd;
 386     t_pull_group *pgrp0, *pgrp1;
 387     t_pbc         pbc;
 388     int           c, m;
 389     double        t_start, tinvrate;
 390     real          init;
 391     dvec          dr;
 392     double        dev;
 393
 394     init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
 395     md = init_mdatoms(NULL, mtop, ir->efep);
 396     atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
 397     if (ir->efep)
 398     {
 399         update_mdatoms(md, lambda);
 400     }
 401     pull = ir->pull;
 402
 403     set_pbc(&pbc, ir->ePBC, box);
 404
 405     t_start = ir->init_t + ir->init_step*ir->delta_t;
 406
 407     pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL);
 408
 409     fprintf(stderr, "Pull group  natoms  pbc atom  distance at start     reference at t=0\n");
 410     for (c = 0; c < pull->ncoord; c++)
 411     {
 412         pcrd  = &pull->coord[c];
 413
 414         pgrp0 = &pull->group[pcrd->group[0]];
 415         pgrp1 = &pull->group[pcrd->group[1]];
 416         fprintf(stderr, "%8d  %8d  %8d\n",
 417                 pcrd->group[0], pgrp0->nat, pgrp0->pbcatom+1);
 418         fprintf(stderr, "%8d  %8d  %8d ",
 419                 pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);
 420
 421         init       = pcrd->init;
 422         pcrd->init = 0;
 423
 424         if (pcrd->rate == 0)
 425         {
 426             tinvrate = 0;
 427         }
 428         else
 429         {
 430             tinvrate = t_start/pcrd->rate;
 431         }
 432         get_pull_coord_distance(pull, c, &pbc, 0, dr, &dev);
 433         fprintf(stderr, " %6.3f ", dev);
 434
 435         if (bStart)
 436         {
 437             pcrd->init = init + dev - tinvrate;
 438         }
 439         else
 440         {
 441             pcrd->init = init;
 442         }
 443         fprintf(stderr, " %6.3f\n", pcrd->init);
 444     }
 445 }