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