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