src/gromacs/mdlib/pull.c

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