src/mdlib/pull.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
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include "futil.h"
  47 #include "index.h"
  48 #include "statutil.h"
  49 #include "gmxfio.h"
  50 #include "vec.h"
  51 #include "typedefs.h"
  52 #include "network.h"
  53 #include "filenm.h"
  54 #include <string.h>
  55 #include "smalloc.h"
  56 #include "pull.h"
  57 #include "xvgr.h"
  58 #include "names.h"
  59 #include "partdec.h"
  60 #include "pbc.h"
  61 #include "mtop_util.h"
  62 #include "mdrun.h"
  63 #include "gmx_ga2la.h"
  64 #include "copyrite.h"
  65
  66 static void pull_print_x_grp(FILE *out, gmx_bool bRef, ivec dim, t_pullgrp *pgrp)
  67 {
  68     int m;
  69
  70     for (m = 0; m < DIM; m++)
  71     {
  72         if (dim[m])
  73         {
  74             fprintf(out, "\t%g", bRef ? pgrp->x[m] : pgrp->dr[m]);
  75         }
  76     }
  77 }
  78
  79 static void pull_print_x(FILE *out, t_pull *pull, double t)
  80 {
  81     int g;
  82
  83     fprintf(out, "%.4f", t);
  84
  85     if (PULL_CYL(pull))
  86     {
  87         for (g = 1; g < 1+pull->ngrp; g++)
  88         {
  89             pull_print_x_grp(out, TRUE, pull->dim, &pull->dyna[g]);
  90             pull_print_x_grp(out, FALSE, pull->dim, &pull->grp[g]);
  91         }
  92     }
  93     else
  94     {
  95         for (g = 0; g < 1+pull->ngrp; g++)
  96         {
  97             if (pull->grp[g].nat > 0)
  98             {
  99                 pull_print_x_grp(out, g == 0, pull->dim, &pull->grp[g]);
 100             }
 101         }
 102     }
 103     fprintf(out, "\n");
 104 }
 105
 106 static void pull_print_f(FILE *out, t_pull *pull, double t)
 107 {
 108     int g, d;
 109
 110     fprintf(out, "%.4f", t);
 111
 112     for (g = 1; g < 1+pull->ngrp; g++)
 113     {
 114         if (pull->eGeom == epullgPOS)
 115         {
 116             for (d = 0; d < DIM; d++)
 117             {
 118                 if (pull->dim[d])
 119                 {
 120                     fprintf(out, "\t%g", pull->grp[g].f[d]);
 121                 }
 122             }
 123         }
 124         else
 125         {
 126             fprintf(out, "\t%g", pull->grp[g].f_scal);
 127         }
 128     }
 129     fprintf(out, "\n");
 130 }
 131
 132 void pull_print_output(t_pull *pull, gmx_large_int_t step, double time)
 133 {
 134     if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
 135     {
 136         pull_print_x(pull->out_x, pull, time);
 137     }
 138
 139     if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
 140     {
 141         pull_print_f(pull->out_f, pull, time);
 142     }
 143 }
 144
 145 static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
 146                            gmx_bool bCoord, unsigned long Flags)
 147 {
 148     FILE  *fp;
 149     int    nsets, g, m;
 150     char **setname, buf[10];
 151
 152     if (Flags & MD_APPENDFILES)
 153     {
 154         fp = gmx_fio_fopen(fn, "a+");
 155     }
 156     else
 157     {
 158         fp = gmx_fio_fopen(fn, "w+");
 159         if (bCoord)
 160         {
 161             xvgr_header(fp, "Pull COM",  "Time (ps)", "Position (nm)",
 162                         exvggtXNY, oenv);
 163         }
 164         else
 165         {
 166             xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
 167                         exvggtXNY, oenv);
 168         }
 169
 170         snew(setname, (1+pull->ngrp)*DIM);
 171         nsets = 0;
 172         for (g = 0; g < 1+pull->ngrp; g++)
 173         {
 174             if (pull->grp[g].nat > 0 &&
 175                 (g > 0 || (bCoord && !PULL_CYL(pull))))
 176             {
 177                 if (bCoord || pull->eGeom == epullgPOS)
 178                 {
 179                     if (PULL_CYL(pull))
 180                     {
 181                         for (m = 0; m < DIM; m++)
 182                         {
 183                             if (pull->dim[m])
 184                             {
 185                                 sprintf(buf, "%d %s%c", g, "c", 'X'+m);
 186                                 setname[nsets] = strdup(buf);
 187                                 nsets++;
 188                             }
 189                         }
 190                     }
 191                     for (m = 0; m < DIM; m++)
 192                     {
 193                         if (pull->dim[m])
 194                         {
 195                             sprintf(buf, "%d %s%c",
 196                                     g, (bCoord && g > 0) ? "d" : "", 'X'+m);
 197                             setname[nsets] = strdup(buf);
 198                             nsets++;
 199                         }
 200                     }
 201                 }
 202                 else
 203                 {
 204                     sprintf(buf, "%d", g);
 205                     setname[nsets] = strdup(buf);
 206                     nsets++;
 207                 }
 208             }
 209         }
 210         if (bCoord || nsets > 1)
 211         {
 212             xvgr_legend(fp, nsets, (const char**)setname, oenv);
 213         }
 214         for (g = 0; g < nsets; g++)
 215         {
 216             sfree(setname[g]);
 217         }
 218         sfree(setname);
 219     }
 220
 221     return fp;
 222 }
 223
 224 /* Apply forces in a mass weighted fashion */
 225 static void apply_forces_grp(t_pullgrp *pgrp, t_mdatoms * md,
 226                              gmx_ga2la_t ga2la,
 227                              dvec f_pull, int sign, rvec *f)
 228 {
 229     int    i, ii, m, start, end;
 230     double wmass, inv_wm;
 231
 232     start = md->start;
 233     end   = md->homenr + start;
 234
 235     inv_wm = pgrp->wscale*pgrp->invtm;
 236
 237     for (i = 0; i < pgrp->nat_loc; i++)
 238     {
 239         ii    = pgrp->ind_loc[i];
 240         wmass = md->massT[ii];
 241         if (pgrp->weight_loc)
 242         {
 243             wmass *= pgrp->weight_loc[i];
 244         }
 245
 246         for (m = 0; m < DIM; m++)
 247         {
 248             f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
 249         }
 250     }
 251 }
 252
 253 /* Apply forces in a mass weighted fashion */
 254 static void apply_forces(t_pull * pull, t_mdatoms * md, gmx_ga2la_t ga2la,
 255                          rvec *f)
 256 {
 257     int        i;
 258     t_pullgrp *pgrp;
 259
 260     for (i = 1; i < pull->ngrp+1; i++)
 261     {
 262         pgrp = &(pull->grp[i]);
 263         apply_forces_grp(pgrp, md, ga2la, pgrp->f, 1, f);
 264         if (pull->grp[0].nat)
 265         {
 266             if (PULL_CYL(pull))
 267             {
 268                 apply_forces_grp(&(pull->dyna[i]), md, ga2la, pgrp->f, -1, f);
 269             }
 270             else
 271             {
 272                 apply_forces_grp(&(pull->grp[0]), md, ga2la, pgrp->f, -1, f);
 273             }
 274         }
 275     }
 276 }
 277
 278 static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc)
 279 {
 280     double max_d2;
 281     int    m;
 282
 283     max_d2 = GMX_DOUBLE_MAX;
 284
 285     if (pull->eGeom != epullgDIRPBC)
 286     {
 287         for (m = 0; m < pbc->ndim_ePBC; m++)
 288         {
 289             if (pull->dim[m] != 0)
 290             {
 291                 max_d2 = min(max_d2, norm2(pbc->box[m]));
 292             }
 293         }
 294     }
 295
 296     return 0.25*max_d2;
 297 }
 298
 299 static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double t,
 300                             dvec xg, dvec xref, double max_dist2,
 301                             dvec dr)
 302 {
 303     t_pullgrp *pref, *pgrp;
 304     int        m;
 305     dvec       xrefr, dref = {0, 0, 0};
 306     double     dr2;
 307
 308     pgrp = &pull->grp[g];
 309
 310     copy_dvec(xref, xrefr);
 311
 312     if (pull->eGeom == epullgDIRPBC)
 313     {
 314         for (m = 0; m < DIM; m++)
 315         {
 316             dref[m] = (pgrp->init[0] + pgrp->rate*t)*pull->grp[g].vec[m];
 317         }
 318         /* Add the reference position, so we use the correct periodic image */
 319         dvec_inc(xrefr, dref);
 320     }
 321
 322     pbc_dx_d(pbc, xg, xrefr, dr);
 323     dr2 = 0;
 324     for (m = 0; m < DIM; m++)
 325     {
 326         dr[m] *= pull->dim[m];
 327         dr2   += dr[m]*dr[m];
 328     }
 329     if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
 330     {
 331         gmx_fatal(FARGS, "Distance of pull group %d (%f nm) is larger than 0.49 times the box size (%f).\nYou might want to consider using \"pull-geometry = direction-periodic\" instead.\n", g, sqrt(dr2), sqrt(max_dist2));
 332
 333     }
 334
 335     if (pull->eGeom == epullgDIRPBC)
 336     {
 337         dvec_inc(dr, dref);
 338     }
 339 }
 340
 341 static void get_pullgrp_dr(const t_pull *pull, const t_pbc *pbc, int g, double t,
 342                            dvec dr)
 343 {
 344     double md2;
 345
 346     if (pull->eGeom == epullgDIRPBC)
 347     {
 348         md2 = -1;
 349     }
 350     else
 351     {
 352         md2 = max_pull_distance2(pull, pbc);
 353     }
 354
 355     get_pullgrps_dr(pull, pbc, g, t,
 356                     pull->grp[g].x,
 357                     PULL_CYL(pull) ? pull->dyna[g].x : pull->grp[0].x,
 358                     md2,
 359                     dr);
 360 }
 361
 362 void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t,
 363                           dvec dr, dvec dev)
 364 {
 365     static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
 366                                         but is fairly benign */
 367     t_pullgrp      *pgrp;
 368     int             m;
 369     dvec            ref;
 370     double          drs, inpr;
 371
 372     pgrp = &pull->grp[g];
 373
 374     get_pullgrp_dr(pull, pbc, g, t, dr);
 375
 376     if (pull->eGeom == epullgPOS)
 377     {
 378         for (m = 0; m < DIM; m++)
 379         {
 380             ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
 381         }
 382     }
 383     else
 384     {
 385         ref[0] = pgrp->init[0] + pgrp->rate*t;
 386     }
 387
 388     switch (pull->eGeom)
 389     {
 390         case epullgDIST:
 391             /* Pull along the vector between the com's */
 392             if (ref[0] < 0 && !bWarned)
 393             {
 394                 fprintf(stderr, "\nPull reference distance for group %d is negative (%f)\n", g, ref[0]);
 395                 bWarned = TRUE;
 396             }
 397             drs = dnorm(dr);
 398             if (drs == 0)
 399             {
 400                 /* With no vector we can not determine the direction for the force,
 401                  * so we set the force to zero.
 402                  */
 403                 dev[0] = 0;
 404             }
 405             else
 406             {
 407                 /* Determine the deviation */
 408                 dev[0] = drs - ref[0];
 409             }
 410             break;
 411         case epullgDIR:
 412         case epullgDIRPBC:
 413         case epullgCYL:
 414             /* Pull along vec */
 415             inpr = 0;
 416             for (m = 0; m < DIM; m++)
 417             {
 418                 inpr += pgrp->vec[m]*dr[m];
 419             }
 420             dev[0] = inpr - ref[0];
 421             break;
 422         case epullgPOS:
 423             /* Determine the difference of dr and ref along each dimension */
 424             for (m = 0; m < DIM; m++)
 425             {
 426                 dev[m] = (dr[m] - ref[m])*pull->dim[m];
 427             }
 428             break;
 429     }
 430 }
 431
 432 void clear_pull_forces(t_pull *pull)
 433 {
 434     int i;
 435
 436     /* Zeroing the forces is only required for constraint pulling.
 437      * It can happen that multiple constraint steps need to be applied
 438      * and therefore the constraint forces need to be accumulated.
 439      */
 440     for (i = 0; i < 1+pull->ngrp; i++)
 441     {
 442         clear_dvec(pull->grp[i].f);
 443         pull->grp[i].f_scal = 0;
 444     }
 445 }
 446
 447 /* Apply constraint using SHAKE */
 448 static void do_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
 449                           rvec *x, rvec *v,
 450                           gmx_bool bMaster, tensor vir,
 451                           double dt, double t)
 452 {
 453
 454     dvec      *r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
 455     dvec       unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
 456
 457     dvec      *rinew;  /* current 'new' position of group i */
 458     dvec      *rjnew;  /* current 'new' position of group j */
 459     dvec       ref, vec;
 460     double     d0, inpr;
 461     double     lambda, rm, mass, invdt = 0;
 462     gmx_bool   bConverged_all, bConverged = FALSE;
 463     int        niter = 0, g, ii, j, m, max_iter = 100;
 464     double     q, a, b, c; /* for solving the quadratic equation,
 465                               see Num. Recipes in C ed 2 p. 184 */
 466     dvec      *dr;         /* correction for group i */
 467     dvec       ref_dr;     /* correction for group j */
 468     dvec       f;          /* the pull force */
 469     dvec       tmp, tmp3;
 470     t_pullgrp *pdyna, *pgrp, *pref;
 471
 472     snew(r_ij, pull->ngrp+1);
 473     if (PULL_CYL(pull))
 474     {
 475         snew(rjnew, pull->ngrp+1);
 476     }
 477     else
 478     {
 479         snew(rjnew, 1);
 480     }
 481     snew(dr, pull->ngrp+1);
 482     snew(rinew, pull->ngrp+1);
 483
 484     /* copy the current unconstrained positions for use in iterations. We
 485        iterate until rinew[i] and rjnew[j] obey the constraints. Then
 486        rinew - pull.x_unc[i] is the correction dr to group i */
 487     for (g = 1; g < 1+pull->ngrp; g++)
 488     {
 489         copy_dvec(pull->grp[g].xp, rinew[g]);
 490     }
 491     if (PULL_CYL(pull))
 492     {
 493         for (g = 1; g < 1+pull->ngrp; g++)
 494         {
 495             copy_dvec(pull->dyna[g].xp, rjnew[g]);
 496         }
 497     }
 498     else
 499     {
 500         copy_dvec(pull->grp[0].xp, rjnew[0]);
 501     }
 502
 503     /* Determine the constraint directions from the old positions */
 504     for (g = 1; g < 1+pull->ngrp; g++)
 505     {
 506         get_pullgrp_dr(pull, pbc, g, t, r_ij[g]);
 507         /* Store the difference vector at time t for printing */
 508         copy_dvec(r_ij[g], pull->grp[g].dr);
 509         if (debug)
 510         {
 511             fprintf(debug, "Pull group %d dr %f %f %f\n",
 512                     g, r_ij[g][XX], r_ij[g][YY], r_ij[g][ZZ]);
 513         }
 514
 515         if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
 516         {
 517             /* Select the component along vec */
 518             a = 0;
 519             for (m = 0; m < DIM; m++)
 520             {
 521                 a += pull->grp[g].vec[m]*r_ij[g][m];
 522             }
 523             for (m = 0; m < DIM; m++)
 524             {
 525                 r_ij[g][m] = a*pull->grp[g].vec[m];
 526             }
 527         }
 528     }
 529
 530     bConverged_all = FALSE;
 531     while (!bConverged_all && niter < max_iter)
 532     {
 533         bConverged_all = TRUE;
 534
 535         /* loop over all constraints */
 536         for (g = 1; g < 1+pull->ngrp; g++)
 537         {
 538             pgrp = &pull->grp[g];
 539             if (PULL_CYL(pull))
 540             {
 541                 pref = &pull->dyna[g];
 542             }
 543             else
 544             {
 545                 pref = &pull->grp[0];
 546             }
 547
 548             /* Get the current difference vector */
 549             get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0],
 550                             -1, unc_ij);
 551
 552             if (pull->eGeom == epullgPOS)
 553             {
 554                 for (m = 0; m < DIM; m++)
 555                 {
 556                     ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
 557                 }
 558             }
 559             else
 560             {
 561                 ref[0] = pgrp->init[0] + pgrp->rate*t;
 562                 /* Keep the compiler happy */
 563                 ref[1] = 0;
 564                 ref[2] = 0;
 565             }
 566
 567             if (debug)
 568             {
 569                 fprintf(debug, "Pull group %d, iteration %d\n", g, niter);
 570             }
 571
 572             rm = 1.0/(pull->grp[g].invtm + pref->invtm);
 573
 574             switch (pull->eGeom)
 575             {
 576                 case epullgDIST:
 577                     if (ref[0] <= 0)
 578                     {
 579                         gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", g, ref[0]);
 580                     }
 581
 582                     a = diprod(r_ij[g], r_ij[g]);
 583                     b = diprod(unc_ij, r_ij[g])*2;
 584                     c = diprod(unc_ij, unc_ij) - dsqr(ref[0]);
 585
 586                     if (b < 0)
 587                     {
 588                         q      = -0.5*(b - sqrt(b*b - 4*a*c));
 589                         lambda = -q/a;
 590                     }
 591                     else
 592                     {
 593                         q      = -0.5*(b + sqrt(b*b - 4*a*c));
 594                         lambda = -c/q;
 595                     }
 596
 597                     if (debug)
 598                     {
 599                         fprintf(debug,
 600                                 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
 601                                 a, b, c, lambda);
 602                     }
 603
 604                     /* The position corrections dr due to the constraints */
 605                     dsvmul(-lambda*rm*pgrp->invtm, r_ij[g],  dr[g]);
 606                     dsvmul( lambda*rm*pref->invtm, r_ij[g], ref_dr);
 607                     break;
 608                 case epullgDIR:
 609                 case epullgDIRPBC:
 610                 case epullgCYL:
 611                     /* A 1-dimensional constraint along a vector */
 612                     a = 0;
 613                     for (m = 0; m < DIM; m++)
 614                     {
 615                         vec[m] = pgrp->vec[m];
 616                         a     += unc_ij[m]*vec[m];
 617                     }
 618                     /* Select only the component along the vector */
 619                     dsvmul(a, vec, unc_ij);
 620                     lambda = a - ref[0];
 621                     if (debug)
 622                     {
 623                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
 624                     }
 625
 626                     /* The position corrections dr due to the constraints */
 627                     dsvmul(-lambda*rm*pull->grp[g].invtm, vec, dr[g]);
 628                     dsvmul( lambda*rm*       pref->invtm, vec, ref_dr);
 629                     break;
 630                 case epullgPOS:
 631                     for (m = 0; m < DIM; m++)
 632                     {
 633                         if (pull->dim[m])
 634                         {
 635                             lambda = r_ij[g][m] - ref[m];
 636                             /* The position corrections dr due to the constraints */
 637                             dr[g][m]  = -lambda*rm*pull->grp[g].invtm;
 638                             ref_dr[m] =  lambda*rm*pref->invtm;
 639                         }
 640                         else
 641                         {
 642                             dr[g][m]  = 0;
 643                             ref_dr[m] = 0;
 644                         }
 645                     }
 646                     break;
 647             }
 648
 649             /* DEBUG */
 650             if (debug)
 651             {
 652                 j = (PULL_CYL(pull) ? g : 0);
 653                 get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[j], -1, tmp);
 654                 get_pullgrps_dr(pull, pbc, g, t, dr[g], ref_dr, -1, tmp3);
 655                 fprintf(debug,
 656                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 657                         rinew[g][0], rinew[g][1], rinew[g][2],
 658                         rjnew[j][0], rjnew[j][1], rjnew[j][2], dnorm(tmp));
 659                 if (pull->eGeom == epullgPOS)
 660                 {
 661                     fprintf(debug,
 662                             "Pull ref %8.5f %8.5f %8.5f\n",
 663                             pgrp->vec[0], pgrp->vec[1], pgrp->vec[2]);
 664                 }
 665                 else
 666                 {
 667                     fprintf(debug,
 668                             "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f %8.5f %8.5f\n",
 669                             "", "", "", "", "", "", ref[0], ref[1], ref[2]);
 670                 }
 671                 fprintf(debug,
 672                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 673                         dr[g][0], dr[g][1], dr[g][2],
 674                         ref_dr[0], ref_dr[1], ref_dr[2],
 675                         dnorm(tmp3));
 676                 fprintf(debug,
 677                         "Pull cor %10.7f %10.7f %10.7f\n",
 678                         dr[g][0], dr[g][1], dr[g][2]);
 679             } /* END DEBUG */
 680
 681             /* Update the COMs with dr */
 682             dvec_inc(rinew[g],                     dr[g]);
 683             dvec_inc(rjnew[PULL_CYL(pull) ? g : 0], ref_dr);
 684         }
 685
 686         /* Check if all constraints are fullfilled now */
 687         for (g = 1; g < 1+pull->ngrp; g++)
 688         {
 689             pgrp = &pull->grp[g];
 690
 691             if (pull->eGeom == epullgPOS)
 692             {
 693                 for (m = 0; m < DIM; m++)
 694                 {
 695                     ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
 696                 }
 697             }
 698             else
 699             {
 700                 ref[0] = pgrp->init[0] + pgrp->rate*t;
 701                 /* Keep the compiler happy */
 702                 ref[1] = 0;
 703                 ref[2] = 0;
 704             }
 705
 706             get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0],
 707                             -1, unc_ij);
 708
 709             switch (pull->eGeom)
 710             {
 711                 case epullgDIST:
 712                     bConverged = fabs(dnorm(unc_ij) - ref[0]) < pull->constr_tol;
 713                     break;
 714                 case epullgDIR:
 715                 case epullgDIRPBC:
 716                 case epullgCYL:
 717                     for (m = 0; m < DIM; m++)
 718                     {
 719                         vec[m] = pgrp->vec[m];
 720                     }
 721                     inpr = diprod(unc_ij, vec);
 722                     dsvmul(inpr, vec, unc_ij);
 723                     bConverged =
 724                         fabs(diprod(unc_ij, vec) - ref[0]) < pull->constr_tol;
 725                     break;
 726                 case epullgPOS:
 727                     bConverged = TRUE;
 728                     for (m = 0; m < DIM; m++)
 729                     {
 730                         if (pull->dim[m] &&
 731                             fabs(unc_ij[m] - ref[m]) >= pull->constr_tol)
 732                         {
 733                             bConverged = FALSE;
 734                         }
 735                     }
 736                     break;
 737             }
 738
 739             if (!bConverged)
 740             {
 741                 if (debug)
 742                 {
 743                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
 744                             "d_ref = %f %f %f, current d = %f\n",
 745                             g, ref[0], ref[1], ref[2], dnorm(unc_ij));
 746                 }
 747
 748                 bConverged_all = FALSE;
 749             }
 750         }
 751
 752         niter++;
 753         /* if after all constraints are dealt with and bConverged is still TRUE
 754            we're finished, if not we do another iteration */
 755     }
 756     if (niter > max_iter)
 757     {
 758         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
 759     }
 760
 761     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 762
 763     if (v)
 764     {
 765         invdt = 1/dt;
 766     }
 767
 768     /* update the normal groups */
 769     for (g = 1; g < 1+pull->ngrp; g++)
 770     {
 771         pgrp = &pull->grp[g];
 772         /* get the final dr and constraint force for group i */
 773         dvec_sub(rinew[g], pgrp->xp, dr[g]);
 774         /* select components of dr */
 775         for (m = 0; m < DIM; m++)
 776         {
 777             dr[g][m] *= pull->dim[m];
 778         }
 779         dsvmul(1.0/(pgrp->invtm*dt*dt), dr[g], f);
 780         dvec_inc(pgrp->f, f);
 781         switch (pull->eGeom)
 782         {
 783             case epullgDIST:
 784                 for (m = 0; m < DIM; m++)
 785                 {
 786                     pgrp->f_scal += r_ij[g][m]*f[m]/dnorm(r_ij[g]);
 787                 }
 788                 break;
 789             case epullgDIR:
 790             case epullgDIRPBC:
 791             case epullgCYL:
 792                 for (m = 0; m < DIM; m++)
 793                 {
 794                     pgrp->f_scal += pgrp->vec[m]*f[m];
 795                 }
 796                 break;
 797             case epullgPOS:
 798                 break;
 799         }
 800
 801         if (vir && bMaster)
 802         {
 803             /* Add the pull contribution to the virial */
 804             for (j = 0; j < DIM; j++)
 805             {
 806                 for (m = 0; m < DIM; m++)
 807                 {
 808                     vir[j][m] -= 0.5*f[j]*r_ij[g][m];
 809                 }
 810             }
 811         }
 812
 813         /* update the atom positions */
 814         copy_dvec(dr[g], tmp);
 815         for (j = 0; j < pgrp->nat_loc; j++)
 816         {
 817             ii = pgrp->ind_loc[j];
 818             if (pgrp->weight_loc)
 819             {
 820                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr[g], tmp);
 821             }
 822             for (m = 0; m < DIM; m++)
 823             {
 824                 x[ii][m] += tmp[m];
 825             }
 826             if (v)
 827             {
 828                 for (m = 0; m < DIM; m++)
 829                 {
 830                     v[ii][m] += invdt*tmp[m];
 831                 }
 832             }
 833         }
 834     }
 835
 836     /* update the reference groups */
 837     if (PULL_CYL(pull))
 838     {
 839         /* update the dynamic reference groups */
 840         for (g = 1; g < 1+pull->ngrp; g++)
 841         {
 842             pdyna = &pull->dyna[g];
 843             dvec_sub(rjnew[g], pdyna->xp, ref_dr);
 844             /* select components of ref_dr */
 845             for (m = 0; m < DIM; m++)
 846             {
 847                 ref_dr[m] *= pull->dim[m];
 848             }
 849
 850             for (j = 0; j < pdyna->nat_loc; j++)
 851             {
 852                 /* reset the atoms with dr, weighted by w_i */
 853                 dsvmul(pdyna->wscale*pdyna->weight_loc[j], ref_dr, tmp);
 854                 ii = pdyna->ind_loc[j];
 855                 for (m = 0; m < DIM; m++)
 856                 {
 857                     x[ii][m] += tmp[m];
 858                 }
 859                 if (v)
 860                 {
 861                     for (m = 0; m < DIM; m++)
 862                     {
 863                         v[ii][m] += invdt*tmp[m];
 864                     }
 865                 }
 866             }
 867         }
 868     }
 869     else
 870     {
 871         pgrp = &pull->grp[0];
 872         /* update the reference group */
 873         dvec_sub(rjnew[0], pgrp->xp, ref_dr);
 874         /* select components of ref_dr */
 875         for (m = 0; m < DIM; m++)
 876         {
 877             ref_dr[m] *= pull->dim[m];
 878         }
 879
 880         copy_dvec(ref_dr, tmp);
 881         for (j = 0; j < pgrp->nat_loc; j++)
 882         {
 883             ii = pgrp->ind_loc[j];
 884             if (pgrp->weight_loc)
 885             {
 886                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], ref_dr, tmp);
 887             }
 888             for (m = 0; m < DIM; m++)
 889             {
 890                 x[ii][m] += tmp[m];
 891             }
 892             if (v)
 893             {
 894                 for (m = 0; m < DIM; m++)
 895                 {
 896                     v[ii][m] += invdt*tmp[m];
 897                 }
 898             }
 899         }
 900     }
 901
 902     /* finished! I hope. Give back some memory */
 903     sfree(r_ij);
 904     sfree(rinew);
 905     sfree(rjnew);
 906     sfree(dr);
 907 }
 908
 909 /* Pulling with a harmonic umbrella potential or constant force */
 910 static void do_pull_pot(int ePull,
 911                         t_pull *pull, t_pbc *pbc, double t, real lambda,
 912                         real *V, tensor vir, real *dVdl)
 913 {
 914     int        g, j, m;
 915     dvec       dev;
 916     double     ndr, invdr;
 917     real       k, dkdl;
 918     t_pullgrp *pgrp;
 919
 920     /* loop over the groups that are being pulled */
 921     *V    = 0;
 922     *dVdl = 0;
 923     for (g = 1; g < 1+pull->ngrp; g++)
 924     {
 925         pgrp = &pull->grp[g];
 926         get_pullgrp_distance(pull, pbc, g, t, pgrp->dr, dev);
 927
 928         k    = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB;
 929         dkdl = pgrp->kB - pgrp->k;
 930
 931         switch (pull->eGeom)
 932         {
 933             case epullgDIST:
 934                 ndr   = dnorm(pgrp->dr);
 935                 invdr = 1/ndr;
 936                 if (ePull == epullUMBRELLA)
 937                 {
 938                     pgrp->f_scal  =       -k*dev[0];
 939                     *V           += 0.5*   k*dsqr(dev[0]);
 940                     *dVdl        += 0.5*dkdl*dsqr(dev[0]);
 941                 }
 942                 else
 943                 {
 944                     pgrp->f_scal  =   -k;
 945                     *V           +=    k*ndr;
 946                     *dVdl        += dkdl*ndr;
 947                 }
 948                 for (m = 0; m < DIM; m++)
 949                 {
 950                     pgrp->f[m]    = pgrp->f_scal*pgrp->dr[m]*invdr;
 951                 }
 952                 break;
 953             case epullgDIR:
 954             case epullgDIRPBC:
 955             case epullgCYL:
 956                 if (ePull == epullUMBRELLA)
 957                 {
 958                     pgrp->f_scal  =       -k*dev[0];
 959                     *V           += 0.5*   k*dsqr(dev[0]);
 960                     *dVdl        += 0.5*dkdl*dsqr(dev[0]);
 961                 }
 962                 else
 963                 {
 964                     ndr = 0;
 965                     for (m = 0; m < DIM; m++)
 966                     {
 967                         ndr += pgrp->vec[m]*pgrp->dr[m];
 968                     }
 969                     pgrp->f_scal  =   -k;
 970                     *V           +=    k*ndr;
 971                     *dVdl        += dkdl*ndr;
 972                 }
 973                 for (m = 0; m < DIM; m++)
 974                 {
 975                     pgrp->f[m]    = pgrp->f_scal*pgrp->vec[m];
 976                 }
 977                 break;
 978             case epullgPOS:
 979                 for (m = 0; m < DIM; m++)
 980                 {
 981                     if (ePull == epullUMBRELLA)
 982                     {
 983                         pgrp->f[m]  =       -k*dev[m];
 984                         *V         += 0.5*   k*dsqr(dev[m]);
 985                         *dVdl      += 0.5*dkdl*dsqr(dev[m]);
 986                     }
 987                     else
 988                     {
 989                         pgrp->f[m]  =   -k*pull->dim[m];
 990                         *V         +=    k*pgrp->dr[m]*pull->dim[m];
 991                         *dVdl      += dkdl*pgrp->dr[m]*pull->dim[m];
 992                     }
 993                 }
 994                 break;
 995         }
 996
 997         if (vir)
 998         {
 999             /* Add the pull contribution to the virial */
1000             for (j = 0; j < DIM; j++)
1001             {
1002                 for (m = 0; m < DIM; m++)
1003                 {
1004                     vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m];
1005                 }
1006             }
1007         }
1008     }
1009 }
1010
1011 real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1012                     t_commrec *cr, double t, real lambda,
1013                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
1014 {
1015     real V, dVdl;
1016
1017     pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1018
1019     do_pull_pot(ePull, pull, pbc, t, lambda,
1020                 &V, pull->bVirial && MASTER(cr) ? vir : NULL, &dVdl);
1021
1022     /* Distribute forces over pulled groups */
1023     apply_forces(pull, md, DOMAINDECOMP(cr) ? cr->dd->ga2la : NULL, f);
1024
1025     if (MASTER(cr))
1026     {
1027         *dvdlambda += dVdl;
1028     }
1029
1030     return (MASTER(cr) ? V : 0.0);
1031 }
1032
1033 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1034                      t_commrec *cr, double dt, double t,
1035                      rvec *x, rvec *xp, rvec *v, tensor vir)
1036 {
1037     pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1038
1039     do_constraint(pull, md, pbc, xp, v, pull->bVirial && MASTER(cr), vir, dt, t);
1040 }
1041
1042 static void make_local_pull_group(gmx_ga2la_t ga2la,
1043                                   t_pullgrp *pg, int start, int end)
1044 {
1045     int i, ii;
1046
1047     pg->nat_loc = 0;
1048     for (i = 0; i < pg->nat; i++)
1049     {
1050         ii = pg->ind[i];
1051         if (ga2la)
1052         {
1053             if (!ga2la_get_home(ga2la, ii, &ii))
1054             {
1055                 ii = -1;
1056             }
1057         }
1058         if (ii >= start && ii < end)
1059         {
1060             /* This is a home atom, add it to the local pull group */
1061             if (pg->nat_loc >= pg->nalloc_loc)
1062             {
1063                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1064                 srenew(pg->ind_loc, pg->nalloc_loc);
1065                 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
1066                 {
1067                     srenew(pg->weight_loc, pg->nalloc_loc);
1068                 }
1069             }
1070             pg->ind_loc[pg->nat_loc] = ii;
1071             if (pg->weight)
1072             {
1073                 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1074             }
1075             pg->nat_loc++;
1076         }
1077     }
1078 }
1079
1080 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
1081 {
1082     gmx_ga2la_t ga2la;
1083     int         g;
1084
1085     if (dd)
1086     {
1087         ga2la = dd->ga2la;
1088     }
1089     else
1090     {
1091         ga2la = NULL;
1092     }
1093
1094     if (pull->grp[0].nat > 0)
1095     {
1096         make_local_pull_group(ga2la, &pull->grp[0], md->start, md->start+md->homenr);
1097     }
1098     for (g = 1; g < 1+pull->ngrp; g++)
1099     {
1100         make_local_pull_group(ga2la, &pull->grp[g], md->start, md->start+md->homenr);
1101     }
1102 }
1103
1104 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1105                                   int start, int end,
1106                                   int g, t_pullgrp *pg, ivec pulldims,
1107                                   gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
1108 {
1109     int                   i, ii, d, nfrozen, ndim;
1110     real                  m, w, mbd;
1111     double                tmass, wmass, wwmass;
1112     gmx_bool              bDomDec;
1113     gmx_ga2la_t           ga2la = NULL;
1114     gmx_groups_t         *groups;
1115     gmx_mtop_atomlookup_t alook;
1116     t_atom               *atom;
1117
1118     bDomDec = (cr && DOMAINDECOMP(cr));
1119     if (bDomDec)
1120     {
1121         ga2la = cr->dd->ga2la;
1122     }
1123
1124     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1125     {
1126         /* There are no masses in the integrator.
1127          * But we still want to have the correct mass-weighted COMs.
1128          * So we store the real masses in the weights.
1129          * We do not set nweight, so these weights do not end up in the tpx file.
1130          */
1131         if (pg->nweight == 0)
1132         {
1133             snew(pg->weight, pg->nat);
1134         }
1135     }
1136
1137     if (cr && PAR(cr))
1138     {
1139         pg->nat_loc    = 0;
1140         pg->nalloc_loc = 0;
1141         pg->ind_loc    = NULL;
1142         pg->weight_loc = NULL;
1143     }
1144     else
1145     {
1146         pg->nat_loc = pg->nat;
1147         pg->ind_loc = pg->ind;
1148         if (pg->epgrppbc == epgrppbcCOS)
1149         {
1150             snew(pg->weight_loc, pg->nat);
1151         }
1152         else
1153         {
1154             pg->weight_loc = pg->weight;
1155         }
1156     }
1157
1158     groups = &mtop->groups;
1159
1160     alook = gmx_mtop_atomlookup_init(mtop);
1161
1162     nfrozen = 0;
1163     tmass   = 0;
1164     wmass   = 0;
1165     wwmass  = 0;
1166     for (i = 0; i < pg->nat; i++)
1167     {
1168         ii = pg->ind[i];
1169         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1170         if (cr && PAR(cr) && !bDomDec && ii >= start && ii < end)
1171         {
1172             pg->ind_loc[pg->nat_loc++] = ii;
1173         }
1174         if (ir->opts.nFreeze)
1175         {
1176             for (d = 0; d < DIM; d++)
1177             {
1178                 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1179                 {
1180                     nfrozen++;
1181                 }
1182             }
1183         }
1184         if (ir->efep == efepNO)
1185         {
1186             m = atom->m;
1187         }
1188         else
1189         {
1190             m = (1 - lambda)*atom->m + lambda*atom->mB;
1191         }
1192         if (pg->nweight > 0)
1193         {
1194             w = pg->weight[i];
1195         }
1196         else
1197         {
1198             w = 1;
1199         }
1200         if (EI_ENERGY_MINIMIZATION(ir->eI))
1201         {
1202             /* Move the mass to the weight */
1203             w            *= m;
1204             m             = 1;
1205             pg->weight[i] = w;
1206         }
1207         else if (ir->eI == eiBD)
1208         {
1209             if (ir->bd_fric)
1210             {
1211                 mbd = ir->bd_fric*ir->delta_t;
1212             }
1213             else
1214             {
1215                 if (groups->grpnr[egcTC] == NULL)
1216                 {
1217                     mbd = ir->delta_t/ir->opts.tau_t[0];
1218                 }
1219                 else
1220                 {
1221                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1222                 }
1223             }
1224             w            *= m/mbd;
1225             m             = mbd;
1226             pg->weight[i] = w;
1227         }
1228         tmass  += m;
1229         wmass  += m*w;
1230         wwmass += m*w*w;
1231     }
1232
1233     gmx_mtop_atomlookup_destroy(alook);
1234
1235     if (wmass == 0)
1236     {
1237         gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1238                   pg->weight ? " weighted" : "", g);
1239     }
1240     if (fplog)
1241     {
1242         fprintf(fplog,
1243                 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1244         if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1245         {
1246             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1247         }
1248         if (pg->epgrppbc == epgrppbcCOS)
1249         {
1250             fprintf(fplog, ", cosine weighting will be used");
1251         }
1252         fprintf(fplog, "\n");
1253     }
1254
1255     if (nfrozen == 0)
1256     {
1257         /* A value > 0 signals not frozen, it is updated later */
1258         pg->invtm  = 1.0;
1259     }
1260     else
1261     {
1262         ndim = 0;
1263         for (d = 0; d < DIM; d++)
1264         {
1265             ndim += pulldims[d]*pg->nat;
1266         }
1267         if (fplog && nfrozen > 0 && nfrozen < ndim)
1268         {
1269             fprintf(fplog,
1270                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1271                     "         that are subject to pulling are frozen.\n"
1272                     "         For pulling the whole group will be frozen.\n\n",
1273                     g);
1274         }
1275         pg->invtm  = 0.0;
1276         pg->wscale = 1.0;
1277     }
1278 }
1279
1280 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1281                gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1282                gmx_bool bOutFile, unsigned long Flags)
1283 {
1284     t_pull       *pull;
1285     t_pullgrp    *pgrp;
1286     int           g, start = 0, end = 0, m;
1287     gmx_bool      bCite;
1288
1289     pull = ir->pull;
1290
1291     pull->ePBC = ir->ePBC;
1292     switch (pull->ePBC)
1293     {
1294         case epbcNONE: pull->npbcdim = 0; break;
1295         case epbcXY:   pull->npbcdim = 2; break;
1296         default:       pull->npbcdim = 3; break;
1297     }
1298
1299     if (fplog)
1300     {
1301         fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
1302                 EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
1303         if (pull->grp[0].nat > 0)
1304         {
1305             fprintf(fplog, "between a reference group and %d group%s\n",
1306                     pull->ngrp, pull->ngrp == 1 ? "" : "s");
1307         }
1308         else
1309         {
1310             fprintf(fplog, "with an absolute reference on %d group%s\n",
1311                     pull->ngrp, pull->ngrp == 1 ? "" : "s");
1312         }
1313         bCite = FALSE;
1314         for (g = 0; g < pull->ngrp+1; g++)
1315         {
1316             if (pull->grp[g].nat > 1 &&
1317                 pull->grp[g].pbcatom < 0)
1318             {
1319                 /* We are using cosine weighting */
1320                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1321                 bCite = TRUE;
1322             }
1323         }
1324         if (bCite)
1325         {
1326             please_cite(fplog, "Engin2010");
1327         }
1328     }
1329
1330     /* We always add the virial contribution,
1331      * except for geometry = direction_periodic where this is impossible.
1332      */
1333     pull->bVirial = (pull->eGeom != epullgDIRPBC);
1334     if (getenv("GMX_NO_PULLVIR") != NULL)
1335     {
1336         if (fplog)
1337         {
1338             fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
1339         }
1340         pull->bVirial = FALSE;
1341     }
1342
1343     if (cr && PARTDECOMP(cr))
1344     {
1345         pd_at_range(cr, &start, &end);
1346     }
1347     pull->rbuf     = NULL;
1348     pull->dbuf     = NULL;
1349     pull->dbuf_cyl = NULL;
1350     pull->bRefAt   = FALSE;
1351     pull->cosdim   = -1;
1352     for (g = 0; g < pull->ngrp+1; g++)
1353     {
1354         pgrp           = &pull->grp[g];
1355         pgrp->epgrppbc = epgrppbcNONE;
1356         if (pgrp->nat > 0)
1357         {
1358             /* Determine if we need to take PBC into account for calculating
1359              * the COM's of the pull groups.
1360              */
1361             for (m = 0; m < pull->npbcdim; m++)
1362             {
1363                 if (pull->dim[m] && pgrp->nat > 1)
1364                 {
1365                     if (pgrp->pbcatom >= 0)
1366                     {
1367                         pgrp->epgrppbc = epgrppbcREFAT;
1368                         pull->bRefAt   = TRUE;
1369                     }
1370                     else
1371                     {
1372                         if (pgrp->weight)
1373                         {
1374                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1375                         }
1376                         pgrp->epgrppbc = epgrppbcCOS;
1377                         if (pull->cosdim >= 0 && pull->cosdim != m)
1378                         {
1379                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1380                         }
1381                         pull->cosdim = m;
1382                     }
1383                 }
1384             }
1385             /* Set the indices */
1386             init_pull_group_index(fplog, cr, start, end, g, pgrp, pull->dim, mtop, ir, lambda);
1387             if (PULL_CYL(pull) && pgrp->invtm == 0)
1388             {
1389                 gmx_fatal(FARGS, "Can not have frozen atoms in a cylinder pull group");
1390             }
1391         }
1392         else
1393         {
1394             /* Absolute reference, set the inverse mass to zero */
1395             pgrp->invtm  = 0;
1396             pgrp->wscale = 1;
1397         }
1398     }
1399
1400     /* if we use dynamic reference groups, do some initialising for them */
1401     if (PULL_CYL(pull))
1402     {
1403         if (pull->grp[0].nat == 0)
1404         {
1405             gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
1406         }
1407         snew(pull->dyna, pull->ngrp+1);
1408     }
1409
1410     /* Only do I/O when we are doing dynamics and if we are the MASTER */
1411     pull->out_x = NULL;
1412     pull->out_f = NULL;
1413     if (bOutFile)
1414     {
1415         if (pull->nstxout > 0)
1416         {
1417             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags);
1418         }
1419         if (pull->nstfout > 0)
1420         {
1421             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1422                                         FALSE, Flags);
1423         }
1424     }
1425 }
1426
1427 void finish_pull(FILE *fplog, t_pull *pull)
1428 {
1429     if (pull->out_x)
1430     {
1431         gmx_fio_fclose(pull->out_x);
1432     }
1433     if (pull->out_f)
1434     {
1435         gmx_fio_fclose(pull->out_f);
1436     }
1437 }