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             get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0],
 692                             -1, unc_ij);
 693
 694             switch (pull->eGeom)
 695             {
 696                 case epullgDIST:
 697                     bConverged = fabs(dnorm(unc_ij) - ref[0]) < pull->constr_tol;
 698                     break;
 699                 case epullgDIR:
 700                 case epullgDIRPBC:
 701                 case epullgCYL:
 702                     for (m = 0; m < DIM; m++)
 703                     {
 704                         vec[m] = pgrp->vec[m];
 705                     }
 706                     inpr = diprod(unc_ij, vec);
 707                     dsvmul(inpr, vec, unc_ij);
 708                     bConverged =
 709                         fabs(diprod(unc_ij, vec) - ref[0]) < pull->constr_tol;
 710                     break;
 711                 case epullgPOS:
 712                     bConverged = TRUE;
 713                     for (m = 0; m < DIM; m++)
 714                     {
 715                         if (pull->dim[m] &&
 716                             fabs(unc_ij[m] - ref[m]) >= pull->constr_tol)
 717                         {
 718                             bConverged = FALSE;
 719                         }
 720                     }
 721                     break;
 722             }
 723
 724             if (!bConverged)
 725             {
 726                 if (debug)
 727                 {
 728                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
 729                             "d_ref = %f %f %f, current d = %f\n",
 730                             g, ref[0], ref[1], ref[2], dnorm(unc_ij));
 731                 }
 732
 733                 bConverged_all = FALSE;
 734             }
 735         }
 736
 737         niter++;
 738         /* if after all constraints are dealt with and bConverged is still TRUE
 739            we're finished, if not we do another iteration */
 740     }
 741     if (niter > max_iter)
 742     {
 743         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
 744     }
 745
 746     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 747
 748     if (v)
 749     {
 750         invdt = 1/dt;
 751     }
 752
 753     /* update the normal groups */
 754     for (g = 1; g < 1+pull->ngrp; g++)
 755     {
 756         pgrp = &pull->grp[g];
 757         /* get the final dr and constraint force for group i */
 758         dvec_sub(rinew[g], pgrp->xp, dr[g]);
 759         /* select components of dr */
 760         for (m = 0; m < DIM; m++)
 761         {
 762             dr[g][m] *= pull->dim[m];
 763         }
 764         dsvmul(1.0/(pgrp->invtm*dt*dt), dr[g], f);
 765         dvec_inc(pgrp->f, f);
 766         switch (pull->eGeom)
 767         {
 768             case epullgDIST:
 769                 for (m = 0; m < DIM; m++)
 770                 {
 771                     pgrp->f_scal += r_ij[g][m]*f[m]/dnorm(r_ij[g]);
 772                 }
 773                 break;
 774             case epullgDIR:
 775             case epullgDIRPBC:
 776             case epullgCYL:
 777                 for (m = 0; m < DIM; m++)
 778                 {
 779                     pgrp->f_scal += pgrp->vec[m]*f[m];
 780                 }
 781                 break;
 782             case epullgPOS:
 783                 break;
 784         }
 785
 786         if (vir && bMaster)
 787         {
 788             /* Add the pull contribution to the virial */
 789             for (j = 0; j < DIM; j++)
 790             {
 791                 for (m = 0; m < DIM; m++)
 792                 {
 793                     vir[j][m] -= 0.5*f[j]*r_ij[g][m];
 794                 }
 795             }
 796         }
 797
 798         /* update the atom positions */
 799         copy_dvec(dr[g], tmp);
 800         for (j = 0; j < pgrp->nat_loc; j++)
 801         {
 802             ii = pgrp->ind_loc[j];
 803             if (pgrp->weight_loc)
 804             {
 805                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr[g], tmp);
 806             }
 807             for (m = 0; m < DIM; m++)
 808             {
 809                 x[ii][m] += tmp[m];
 810             }
 811             if (v)
 812             {
 813                 for (m = 0; m < DIM; m++)
 814                 {
 815                     v[ii][m] += invdt*tmp[m];
 816                 }
 817             }
 818         }
 819     }
 820
 821     /* update the reference groups */
 822     if (PULL_CYL(pull))
 823     {
 824         /* update the dynamic reference groups */
 825         for (g = 1; g < 1+pull->ngrp; g++)
 826         {
 827             pdyna = &pull->dyna[g];
 828             dvec_sub(rjnew[g], pdyna->xp, ref_dr);
 829             /* select components of ref_dr */
 830             for (m = 0; m < DIM; m++)
 831             {
 832                 ref_dr[m] *= pull->dim[m];
 833             }
 834
 835             for (j = 0; j < pdyna->nat_loc; j++)
 836             {
 837                 /* reset the atoms with dr, weighted by w_i */
 838                 dsvmul(pdyna->wscale*pdyna->weight_loc[j], ref_dr, tmp);
 839                 ii = pdyna->ind_loc[j];
 840                 for (m = 0; m < DIM; m++)
 841                 {
 842                     x[ii][m] += tmp[m];
 843                 }
 844                 if (v)
 845                 {
 846                     for (m = 0; m < DIM; m++)
 847                     {
 848                         v[ii][m] += invdt*tmp[m];
 849                     }
 850                 }
 851             }
 852         }
 853     }
 854     else
 855     {
 856         pgrp = &pull->grp[0];
 857         /* update the reference group */
 858         dvec_sub(rjnew[0], pgrp->xp, ref_dr);
 859         /* select components of ref_dr */
 860         for (m = 0; m < DIM; m++)
 861         {
 862             ref_dr[m] *= pull->dim[m];
 863         }
 864
 865         copy_dvec(ref_dr, tmp);
 866         for (j = 0; j < pgrp->nat_loc; j++)
 867         {
 868             ii = pgrp->ind_loc[j];
 869             if (pgrp->weight_loc)
 870             {
 871                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], ref_dr, tmp);
 872             }
 873             for (m = 0; m < DIM; m++)
 874             {
 875                 x[ii][m] += tmp[m];
 876             }
 877             if (v)
 878             {
 879                 for (m = 0; m < DIM; m++)
 880                 {
 881                     v[ii][m] += invdt*tmp[m];
 882                 }
 883             }
 884         }
 885     }
 886
 887     /* finished! I hope. Give back some memory */
 888     sfree(r_ij);
 889     sfree(rinew);
 890     sfree(rjnew);
 891     sfree(dr);
 892 }
 893
 894 /* Pulling with a harmonic umbrella potential or constant force */
 895 static void do_pull_pot(int ePull,
 896                         t_pull *pull, t_pbc *pbc, double t, real lambda,
 897                         real *V, tensor vir, real *dVdl)
 898 {
 899     int        g, j, m;
 900     dvec       dev;
 901     double     ndr, invdr;
 902     real       k, dkdl;
 903     t_pullgrp *pgrp;
 904
 905     /* loop over the groups that are being pulled */
 906     *V    = 0;
 907     *dVdl = 0;
 908     for (g = 1; g < 1+pull->ngrp; g++)
 909     {
 910         pgrp = &pull->grp[g];
 911         get_pullgrp_distance(pull, pbc, g, t, pgrp->dr, dev);
 912
 913         k    = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB;
 914         dkdl = pgrp->kB - pgrp->k;
 915
 916         switch (pull->eGeom)
 917         {
 918             case epullgDIST:
 919                 ndr   = dnorm(pgrp->dr);
 920                 invdr = 1/ndr;
 921                 if (ePull == epullUMBRELLA)
 922                 {
 923                     pgrp->f_scal  =       -k*dev[0];
 924                     *V           += 0.5*   k*dsqr(dev[0]);
 925                     *dVdl        += 0.5*dkdl*dsqr(dev[0]);
 926                 }
 927                 else
 928                 {
 929                     pgrp->f_scal  =   -k;
 930                     *V           +=    k*ndr;
 931                     *dVdl        += dkdl*ndr;
 932                 }
 933                 for (m = 0; m < DIM; m++)
 934                 {
 935                     pgrp->f[m]    = pgrp->f_scal*pgrp->dr[m]*invdr;
 936                 }
 937                 break;
 938             case epullgDIR:
 939             case epullgDIRPBC:
 940             case epullgCYL:
 941                 if (ePull == epullUMBRELLA)
 942                 {
 943                     pgrp->f_scal  =       -k*dev[0];
 944                     *V           += 0.5*   k*dsqr(dev[0]);
 945                     *dVdl        += 0.5*dkdl*dsqr(dev[0]);
 946                 }
 947                 else
 948                 {
 949                     ndr = 0;
 950                     for (m = 0; m < DIM; m++)
 951                     {
 952                         ndr += pgrp->vec[m]*pgrp->dr[m];
 953                     }
 954                     pgrp->f_scal  =   -k;
 955                     *V           +=    k*ndr;
 956                     *dVdl        += dkdl*ndr;
 957                 }
 958                 for (m = 0; m < DIM; m++)
 959                 {
 960                     pgrp->f[m]    = pgrp->f_scal*pgrp->vec[m];
 961                 }
 962                 break;
 963             case epullgPOS:
 964                 for (m = 0; m < DIM; m++)
 965                 {
 966                     if (ePull == epullUMBRELLA)
 967                     {
 968                         pgrp->f[m]  =       -k*dev[m];
 969                         *V         += 0.5*   k*dsqr(dev[m]);
 970                         *dVdl      += 0.5*dkdl*dsqr(dev[m]);
 971                     }
 972                     else
 973                     {
 974                         pgrp->f[m]  =   -k*pull->dim[m];
 975                         *V         +=    k*pgrp->dr[m]*pull->dim[m];
 976                         *dVdl      += dkdl*pgrp->dr[m]*pull->dim[m];
 977                     }
 978                 }
 979                 break;
 980         }
 981
 982         if (vir)
 983         {
 984             /* Add the pull contribution to the virial */
 985             for (j = 0; j < DIM; j++)
 986             {
 987                 for (m = 0; m < DIM; m++)
 988                 {
 989                     vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m];
 990                 }
 991             }
 992         }
 993     }
 994 }
 995
 996 real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
 997                     t_commrec *cr, double t, real lambda,
 998                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
 999 {
1000     real V, dVdl;
1001
1002     pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1003
1004     do_pull_pot(ePull, pull, pbc, t, lambda,
1005                 &V, pull->bVirial && MASTER(cr) ? vir : NULL, &dVdl);
1006
1007     /* Distribute forces over pulled groups */
1008     apply_forces(pull, md, DOMAINDECOMP(cr) ? cr->dd->ga2la : NULL, f);
1009
1010     if (MASTER(cr))
1011     {
1012         *dvdlambda += dVdl;
1013     }
1014
1015     return (MASTER(cr) ? V : 0.0);
1016 }
1017
1018 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1019                      t_commrec *cr, double dt, double t,
1020                      rvec *x, rvec *xp, rvec *v, tensor vir)
1021 {
1022     pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1023
1024     do_constraint(pull, md, pbc, xp, v, pull->bVirial && MASTER(cr), vir, dt, t);
1025 }
1026
1027 static void make_local_pull_group(gmx_ga2la_t ga2la,
1028                                   t_pullgrp *pg, int start, int end)
1029 {
1030     int i, ii;
1031
1032     pg->nat_loc = 0;
1033     for (i = 0; i < pg->nat; i++)
1034     {
1035         ii = pg->ind[i];
1036         if (ga2la)
1037         {
1038             if (!ga2la_get_home(ga2la, ii, &ii))
1039             {
1040                 ii = -1;
1041             }
1042         }
1043         if (ii >= start && ii < end)
1044         {
1045             /* This is a home atom, add it to the local pull group */
1046             if (pg->nat_loc >= pg->nalloc_loc)
1047             {
1048                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1049                 srenew(pg->ind_loc, pg->nalloc_loc);
1050                 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
1051                 {
1052                     srenew(pg->weight_loc, pg->nalloc_loc);
1053                 }
1054             }
1055             pg->ind_loc[pg->nat_loc] = ii;
1056             if (pg->weight)
1057             {
1058                 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1059             }
1060             pg->nat_loc++;
1061         }
1062     }
1063 }
1064
1065 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
1066 {
1067     gmx_ga2la_t ga2la;
1068     int         g;
1069
1070     if (dd)
1071     {
1072         ga2la = dd->ga2la;
1073     }
1074     else
1075     {
1076         ga2la = NULL;
1077     }
1078
1079     if (pull->grp[0].nat > 0)
1080     {
1081         make_local_pull_group(ga2la, &pull->grp[0], md->start, md->start+md->homenr);
1082     }
1083     for (g = 1; g < 1+pull->ngrp; g++)
1084     {
1085         make_local_pull_group(ga2la, &pull->grp[g], md->start, md->start+md->homenr);
1086     }
1087 }
1088
1089 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1090                                   int start, int end,
1091                                   int g, t_pullgrp *pg, ivec pulldims,
1092                                   gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
1093 {
1094     int                   i, ii, d, nfrozen, ndim;
1095     real                  m, w, mbd;
1096     double                tmass, wmass, wwmass;
1097     gmx_bool              bDomDec;
1098     gmx_ga2la_t           ga2la = NULL;
1099     gmx_groups_t         *groups;
1100     gmx_mtop_atomlookup_t alook;
1101     t_atom               *atom;
1102
1103     bDomDec = (cr && DOMAINDECOMP(cr));
1104     if (bDomDec)
1105     {
1106         ga2la = cr->dd->ga2la;
1107     }
1108
1109     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1110     {
1111         /* There are no masses in the integrator.
1112          * But we still want to have the correct mass-weighted COMs.
1113          * So we store the real masses in the weights.
1114          * We do not set nweight, so these weights do not end up in the tpx file.
1115          */
1116         if (pg->nweight == 0)
1117         {
1118             snew(pg->weight, pg->nat);
1119         }
1120     }
1121
1122     if (cr && PAR(cr))
1123     {
1124         pg->nat_loc    = 0;
1125         pg->nalloc_loc = 0;
1126         pg->ind_loc    = NULL;
1127         pg->weight_loc = NULL;
1128     }
1129     else
1130     {
1131         pg->nat_loc = pg->nat;
1132         pg->ind_loc = pg->ind;
1133         if (pg->epgrppbc == epgrppbcCOS)
1134         {
1135             snew(pg->weight_loc, pg->nat);
1136         }
1137         else
1138         {
1139             pg->weight_loc = pg->weight;
1140         }
1141     }
1142
1143     groups = &mtop->groups;
1144
1145     alook = gmx_mtop_atomlookup_init(mtop);
1146
1147     nfrozen = 0;
1148     tmass   = 0;
1149     wmass   = 0;
1150     wwmass  = 0;
1151     for (i = 0; i < pg->nat; i++)
1152     {
1153         ii = pg->ind[i];
1154         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1155         if (cr && PAR(cr) && !bDomDec && ii >= start && ii < end)
1156         {
1157             pg->ind_loc[pg->nat_loc++] = ii;
1158         }
1159         if (ir->opts.nFreeze)
1160         {
1161             for (d = 0; d < DIM; d++)
1162             {
1163                 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1164                 {
1165                     nfrozen++;
1166                 }
1167             }
1168         }
1169         if (ir->efep == efepNO)
1170         {
1171             m = atom->m;
1172         }
1173         else
1174         {
1175             m = (1 - lambda)*atom->m + lambda*atom->mB;
1176         }
1177         if (pg->nweight > 0)
1178         {
1179             w = pg->weight[i];
1180         }
1181         else
1182         {
1183             w = 1;
1184         }
1185         if (EI_ENERGY_MINIMIZATION(ir->eI))
1186         {
1187             /* Move the mass to the weight */
1188             w            *= m;
1189             m             = 1;
1190             pg->weight[i] = w;
1191         }
1192         else if (ir->eI == eiBD)
1193         {
1194             if (ir->bd_fric)
1195             {
1196                 mbd = ir->bd_fric*ir->delta_t;
1197             }
1198             else
1199             {
1200                 if (groups->grpnr[egcTC] == NULL)
1201                 {
1202                     mbd = ir->delta_t/ir->opts.tau_t[0];
1203                 }
1204                 else
1205                 {
1206                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1207                 }
1208             }
1209             w            *= m/mbd;
1210             m             = mbd;
1211             pg->weight[i] = w;
1212         }
1213         tmass  += m;
1214         wmass  += m*w;
1215         wwmass += m*w*w;
1216     }
1217
1218     gmx_mtop_atomlookup_destroy(alook);
1219
1220     if (wmass == 0)
1221     {
1222         gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1223                   pg->weight ? " weighted" : "", g);
1224     }
1225     if (fplog)
1226     {
1227         fprintf(fplog,
1228                 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1229         if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1230         {
1231             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1232         }
1233         if (pg->epgrppbc == epgrppbcCOS)
1234         {
1235             fprintf(fplog, ", cosine weighting will be used");
1236         }
1237         fprintf(fplog, "\n");
1238     }
1239
1240     if (nfrozen == 0)
1241     {
1242         /* A value > 0 signals not frozen, it is updated later */
1243         pg->invtm  = 1.0;
1244     }
1245     else
1246     {
1247         ndim = 0;
1248         for (d = 0; d < DIM; d++)
1249         {
1250             ndim += pulldims[d]*pg->nat;
1251         }
1252         if (fplog && nfrozen > 0 && nfrozen < ndim)
1253         {
1254             fprintf(fplog,
1255                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1256                     "         that are subject to pulling are frozen.\n"
1257                     "         For pulling the whole group will be frozen.\n\n",
1258                     g);
1259         }
1260         pg->invtm  = 0.0;
1261         pg->wscale = 1.0;
1262     }
1263 }
1264
1265 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1266                gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1267                gmx_bool bOutFile, unsigned long Flags)
1268 {
1269     t_pull       *pull;
1270     t_pullgrp    *pgrp;
1271     int           g, start = 0, end = 0, m;
1272     gmx_bool      bCite;
1273
1274     pull = ir->pull;
1275
1276     pull->ePBC = ir->ePBC;
1277     switch (pull->ePBC)
1278     {
1279         case epbcNONE: pull->npbcdim = 0; break;
1280         case epbcXY:   pull->npbcdim = 2; break;
1281         default:       pull->npbcdim = 3; break;
1282     }
1283
1284     if (fplog)
1285     {
1286         fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
1287                 EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
1288         if (pull->grp[0].nat > 0)
1289         {
1290             fprintf(fplog, "between a reference group and %d group%s\n",
1291                     pull->ngrp, pull->ngrp == 1 ? "" : "s");
1292         }
1293         else
1294         {
1295             fprintf(fplog, "with an absolute reference on %d group%s\n",
1296                     pull->ngrp, pull->ngrp == 1 ? "" : "s");
1297         }
1298         bCite = FALSE;
1299         for (g = 0; g < pull->ngrp+1; g++)
1300         {
1301             if (pull->grp[g].nat > 1 &&
1302                 pull->grp[g].pbcatom < 0)
1303             {
1304                 /* We are using cosine weighting */
1305                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1306                 bCite = TRUE;
1307             }
1308         }
1309         if (bCite)
1310         {
1311             please_cite(fplog, "Engin2010");
1312         }
1313     }
1314
1315     /* We always add the virial contribution,
1316      * except for geometry = direction_periodic where this is impossible.
1317      */
1318     pull->bVirial = (pull->eGeom != epullgDIRPBC);
1319     if (getenv("GMX_NO_PULLVIR") != NULL)
1320     {
1321         if (fplog)
1322         {
1323             fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
1324         }
1325         pull->bVirial = FALSE;
1326     }
1327
1328     if (cr && PARTDECOMP(cr))
1329     {
1330         pd_at_range(cr, &start, &end);
1331     }
1332     pull->rbuf     = NULL;
1333     pull->dbuf     = NULL;
1334     pull->dbuf_cyl = NULL;
1335     pull->bRefAt   = FALSE;
1336     pull->cosdim   = -1;
1337     for (g = 0; g < pull->ngrp+1; g++)
1338     {
1339         pgrp           = &pull->grp[g];
1340         pgrp->epgrppbc = epgrppbcNONE;
1341         if (pgrp->nat > 0)
1342         {
1343             /* Determine if we need to take PBC into account for calculating
1344              * the COM's of the pull groups.
1345              */
1346             for (m = 0; m < pull->npbcdim; m++)
1347             {
1348                 if (pull->dim[m] && pgrp->nat > 1)
1349                 {
1350                     if (pgrp->pbcatom >= 0)
1351                     {
1352                         pgrp->epgrppbc = epgrppbcREFAT;
1353                         pull->bRefAt   = TRUE;
1354                     }
1355                     else
1356                     {
1357                         if (pgrp->weight)
1358                         {
1359                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1360                         }
1361                         pgrp->epgrppbc = epgrppbcCOS;
1362                         if (pull->cosdim >= 0 && pull->cosdim != m)
1363                         {
1364                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1365                         }
1366                         pull->cosdim = m;
1367                     }
1368                 }
1369             }
1370             /* Set the indices */
1371             init_pull_group_index(fplog, cr, start, end, g, pgrp, pull->dim, mtop, ir, lambda);
1372             if (PULL_CYL(pull) && pgrp->invtm == 0)
1373             {
1374                 gmx_fatal(FARGS, "Can not have frozen atoms in a cylinder pull group");
1375             }
1376         }
1377         else
1378         {
1379             /* Absolute reference, set the inverse mass to zero */
1380             pgrp->invtm  = 0;
1381             pgrp->wscale = 1;
1382         }
1383     }
1384
1385     /* if we use dynamic reference groups, do some initialising for them */
1386     if (PULL_CYL(pull))
1387     {
1388         if (pull->grp[0].nat == 0)
1389         {
1390             gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
1391         }
1392         snew(pull->dyna, pull->ngrp+1);
1393     }
1394
1395     /* Only do I/O when we are doing dynamics and if we are the MASTER */
1396     pull->out_x = NULL;
1397     pull->out_f = NULL;
1398     if (bOutFile)
1399     {
1400         if (pull->nstxout > 0)
1401         {
1402             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags);
1403         }
1404         if (pull->nstfout > 0)
1405         {
1406             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1407                                         FALSE, Flags);
1408         }
1409     }
1410 }
1411
1412 void finish_pull(FILE *fplog, t_pull *pull)
1413 {
1414     if (pull->out_x)
1415     {
1416         gmx_fio_fclose(pull->out_x);
1417     }
1418     if (pull->out_f)
1419     {
1420         gmx_fio_fclose(pull->out_f);
1421     }
1422 }