src/gromacs/gmxpreprocess/readpull.c

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