src/gromacs/pulling/pull.cpp

   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 "config.h"
  42
  43 #include <assert.h>
  44 #include <math.h>
  45 #include <stdio.h>
  46 #include <stdlib.h>
  47 #include <string.h>
  48
  49 #include <algorithm>
  50
  51 #include "gromacs/fileio/filenm.h"
  52 #include "gromacs/fileio/gmxfio.h"
  53 #include "gromacs/fileio/xvgr.h"
  54 #include "gromacs/legacyheaders/copyrite.h"
  55 #include "gromacs/legacyheaders/gmx_ga2la.h"
  56 #include "gromacs/legacyheaders/macros.h"
  57 #include "gromacs/legacyheaders/mdrun.h"
  58 #include "gromacs/legacyheaders/names.h"
  59 #include "gromacs/legacyheaders/network.h"
  60 #include "gromacs/legacyheaders/typedefs.h"
  61 #include "gromacs/legacyheaders/types/commrec.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/pbcutil/pbc.h"
  64 #include "gromacs/pulling/pull_internal.h"
  65 #include "gromacs/topology/mtop_util.h"
  66 #include "gromacs/utility/futil.h"
  67 #include "gromacs/utility/gmxassert.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 static void pull_print_group_x(FILE *out, ivec dim,
  71                                const pull_group_work_t *pgrp)
  72 {
  73     int m;
  74
  75     for (m = 0; m < DIM; m++)
  76     {
  77         if (dim[m])
  78         {
  79             fprintf(out, "\t%g", pgrp->x[m]);
  80         }
  81     }
  82 }
  83
  84 static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
  85                                 gmx_bool bPrintRefValue,
  86                                 gmx_bool bPrintComponents)
  87 {
  88     fprintf(out, "\t%g", pcrd->value);
  89
  90     if (bPrintRefValue)
  91     {
  92         fprintf(out, "\t%g", pcrd->value_ref);
  93     }
  94
  95     if (bPrintComponents)
  96     {
  97         int m;
  98
  99         for (m = 0; m < DIM; m++)
 100         {
 101             if (pcrd->params.dim[m])
 102             {
 103                 fprintf(out, "\t%g", pcrd->dr[m]);
 104             }
 105         }
 106     }
 107 }
 108
 109 static void pull_print_x(FILE *out, struct pull_t *pull, double t)
 110 {
 111     int c;
 112
 113     fprintf(out, "%.4f", t);
 114
 115     for (c = 0; c < pull->ncoord; c++)
 116     {
 117         pull_coord_work_t *pcrd;
 118
 119         pcrd = &pull->coord[c];
 120
 121         if (pull->params.bPrintCOM1)
 122         {
 123             if (pcrd->params.eGeom == epullgCYL)
 124             {
 125                 pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
 126             }
 127             else
 128             {
 129                 pull_print_group_x(out, pcrd->params.dim,
 130                                    &pull->group[pcrd->params.group[0]]);
 131             }
 132         }
 133         if (pull->params.bPrintCOM2)
 134         {
 135             pull_print_group_x(out, pcrd->params.dim,
 136                                &pull->group[pcrd->params.group[1]]);
 137         }
 138         pull_print_coord_dr(out, pcrd,
 139                             pull->params.bPrintRefValue,
 140                             pull->params.bPrintComp);
 141     }
 142     fprintf(out, "\n");
 143 }
 144
 145 static void pull_print_f(FILE *out, struct pull_t *pull, double t)
 146 {
 147     int c;
 148
 149     fprintf(out, "%.4f", t);
 150
 151     for (c = 0; c < pull->ncoord; c++)
 152     {
 153         fprintf(out, "\t%g", pull->coord[c].f_scal);
 154     }
 155     fprintf(out, "\n");
 156 }
 157
 158 void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time)
 159 {
 160     if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
 161     {
 162         pull_print_x(pull->out_x, pull, time);
 163     }
 164
 165     if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
 166     {
 167         pull_print_f(pull->out_f, pull, time);
 168     }
 169 }
 170
 171 static FILE *open_pull_out(const char *fn, struct pull_t *pull, const output_env_t oenv,
 172                            gmx_bool bCoord, unsigned long Flags)
 173 {
 174     FILE  *fp;
 175     int    nsets, c, m;
 176     char **setname, buf[10];
 177
 178     if (Flags & MD_APPENDFILES)
 179     {
 180         fp = gmx_fio_fopen(fn, "a+");
 181     }
 182     else
 183     {
 184         fp = gmx_fio_fopen(fn, "w+");
 185         if (bCoord)
 186         {
 187             xvgr_header(fp, "Pull COM",  "Time (ps)", "Position (nm)",
 188                         exvggtXNY, oenv);
 189         }
 190         else
 191         {
 192             xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
 193                         exvggtXNY, oenv);
 194         }
 195
 196         /* With default mdp options only the actual distance is printed,
 197          * but optionally 2 COMs, the reference distance and distance
 198          * components can also be printed.
 199          */
 200         snew(setname, pull->ncoord*(DIM + DIM + 1 + 1 + DIM));
 201         nsets = 0;
 202         for (c = 0; c < pull->ncoord; c++)
 203         {
 204             if (bCoord)
 205             {
 206                 /* The order of this legend should match the order of printing
 207                  * the data in print_pull_x above.
 208                  */
 209
 210                 if (pull->params.bPrintCOM1)
 211                 {
 212                     /* Legend for reference group position */
 213                     for (m = 0; m < DIM; m++)
 214                     {
 215                         if (pull->coord[c].params.dim[m])
 216                         {
 217                             sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
 218                             setname[nsets] = gmx_strdup(buf);
 219                             nsets++;
 220                         }
 221                     }
 222                 }
 223                 if (pull->params.bPrintCOM2)
 224                 {
 225                     /* Legend for reference group position */
 226                     for (m = 0; m < DIM; m++)
 227                     {
 228                         if (pull->coord[c].params.dim[m])
 229                         {
 230                             sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
 231                             setname[nsets] = gmx_strdup(buf);
 232                             nsets++;
 233                         }
 234                     }
 235                 }
 236                 /* The pull coord distance */
 237                 sprintf(buf, "%d", c+1);
 238                 setname[nsets] = gmx_strdup(buf);
 239                 nsets++;
 240                 if (pull->params.bPrintRefValue)
 241                 {
 242                     sprintf(buf, "%c%d", 'r', c+1);
 243                     setname[nsets] = gmx_strdup(buf);
 244                     nsets++;
 245                 }
 246                 if (pull->params.bPrintComp)
 247                 {
 248                     /* The pull coord distance components */
 249                     for (m = 0; m < DIM; m++)
 250                     {
 251                         if (pull->coord[c].params.dim[m])
 252                         {
 253                             sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
 254                             setname[nsets] = gmx_strdup(buf);
 255                             nsets++;
 256                         }
 257                     }
 258                 }
 259             }
 260             else
 261             {
 262                 /* For the pull force we always only use one scalar */
 263                 sprintf(buf, "%d", c+1);
 264                 setname[nsets] = gmx_strdup(buf);
 265                 nsets++;
 266             }
 267         }
 268         if (nsets > 1)
 269         {
 270             xvgr_legend(fp, nsets, (const char**)setname, oenv);
 271         }
 272         for (c = 0; c < nsets; c++)
 273         {
 274             sfree(setname[c]);
 275         }
 276         sfree(setname);
 277     }
 278
 279     return fp;
 280 }
 281
 282 /* Apply forces in a mass weighted fashion */
 283 static void apply_forces_grp(const pull_group_work_t *pgrp,
 284                              const t_mdatoms *md,
 285                              const dvec f_pull, int sign, rvec *f)
 286 {
 287     int    i, ii, m;
 288     double wmass, inv_wm;
 289
 290     inv_wm = pgrp->mwscale;
 291
 292     if (pgrp->params.nat == 1 && pgrp->nat_loc == 1)
 293     {
 294         /* Only one atom and our rank has this atom: we can skip
 295          * the mass weighting, which means that this code also works
 296          * for mass=0, e.g. with a virtual site.
 297          */
 298         for (m = 0; m < DIM; m++)
 299         {
 300             f[pgrp->ind_loc[0]][m] += sign*f_pull[m];
 301         }
 302     }
 303     else
 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
 322 /* Apply forces in a mass weighted fashion to a cylinder group */
 323 static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
 324                                  double dv_corr,
 325                                  const t_mdatoms *md,
 326                                  const dvec f_pull, double f_scal,
 327                                  int sign, rvec *f)
 328 {
 329     int    i, ii, m;
 330     double mass, weight, inv_wm, dv_com;
 331
 332     inv_wm = pgrp->mwscale;
 333
 334     for (i = 0; i < pgrp->nat_loc; i++)
 335     {
 336         ii     = pgrp->ind_loc[i];
 337         mass   = md->massT[ii];
 338         weight = pgrp->weight_loc[i];
 339         /* The stored axial distance from the cylinder center (dv) needs
 340          * to be corrected for an offset (dv_corr), which was unknown when
 341          * we calculated dv.
 342          */
 343         dv_com = pgrp->dv[i] + dv_corr;
 344
 345         /* Here we not only add the pull force working along vec (f_pull),
 346          * but also a radial component, due to the dependence of the weights
 347          * on the radial distance.
 348          */
 349         for (m = 0; m < DIM; m++)
 350         {
 351             f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
 352                                      pgrp->mdw[i][m]*dv_com*f_scal);
 353         }
 354     }
 355 }
 356
 357 /* Apply torque forces in a mass weighted fashion to the groups that define
 358  * the pull vector direction for pull coordinate pcrd.
 359  */
 360 static void apply_forces_vec_torque(const struct pull_t     *pull,
 361                                     const pull_coord_work_t *pcrd,
 362                                     const t_mdatoms         *md,
 363                                     rvec                    *f)
 364 {
 365     double inpr;
 366     int    m;
 367     dvec   f_perp;
 368
 369     /* The component inpr along the pull vector is accounted for in the usual
 370      * way. Here we account for the component perpendicular to vec.
 371      */
 372     inpr = 0;
 373     for (m = 0; m < DIM; m++)
 374     {
 375         inpr += pcrd->dr[m]*pcrd->vec[m];
 376     }
 377     /* The torque force works along the component of the distance vector
 378      * of between the two "usual" pull groups that is perpendicular to
 379      * the pull vector. The magnitude of this force is the "usual" scale force
 380      * multiplied with the ratio of the distance between the two "usual" pull
 381      * groups and the distance between the two groups that define the vector.
 382      */
 383     for (m = 0; m < DIM; m++)
 384     {
 385         f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal;
 386     }
 387
 388     /* Apply the force to the groups defining the vector using opposite signs */
 389     apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f);
 390     apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp,  1, f);
 391 }
 392
 393 /* Apply forces in a mass weighted fashion */
 394 static void apply_forces(struct pull_t * pull, t_mdatoms * md, rvec *f)
 395 {
 396     int c;
 397
 398     for (c = 0; c < pull->ncoord; c++)
 399     {
 400         const pull_coord_work_t *pcrd;
 401
 402         pcrd = &pull->coord[c];
 403
 404         if (pcrd->params.eGeom == epullgCYL)
 405         {
 406             dvec f_tot;
 407             int  m;
 408
 409             apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
 410                                  pcrd->f, pcrd->f_scal, -1, f);
 411
 412             /* Sum the force along the vector and the radial force */
 413             for (m = 0; m < DIM; m++)
 414             {
 415                 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
 416             }
 417             apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f);
 418         }
 419         else
 420         {
 421             if (pcrd->params.eGeom == epullgDIRRELATIVE)
 422             {
 423                 /* We need to apply the torque forces to the pull groups
 424                  * that define the pull vector.
 425                  */
 426                 apply_forces_vec_torque(pull, pcrd, md, f);
 427             }
 428
 429             if (pull->group[pcrd->params.group[0]].params.nat > 0)
 430             {
 431                 apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f);
 432             }
 433             apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f);
 434         }
 435     }
 436 }
 437
 438 static double max_pull_distance2(const pull_coord_work_t *pcrd,
 439                                  const t_pbc             *pbc)
 440 {
 441     double max_d2;
 442     int    m;
 443
 444     max_d2 = GMX_DOUBLE_MAX;
 445
 446     for (m = 0; m < pbc->ndim_ePBC; m++)
 447     {
 448         if (pcrd->params.dim[m] != 0)
 449         {
 450             max_d2 = std::min(max_d2, static_cast<double>(norm2(pbc->box[m])));
 451         }
 452     }
 453
 454     return 0.25*max_d2;
 455 }
 456
 457 /* This function returns the distance based on coordinates xg and xref.
 458  * Note that the pull coordinate struct pcrd is not modified.
 459  */
 460 static void low_get_pull_coord_dr(const struct pull_t *pull,
 461                                   const pull_coord_work_t *pcrd,
 462                                   const t_pbc *pbc,
 463                                   dvec xg, dvec xref, double max_dist2,
 464                                   dvec dr)
 465 {
 466     const pull_group_work_t *pgrp0;
 467     int                      m;
 468     dvec                     xrefr, dref = {0, 0, 0};
 469     double                   dr2;
 470
 471     pgrp0 = &pull->group[pcrd->params.group[0]];
 472
 473     /* Only the first group can be an absolute reference, in that case nat=0 */
 474     if (pgrp0->params.nat == 0)
 475     {
 476         for (m = 0; m < DIM; m++)
 477         {
 478             xref[m] = pcrd->params.origin[m];
 479         }
 480     }
 481
 482     copy_dvec(xref, xrefr);
 483
 484     if (pcrd->params.eGeom == epullgDIRPBC)
 485     {
 486         for (m = 0; m < DIM; m++)
 487         {
 488             dref[m] = pcrd->value_ref*pcrd->vec[m];
 489         }
 490         /* Add the reference position, so we use the correct periodic image */
 491         dvec_inc(xrefr, dref);
 492     }
 493
 494     pbc_dx_d(pbc, xg, xrefr, dr);
 495     dr2 = 0;
 496     for (m = 0; m < DIM; m++)
 497     {
 498         dr[m] *= pcrd->params.dim[m];
 499         dr2   += dr[m]*dr[m];
 500     }
 501     if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
 502     {
 503         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",
 504                   pcrd->params.group[0], pcrd->params.group[1],
 505                   sqrt(dr2), sqrt(max_dist2));
 506     }
 507
 508     if (pcrd->params.eGeom == epullgDIRPBC)
 509     {
 510         dvec_inc(dr, dref);
 511     }
 512 }
 513
 514 /* This function returns the distance based on the contents of the pull struct.
 515  * pull->coord[coord_ind].dr, and potentially vector, are updated.
 516  */
 517 static void get_pull_coord_dr(struct pull_t *pull,
 518                               int            coord_ind,
 519                               const t_pbc   *pbc)
 520 {
 521     double             md2;
 522     pull_coord_work_t *pcrd;
 523
 524     pcrd = &pull->coord[coord_ind];
 525
 526     if (pcrd->params.eGeom == epullgDIRPBC)
 527     {
 528         md2 = -1;
 529     }
 530     else
 531     {
 532         md2 = max_pull_distance2(pcrd, pbc);
 533     }
 534
 535     if (pcrd->params.eGeom == epullgDIRRELATIVE)
 536     {
 537         /* We need to determine the pull vector */
 538         const pull_group_work_t *pgrp2, *pgrp3;
 539         dvec                     vec;
 540         int                      m;
 541
 542         pgrp2 = &pull->group[pcrd->params.group[2]];
 543         pgrp3 = &pull->group[pcrd->params.group[3]];
 544
 545         pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
 546
 547         for (m = 0; m < DIM; m++)
 548         {
 549             vec[m] *= pcrd->params.dim[m];
 550         }
 551         pcrd->vec_len = dnorm(vec);
 552         for (m = 0; m < DIM; m++)
 553         {
 554             pcrd->vec[m] = vec[m]/pcrd->vec_len;
 555         }
 556         if (debug)
 557         {
 558             fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
 559                     coord_ind,
 560                     vec[XX], vec[YY], vec[ZZ],
 561                     pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]);
 562         }
 563     }
 564
 565     low_get_pull_coord_dr(pull, pcrd, pbc,
 566                           pull->group[pcrd->params.group[1]].x,
 567                           pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->params.group[0]].x,
 568                           md2,
 569                           pcrd->dr);
 570 }
 571
 572 static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, double t)
 573 {
 574     /* With zero rate the reference value is set initially and doesn't change */
 575     if (pcrd->params.rate != 0)
 576     {
 577         pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
 578     }
 579 }
 580
 581 /* Calculates pull->coord[coord_ind].value.
 582  * This function also updates pull->coord[coord_ind].dr.
 583  */
 584 static void get_pull_coord_distance(struct pull_t *pull,
 585                                     int            coord_ind,
 586                                     const t_pbc   *pbc)
 587 {
 588     pull_coord_work_t *pcrd;
 589     int                m;
 590
 591     pcrd = &pull->coord[coord_ind];
 592
 593     get_pull_coord_dr(pull, coord_ind, pbc);
 594
 595     switch (pcrd->params.eGeom)
 596     {
 597         case epullgDIST:
 598             /* Pull along the vector between the com's */
 599             pcrd->value = dnorm(pcrd->dr);
 600             break;
 601         case epullgDIR:
 602         case epullgDIRPBC:
 603         case epullgDIRRELATIVE:
 604         case epullgCYL:
 605             /* Pull along vec */
 606             pcrd->value = 0;
 607             for (m = 0; m < DIM; m++)
 608             {
 609                 pcrd->value += pcrd->vec[m]*pcrd->dr[m];
 610             }
 611             break;
 612     }
 613 }
 614
 615 /* Returns the deviation from the reference value.
 616  * Updates pull->coord[coord_ind].dr, .value and .value_ref.
 617  */
 618 static double get_pull_coord_deviation(struct pull_t *pull,
 619                                        int            coord_ind,
 620                                        const t_pbc   *pbc,
 621                                        double         t)
 622 {
 623     static gmx_bool    bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
 624                                            but is fairly benign */
 625     pull_coord_work_t *pcrd;
 626     double             dev = 0;
 627
 628     pcrd = &pull->coord[coord_ind];
 629
 630     get_pull_coord_distance(pull, coord_ind, pbc);
 631
 632     update_pull_coord_reference_value(pcrd, t);
 633
 634     switch (pcrd->params.eGeom)
 635     {
 636         case epullgDIST:
 637             /* Pull along the vector between the com's */
 638             if (pcrd->value_ref < 0 && !bWarned)
 639             {
 640                 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
 641                 bWarned = TRUE;
 642             }
 643             if (pcrd->value == 0)
 644             {
 645                 /* With no vector we can not determine the direction for the force,
 646                  * so we set the force to zero.
 647                  */
 648                 dev = 0;
 649             }
 650             else
 651             {
 652                 /* Determine the deviation */
 653                 dev = pcrd->value - pcrd->value_ref;
 654             }
 655             break;
 656         case epullgDIR:
 657         case epullgDIRPBC:
 658         case epullgDIRRELATIVE:
 659         case epullgCYL:
 660             /* Pull along vec */
 661             dev = pcrd->value - pcrd->value_ref;
 662             break;
 663     }
 664
 665     return dev;
 666 }
 667
 668 void get_pull_coord_value(struct pull_t *pull,
 669                           int            coord_ind,
 670                           const t_pbc   *pbc,
 671                           double        *value)
 672 {
 673     get_pull_coord_distance(pull, coord_ind, pbc);
 674
 675     *value = pull->coord[coord_ind].value;
 676 }
 677
 678 void clear_pull_forces(struct pull_t *pull)
 679 {
 680     int i;
 681
 682     /* Zeroing the forces is only required for constraint pulling.
 683      * It can happen that multiple constraint steps need to be applied
 684      * and therefore the constraint forces need to be accumulated.
 685      */
 686     for (i = 0; i < pull->ncoord; i++)
 687     {
 688         clear_dvec(pull->coord[i].f);
 689         pull->coord[i].f_scal = 0;
 690     }
 691 }
 692
 693 /* Apply constraint using SHAKE */
 694 static void do_constraint(struct pull_t *pull, t_pbc *pbc,
 695                           rvec *x, rvec *v,
 696                           gmx_bool bMaster, tensor vir,
 697                           double dt, double t)
 698 {
 699
 700     dvec         *r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
 701     dvec          unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
 702     dvec         *rnew;   /* current 'new' positions of the groups */
 703     double       *dr_tot; /* the total update of the coords */
 704     dvec          vec;
 705     double        inpr;
 706     double        lambda, rm, invdt = 0;
 707     gmx_bool      bConverged_all, bConverged = FALSE;
 708     int           niter = 0, g, c, ii, j, m, max_iter = 100;
 709     double        a;
 710     dvec          tmp, tmp3;
 711
 712     snew(r_ij,   pull->ncoord);
 713     snew(dr_tot, pull->ncoord);
 714
 715     snew(rnew, pull->ngroup);
 716
 717     /* copy the current unconstrained positions for use in iterations. We
 718        iterate until rinew[i] and rjnew[j] obey the constraints. Then
 719        rinew - pull.x_unc[i] is the correction dr to group i */
 720     for (g = 0; g < pull->ngroup; g++)
 721     {
 722         copy_dvec(pull->group[g].xp, rnew[g]);
 723     }
 724
 725     /* Determine the constraint directions from the old positions */
 726     for (c = 0; c < pull->ncoord; c++)
 727     {
 728         pull_coord_work_t *pcrd;
 729
 730         pcrd = &pull->coord[c];
 731
 732         if (pcrd->params.eType != epullCONSTRAINT)
 733         {
 734             continue;
 735         }
 736
 737         /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
 738          * We don't modify dr and value anymore, so these values are also used
 739          * for printing.
 740          */
 741         get_pull_coord_distance(pull, c, pbc);
 742
 743         if (debug)
 744         {
 745             fprintf(debug, "Pull coord %d dr %f %f %f\n",
 746                     c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
 747         }
 748
 749         if (pcrd->params.eGeom == epullgDIR ||
 750             pcrd->params.eGeom == epullgDIRPBC)
 751         {
 752             /* Select the component along vec */
 753             a = 0;
 754             for (m = 0; m < DIM; m++)
 755             {
 756                 a += pcrd->vec[m]*pcrd->dr[m];
 757             }
 758             for (m = 0; m < DIM; m++)
 759             {
 760                 r_ij[c][m] = a*pcrd->vec[m];
 761             }
 762         }
 763         else
 764         {
 765             copy_dvec(pcrd->dr, r_ij[c]);
 766         }
 767
 768         if (dnorm2(r_ij[c]) == 0)
 769         {
 770             gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
 771         }
 772     }
 773
 774     bConverged_all = FALSE;
 775     while (!bConverged_all && niter < max_iter)
 776     {
 777         bConverged_all = TRUE;
 778
 779         /* loop over all constraints */
 780         for (c = 0; c < pull->ncoord; c++)
 781         {
 782             pull_coord_work_t *pcrd;
 783             pull_group_work_t *pgrp0, *pgrp1;
 784             dvec               dr0, dr1;
 785
 786             pcrd  = &pull->coord[c];
 787
 788             if (pcrd->params.eType != epullCONSTRAINT)
 789             {
 790                 continue;
 791             }
 792
 793             update_pull_coord_reference_value(pcrd, t);
 794
 795             pgrp0 = &pull->group[pcrd->params.group[0]];
 796             pgrp1 = &pull->group[pcrd->params.group[1]];
 797
 798             /* Get the current difference vector */
 799             low_get_pull_coord_dr(pull, pcrd, pbc,
 800                                   rnew[pcrd->params.group[1]],
 801                                   rnew[pcrd->params.group[0]],
 802                                   -1, unc_ij);
 803
 804             if (debug)
 805             {
 806                 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
 807             }
 808
 809             rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
 810
 811             switch (pcrd->params.eGeom)
 812             {
 813                 case epullgDIST:
 814                     if (pcrd->value_ref <= 0)
 815                     {
 816                         gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
 817                     }
 818
 819                     {
 820                         double q, c_a, c_b, c_c;
 821
 822                         c_a = diprod(r_ij[c], r_ij[c]);
 823                         c_b = diprod(unc_ij, r_ij[c])*2;
 824                         c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
 825
 826                         if (c_b < 0)
 827                         {
 828                             q      = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
 829                             lambda = -q/c_a;
 830                         }
 831                         else
 832                         {
 833                             q      = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
 834                             lambda = -c_c/q;
 835                         }
 836
 837                         if (debug)
 838                         {
 839                             fprintf(debug,
 840                                     "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
 841                                     c_a, c_b, c_c, lambda);
 842                         }
 843                     }
 844
 845                     /* The position corrections dr due to the constraints */
 846                     dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
 847                     dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
 848                     dr_tot[c] += -lambda*dnorm(r_ij[c]);
 849                     break;
 850                 case epullgDIR:
 851                 case epullgDIRPBC:
 852                 case epullgCYL:
 853                     /* A 1-dimensional constraint along a vector */
 854                     a = 0;
 855                     for (m = 0; m < DIM; m++)
 856                     {
 857                         vec[m] = pcrd->vec[m];
 858                         a     += unc_ij[m]*vec[m];
 859                     }
 860                     /* Select only the component along the vector */
 861                     dsvmul(a, vec, unc_ij);
 862                     lambda = a - pcrd->value_ref;
 863                     if (debug)
 864                     {
 865                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
 866                     }
 867
 868                     /* The position corrections dr due to the constraints */
 869                     dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
 870                     dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
 871                     dr_tot[c] += -lambda;
 872                     break;
 873                 default:
 874                     gmx_incons("Invalid enumeration value for eGeom");
 875                     /* Keep static analyzer happy */
 876                     clear_dvec(dr0);
 877                     clear_dvec(dr1);
 878             }
 879
 880             /* DEBUG */
 881             if (debug)
 882             {
 883                 int g0, g1;
 884
 885                 g0 = pcrd->params.group[0];
 886                 g1 = pcrd->params.group[1];
 887                 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
 888                 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
 889                 fprintf(debug,
 890                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 891                         rnew[g0][0], rnew[g0][1], rnew[g0][2],
 892                         rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
 893                 fprintf(debug,
 894                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
 895                         "", "", "", "", "", "", pcrd->value_ref);
 896                 fprintf(debug,
 897                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 898                         dr0[0], dr0[1], dr0[2],
 899                         dr1[0], dr1[1], dr1[2],
 900                         dnorm(tmp3));
 901             } /* END DEBUG */
 902
 903             /* Update the COMs with dr */
 904             dvec_inc(rnew[pcrd->params.group[1]], dr1);
 905             dvec_inc(rnew[pcrd->params.group[0]], dr0);
 906         }
 907
 908         /* Check if all constraints are fullfilled now */
 909         for (c = 0; c < pull->ncoord; c++)
 910         {
 911             pull_coord_work_t *pcrd;
 912
 913             pcrd = &pull->coord[c];
 914
 915             if (pcrd->params.eType != epullCONSTRAINT)
 916             {
 917                 continue;
 918             }
 919
 920             low_get_pull_coord_dr(pull, pcrd, pbc,
 921                                   rnew[pcrd->params.group[1]],
 922                                   rnew[pcrd->params.group[0]],
 923                                   -1, unc_ij);
 924
 925             switch (pcrd->params.eGeom)
 926             {
 927                 case epullgDIST:
 928                     bConverged =
 929                         fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->params.constr_tol;
 930                     break;
 931                 case epullgDIR:
 932                 case epullgDIRPBC:
 933                 case epullgCYL:
 934                     for (m = 0; m < DIM; m++)
 935                     {
 936                         vec[m] = pcrd->vec[m];
 937                     }
 938                     inpr = diprod(unc_ij, vec);
 939                     dsvmul(inpr, vec, unc_ij);
 940                     bConverged =
 941                         fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->params.constr_tol;
 942                     break;
 943             }
 944
 945             if (!bConverged)
 946             {
 947                 if (debug)
 948                 {
 949                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
 950                             "d_ref = %f, current d = %f\n",
 951                             g, pcrd->value_ref, dnorm(unc_ij));
 952                 }
 953
 954                 bConverged_all = FALSE;
 955             }
 956         }
 957
 958         niter++;
 959         /* if after all constraints are dealt with and bConverged is still TRUE
 960            we're finished, if not we do another iteration */
 961     }
 962     if (niter > max_iter)
 963     {
 964         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
 965     }
 966
 967     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 968
 969     if (v)
 970     {
 971         invdt = 1/dt;
 972     }
 973
 974     /* update atoms in the groups */
 975     for (g = 0; g < pull->ngroup; g++)
 976     {
 977         const pull_group_work_t *pgrp;
 978         dvec                     dr;
 979
 980         pgrp = &pull->group[g];
 981
 982         /* get the final constraint displacement dr for group g */
 983         dvec_sub(rnew[g], pgrp->xp, dr);
 984
 985         if (dnorm2(dr) == 0)
 986         {
 987             /* No displacement, no update necessary */
 988             continue;
 989         }
 990
 991         /* update the atom positions */
 992         copy_dvec(dr, tmp);
 993         for (j = 0; j < pgrp->nat_loc; j++)
 994         {
 995             ii = pgrp->ind_loc[j];
 996             if (pgrp->weight_loc)
 997             {
 998                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
 999             }
1000             for (m = 0; m < DIM; m++)
1001             {
1002                 x[ii][m] += tmp[m];
1003             }
1004             if (v)
1005             {
1006                 for (m = 0; m < DIM; m++)
1007                 {
1008                     v[ii][m] += invdt*tmp[m];
1009                 }
1010             }
1011         }
1012     }
1013
1014     /* calculate the constraint forces, used for output and virial only */
1015     for (c = 0; c < pull->ncoord; c++)
1016     {
1017         pull_coord_work_t *pcrd;
1018
1019         pcrd = &pull->coord[c];
1020
1021         if (pcrd->params.eType != epullCONSTRAINT)
1022         {
1023             continue;
1024         }
1025
1026         pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
1027
1028         if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1029         {
1030             double f_invr;
1031
1032             /* Add the pull contribution to the virial */
1033             /* We have already checked above that r_ij[c] != 0 */
1034             f_invr = pcrd->f_scal/dnorm(r_ij[c]);
1035
1036             for (j = 0; j < DIM; j++)
1037             {
1038                 for (m = 0; m < DIM; m++)
1039                 {
1040                     vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1041                 }
1042             }
1043         }
1044     }
1045
1046     /* finished! I hope. Give back some memory */
1047     sfree(r_ij);
1048     sfree(dr_tot);
1049     sfree(rnew);
1050 }
1051
1052 /* Pulling with a harmonic umbrella potential or constant force */
1053 static void do_pull_pot(struct pull_t *pull, t_pbc *pbc, double t, real lambda,
1054                         real *V, tensor vir, real *dVdl)
1055 {
1056     int    c, j, m;
1057     double dev, ndr, invdr = 0;
1058     real   k, dkdl;
1059
1060     /* loop over the pull coordinates */
1061     *V    = 0;
1062     *dVdl = 0;
1063     for (c = 0; c < pull->ncoord; c++)
1064     {
1065         pull_coord_work_t *pcrd;
1066
1067         pcrd = &pull->coord[c];
1068
1069         if (pcrd->params.eType == epullCONSTRAINT)
1070         {
1071             continue;
1072         }
1073
1074         dev = get_pull_coord_deviation(pull, c, pbc, t);
1075
1076         k    = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
1077         dkdl = pcrd->params.kB - pcrd->params.k;
1078
1079         if (pcrd->params.eGeom == epullgDIST)
1080         {
1081             ndr   = dnorm(pcrd->dr);
1082             if (ndr > 0)
1083             {
1084                 invdr = 1/ndr;
1085             }
1086             else
1087             {
1088                 /* With an harmonic umbrella, the force is 0 at r=0,
1089                  * so we can set invdr to any value.
1090                  * With a constant force, the force at r=0 is not defined,
1091                  * so we zero it (this is anyhow a very rare event).
1092                  */
1093                 invdr = 0;
1094             }
1095         }
1096         else
1097         {
1098             ndr = 0;
1099             for (m = 0; m < DIM; m++)
1100             {
1101                 ndr += pcrd->vec[m]*pcrd->dr[m];
1102             }
1103         }
1104
1105         switch (pcrd->params.eType)
1106         {
1107             case epullUMBRELLA:
1108             case epullFLATBOTTOM:
1109                 /* The only difference between an umbrella and a flat-bottom
1110                  * potential is that a flat-bottom is zero below zero.
1111                  */
1112                 if (pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1113                 {
1114                     dev = 0;
1115                 }
1116
1117                 pcrd->f_scal  =       -k*dev;
1118                 *V           += 0.5*   k*dsqr(dev);
1119                 *dVdl        += 0.5*dkdl*dsqr(dev);
1120                 break;
1121             case epullCONST_F:
1122                 pcrd->f_scal  =   -k;
1123                 *V           +=    k*ndr;
1124                 *dVdl        += dkdl*ndr;
1125                 break;
1126             default:
1127                 gmx_incons("Unsupported pull type in do_pull_pot");
1128         }
1129
1130         if (pcrd->params.eGeom == epullgDIST)
1131         {
1132             for (m = 0; m < DIM; m++)
1133             {
1134                 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
1135             }
1136         }
1137         else
1138         {
1139             for (m = 0; m < DIM; m++)
1140             {
1141                 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
1142             }
1143         }
1144
1145         if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC)
1146         {
1147             /* Add the pull contribution to the virial */
1148             for (j = 0; j < DIM; j++)
1149             {
1150                 for (m = 0; m < DIM; m++)
1151                 {
1152                     vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
1153                 }
1154             }
1155         }
1156     }
1157 }
1158
1159 real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
1160                     t_commrec *cr, double t, real lambda,
1161                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
1162 {
1163     real V = 0, dVdl;
1164
1165     assert(pull != NULL);
1166
1167     if (pull->comm.bParticipate)
1168     {
1169         pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1170
1171         do_pull_pot(pull, pbc, t, lambda,
1172                     &V, MASTER(cr) ? vir : NULL, &dVdl);
1173
1174         /* Distribute forces over the pull groups */
1175         apply_forces(pull, md, f);
1176
1177         if (MASTER(cr))
1178         {
1179             *dvdlambda += dVdl;
1180         }
1181     }
1182
1183     return (MASTER(cr) ? V : 0.0);
1184 }
1185
1186 void pull_constraint(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
1187                      t_commrec *cr, double dt, double t,
1188                      rvec *x, rvec *xp, rvec *v, tensor vir)
1189 {
1190     assert(pull != NULL);
1191
1192     if (pull->comm.bParticipate)
1193     {
1194         pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1195
1196         do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1197     }
1198 }
1199
1200 static void make_local_pull_group(gmx_ga2la_t ga2la,
1201                                   pull_group_work_t *pg, int start, int end)
1202 {
1203     int i, ii;
1204
1205     pg->nat_loc = 0;
1206     for (i = 0; i < pg->params.nat; i++)
1207     {
1208         ii = pg->params.ind[i];
1209         if (ga2la)
1210         {
1211             if (!ga2la_get_home(ga2la, ii, &ii))
1212             {
1213                 ii = -1;
1214             }
1215         }
1216         if (ii >= start && ii < end)
1217         {
1218             /* This is a home atom, add it to the local pull group */
1219             if (pg->nat_loc >= pg->nalloc_loc)
1220             {
1221                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1222                 srenew(pg->ind_loc, pg->nalloc_loc);
1223                 if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != NULL)
1224                 {
1225                     srenew(pg->weight_loc, pg->nalloc_loc);
1226                 }
1227             }
1228             pg->ind_loc[pg->nat_loc] = ii;
1229             if (pg->params.weight != NULL)
1230             {
1231                 pg->weight_loc[pg->nat_loc] = pg->params.weight[i];
1232             }
1233             pg->nat_loc++;
1234         }
1235     }
1236 }
1237
1238 void dd_make_local_pull_groups(t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
1239 {
1240     gmx_domdec_t *dd;
1241     pull_comm_t  *comm;
1242     gmx_ga2la_t   ga2la;
1243     gmx_bool      bMustParticipate;
1244     int           g;
1245
1246     dd = cr->dd;
1247
1248     comm = &pull->comm;
1249
1250     if (dd)
1251     {
1252         ga2la = dd->ga2la;
1253     }
1254     else
1255     {
1256         ga2la = NULL;
1257     }
1258
1259     /* We always make the master node participate, such that it can do i/o
1260      * and to simplify MC type extensions people might have.
1261      */
1262     bMustParticipate = (comm->bParticipateAll || dd == NULL || DDMASTER(dd));
1263
1264     for (g = 0; g < pull->ngroup; g++)
1265     {
1266         int a;
1267
1268         make_local_pull_group(ga2la, &pull->group[g],
1269                               0, md->homenr);
1270
1271         /* We should participate if we have pull or pbc atoms */
1272         if (!bMustParticipate &&
1273             (pull->group[g].nat_loc > 0 ||
1274              (pull->group[g].params.pbcatom >= 0 &&
1275               ga2la_get_home(dd->ga2la, pull->group[g].params.pbcatom, &a))))
1276         {
1277             bMustParticipate = TRUE;
1278         }
1279     }
1280
1281     if (!comm->bParticipateAll)
1282     {
1283         /* Keep currently not required ranks in the communicator
1284          * if they needed to participate up to 20 decompositions ago.
1285          * This avoids frequent rebuilds due to atoms jumping back and forth.
1286          */
1287         const gmx_int64_t history_count = 20;
1288         gmx_bool          bWillParticipate;
1289         int               count[2];
1290
1291         /* Increase the decomposition counter for the current call */
1292         comm->setup_count++;
1293
1294         if (bMustParticipate)
1295         {
1296             comm->must_count = comm->setup_count;
1297         }
1298
1299         bWillParticipate =
1300             bMustParticipate ||
1301             (comm->bParticipate &&
1302              comm->must_count >= comm->setup_count - history_count);
1303
1304         if (debug && dd != NULL)
1305         {
1306             fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %d, will be part %d\n",
1307                     dd->rank, bMustParticipate, bWillParticipate);
1308         }
1309
1310         if (bWillParticipate)
1311         {
1312             /* Count the number of ranks that we want to have participating */
1313             count[0] = 1;
1314             /* Count the number of ranks that need to be added */
1315             count[1] = comm->bParticipate ? 0 : 1;
1316         }
1317         else
1318         {
1319             count[0] = 0;
1320             count[1] = 0;
1321         }
1322
1323         /* The cost of this global operation will be less that the cost
1324          * of the extra MPI_Comm_split calls that we can avoid.
1325          */
1326         gmx_sumi(2, count, cr);
1327
1328         /* If we are missing ranks or if we have 20% more ranks than needed
1329          * we make a new sub-communicator.
1330          */
1331         if (count[1] > 0 || 6*count[0] < 5*comm->nparticipate)
1332         {
1333             if (debug)
1334             {
1335                 fprintf(debug, "Creating new pull subcommunicator of size %d\n",
1336                         count[0]);
1337             }
1338
1339 #ifdef GMX_MPI
1340             if (comm->mpi_comm_com != MPI_COMM_NULL)
1341             {
1342                 MPI_Comm_free(&comm->mpi_comm_com);
1343             }
1344
1345             /* This might be an extremely expensive operation, so we try
1346              * to avoid this splitting as much as possible.
1347              */
1348             assert(dd != NULL);
1349             MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank,
1350                            &comm->mpi_comm_com);
1351 #endif
1352
1353             comm->bParticipate = bWillParticipate;
1354             comm->nparticipate = count[0];
1355         }
1356     }
1357
1358     /* Since the PBC of atoms might have changed, we need to update the PBC */
1359     pull->bSetPBCatoms = TRUE;
1360 }
1361
1362 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1363                                   int g, pull_group_work_t *pg,
1364                                   gmx_bool bConstraint, ivec pulldim_con,
1365                                   const gmx_mtop_t *mtop,
1366                                   const t_inputrec *ir, real lambda)
1367 {
1368     int                   i, ii, d, nfrozen, ndim;
1369     real                  m, w, mbd;
1370     double                tmass, wmass, wwmass;
1371     const gmx_groups_t   *groups;
1372     gmx_mtop_atomlookup_t alook;
1373     t_atom               *atom;
1374
1375     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1376     {
1377         /* There are no masses in the integrator.
1378          * But we still want to have the correct mass-weighted COMs.
1379          * So we store the real masses in the weights.
1380          * We do not set nweight, so these weights do not end up in the tpx file.
1381          */
1382         if (pg->params.nweight == 0)
1383         {
1384             snew(pg->params.weight, pg->params.nat);
1385         }
1386     }
1387
1388     if (cr && PAR(cr))
1389     {
1390         pg->nat_loc    = 0;
1391         pg->nalloc_loc = 0;
1392         pg->ind_loc    = NULL;
1393         pg->weight_loc = NULL;
1394     }
1395     else
1396     {
1397         pg->nat_loc = pg->params.nat;
1398         pg->ind_loc = pg->params.ind;
1399         if (pg->epgrppbc == epgrppbcCOS)
1400         {
1401             snew(pg->weight_loc, pg->params.nat);
1402         }
1403         else
1404         {
1405             pg->weight_loc = pg->params.weight;
1406         }
1407     }
1408
1409     groups = &mtop->groups;
1410
1411     alook = gmx_mtop_atomlookup_init(mtop);
1412
1413     nfrozen = 0;
1414     tmass   = 0;
1415     wmass   = 0;
1416     wwmass  = 0;
1417     for (i = 0; i < pg->params.nat; i++)
1418     {
1419         ii = pg->params.ind[i];
1420         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1421         if (bConstraint && ir->opts.nFreeze)
1422         {
1423             for (d = 0; d < DIM; d++)
1424             {
1425                 if (pulldim_con[d] == 1 &&
1426                     ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1427                 {
1428                     nfrozen++;
1429                 }
1430             }
1431         }
1432         if (ir->efep == efepNO)
1433         {
1434             m = atom->m;
1435         }
1436         else
1437         {
1438             m = (1 - lambda)*atom->m + lambda*atom->mB;
1439         }
1440         if (pg->params.nweight > 0)
1441         {
1442             w = pg->params.weight[i];
1443         }
1444         else
1445         {
1446             w = 1;
1447         }
1448         if (EI_ENERGY_MINIMIZATION(ir->eI))
1449         {
1450             /* Move the mass to the weight */
1451             w                   *= m;
1452             m                    = 1;
1453             pg->params.weight[i] = w;
1454         }
1455         else if (ir->eI == eiBD)
1456         {
1457             if (ir->bd_fric)
1458             {
1459                 mbd = ir->bd_fric*ir->delta_t;
1460             }
1461             else
1462             {
1463                 if (groups->grpnr[egcTC] == NULL)
1464                 {
1465                     mbd = ir->delta_t/ir->opts.tau_t[0];
1466                 }
1467                 else
1468                 {
1469                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1470                 }
1471             }
1472             w                   *= m/mbd;
1473             m                    = mbd;
1474             pg->params.weight[i] = w;
1475         }
1476         tmass  += m;
1477         wmass  += m*w;
1478         wwmass += m*w*w;
1479     }
1480
1481     gmx_mtop_atomlookup_destroy(alook);
1482
1483     if (wmass == 0)
1484     {
1485         /* We can have single atom groups with zero mass with potential pulling
1486          * without cosine weighting.
1487          */
1488         if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1489         {
1490             /* With one atom the mass doesn't matter */
1491             wwmass = 1;
1492         }
1493         else
1494         {
1495             gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1496                       pg->params.weight ? " weighted" : "", g);
1497         }
1498     }
1499     if (fplog)
1500     {
1501         fprintf(fplog,
1502                 "Pull group %d: %5d atoms, mass %9.3f",
1503                 g, pg->params.nat, tmass);
1504         if (pg->params.weight ||
1505             EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1506         {
1507             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1508         }
1509         if (pg->epgrppbc == epgrppbcCOS)
1510         {
1511             fprintf(fplog, ", cosine weighting will be used");
1512         }
1513         fprintf(fplog, "\n");
1514     }
1515
1516     if (nfrozen == 0)
1517     {
1518         /* A value != 0 signals not frozen, it is updated later */
1519         pg->invtm  = -1.0;
1520     }
1521     else
1522     {
1523         ndim = 0;
1524         for (d = 0; d < DIM; d++)
1525         {
1526             ndim += pulldim_con[d]*pg->params.nat;
1527         }
1528         if (fplog && nfrozen > 0 && nfrozen < ndim)
1529         {
1530             fprintf(fplog,
1531                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1532                     "         that are subject to pulling are frozen.\n"
1533                     "         For constraint pulling the whole group will be frozen.\n\n",
1534                     g);
1535         }
1536         pg->invtm  = 0.0;
1537         pg->wscale = 1.0;
1538     }
1539 }
1540
1541 struct pull_t *
1542 init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
1543           int nfile, const t_filenm fnm[],
1544           gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1545           gmx_bool bOutFile, unsigned long Flags)
1546 {
1547     struct pull_t *pull;
1548     pull_comm_t   *comm;
1549     int            g, c, m;
1550
1551     snew(pull, 1);
1552
1553     /* Copy the pull parameters */
1554     pull->params       = *pull_params;
1555     /* Avoid pointer copies */
1556     pull->params.group = NULL;
1557     pull->params.coord = NULL;
1558
1559     pull->ncoord       = pull_params->ncoord;
1560     pull->ngroup       = pull_params->ngroup;
1561     snew(pull->coord, pull->ncoord);
1562     snew(pull->group, pull->ngroup);
1563
1564     pull->bPotential  = FALSE;
1565     pull->bConstraint = FALSE;
1566     pull->bCylinder   = FALSE;
1567
1568     for (g = 0; g < pull->ngroup; g++)
1569     {
1570         pull_group_work_t *pgrp;
1571         int                i;
1572
1573         pgrp = &pull->group[g];
1574
1575         /* Copy the pull group parameters */
1576         pgrp->params = pull_params->group[g];
1577
1578         /* Avoid pointer copies by allocating and copying arrays */
1579         snew(pgrp->params.ind, pgrp->params.nat);
1580         for (i = 0; i < pgrp->params.nat; i++)
1581         {
1582             pgrp->params.ind[i] = pull_params->group[g].ind[i];
1583         }
1584         if (pgrp->params.nweight > 0)
1585         {
1586             snew(pgrp->params.ind, pgrp->params.nweight);
1587             for (i = 0; i < pgrp->params.nweight; i++)
1588             {
1589                 pgrp->params.weight[i] = pull_params->group[g].weight[i];
1590             }
1591         }
1592     }
1593
1594     for (c = 0; c < pull->ncoord; c++)
1595     {
1596         pull_coord_work_t *pcrd;
1597         int                calc_com_start, calc_com_end, g;
1598
1599         pcrd = &pull->coord[c];
1600
1601         /* Copy all pull coordinate parameters */
1602         pcrd->params = pull_params->coord[c];
1603
1604         switch (pcrd->params.eGeom)
1605         {
1606             case epullgDIST:
1607             case epullgDIRRELATIVE:
1608                 /* Direction vector is determined at each step */
1609                 break;
1610             case epullgDIR:
1611             case epullgDIRPBC:
1612             case epullgCYL:
1613                 copy_rvec(pull_params->coord[c].vec, pcrd->vec);
1614                 break;
1615             default:
1616                 gmx_incons("Pull geometry not handled");
1617         }
1618
1619         if (pcrd->params.eType == epullCONSTRAINT)
1620         {
1621             /* Check restrictions of the constraint pull code */
1622             if (pcrd->params.eGeom == epullgCYL ||
1623                 pcrd->params.eGeom == epullgDIRRELATIVE)
1624             {
1625                 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
1626                           epull_names[pcrd->params.eType],
1627                           epullg_names[pcrd->params.eGeom],
1628                           epull_names[epullUMBRELLA]);
1629             }
1630
1631             pull->bConstraint = TRUE;
1632         }
1633         else
1634         {
1635             pull->bPotential = TRUE;
1636         }
1637
1638         if (pcrd->params.eGeom == epullgCYL)
1639         {
1640             pull->bCylinder = TRUE;
1641         }
1642         /* We only need to calculate the plain COM of a group
1643          * when it is not only used as a cylinder group.
1644          */
1645         calc_com_start = (pcrd->params.eGeom == epullgCYL         ? 1 : 0);
1646         calc_com_end   = (pcrd->params.eGeom == epullgDIRRELATIVE ? 4 : 2);
1647
1648         for (g = calc_com_start; g <= calc_com_end; g++)
1649         {
1650             pull->group[pcrd->params.group[g]].bCalcCOM = TRUE;
1651         }
1652
1653         /* With non-zero rate the reference value is set at every step */
1654         if (pcrd->params.rate == 0)
1655         {
1656             /* Initialize the constant reference value */
1657             pcrd->value_ref = pcrd->params.init;
1658         }
1659     }
1660
1661     pull->ePBC = ir->ePBC;
1662     switch (pull->ePBC)
1663     {
1664         case epbcNONE: pull->npbcdim = 0; break;
1665         case epbcXY:   pull->npbcdim = 2; break;
1666         default:       pull->npbcdim = 3; break;
1667     }
1668
1669     if (fplog)
1670     {
1671         gmx_bool bAbs, bCos;
1672
1673         bAbs = FALSE;
1674         for (c = 0; c < pull->ncoord; c++)
1675         {
1676             if (pull->group[pull->coord[c].params.group[0]].params.nat == 0 ||
1677                 pull->group[pull->coord[c].params.group[1]].params.nat == 0)
1678             {
1679                 bAbs = TRUE;
1680             }
1681         }
1682
1683         fprintf(fplog, "\n");
1684         if (pull->bPotential)
1685         {
1686             fprintf(fplog, "Will apply potential COM pulling\n");
1687         }
1688         if (pull->bConstraint)
1689         {
1690             fprintf(fplog, "Will apply constraint COM pulling\n");
1691         }
1692         fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1693                 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1694                 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1695         if (bAbs)
1696         {
1697             fprintf(fplog, "with an absolute reference\n");
1698         }
1699         bCos = FALSE;
1700         for (g = 0; g < pull->ngroup; g++)
1701         {
1702             if (pull->group[g].params.nat > 1 &&
1703                 pull->group[g].params.pbcatom < 0)
1704             {
1705                 /* We are using cosine weighting */
1706                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1707                 bCos = TRUE;
1708             }
1709         }
1710         if (bCos)
1711         {
1712             please_cite(fplog, "Engin2010");
1713         }
1714     }
1715
1716     pull->bRefAt   = FALSE;
1717     pull->cosdim   = -1;
1718     for (g = 0; g < pull->ngroup; g++)
1719     {
1720         pull_group_work_t *pgrp;
1721
1722         pgrp           = &pull->group[g];
1723         pgrp->epgrppbc = epgrppbcNONE;
1724         if (pgrp->params.nat > 0)
1725         {
1726             /* There is an issue when a group is used in multiple coordinates
1727              * and constraints are applied in different dimensions with atoms
1728              * frozen in some, but not all dimensions.
1729              * Since there is only one mass array per group, we can't have
1730              * frozen/non-frozen atoms for different coords at the same time.
1731              * But since this is a very exotic case, we don't check for this.
1732              * A warning is printed in init_pull_group_index.
1733              */
1734
1735             gmx_bool bConstraint;
1736             ivec     pulldim, pulldim_con;
1737
1738             /* Loop over all pull coordinates to see along which dimensions
1739              * this group is pulled and if it is involved in constraints.
1740              */
1741             bConstraint = FALSE;
1742             clear_ivec(pulldim);
1743             clear_ivec(pulldim_con);
1744             for (c = 0; c < pull->ncoord; c++)
1745             {
1746                 if (pull->coord[c].params.group[0] == g ||
1747                     pull->coord[c].params.group[1] == g ||
1748                     (pull->coord[c].params.eGeom == epullgDIRRELATIVE &&
1749                      (pull->coord[c].params.group[2] == g ||
1750                       pull->coord[c].params.group[3] == g)))
1751                 {
1752                     for (m = 0; m < DIM; m++)
1753                     {
1754                         if (pull->coord[c].params.dim[m] == 1)
1755                         {
1756                             pulldim[m] = 1;
1757
1758                             if (pull->coord[c].params.eType == epullCONSTRAINT)
1759                             {
1760                                 bConstraint    = TRUE;
1761                                 pulldim_con[m] = 1;
1762                             }
1763                         }
1764                     }
1765                 }
1766             }
1767
1768             /* Determine if we need to take PBC into account for calculating
1769              * the COM's of the pull groups.
1770              */
1771             GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
1772             for (m = 0; m < pull->npbcdim; m++)
1773             {
1774                 GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
1775                 if (pulldim[m] == 1 && pgrp->params.nat > 1)
1776                 {
1777                     if (pgrp->params.pbcatom >= 0)
1778                     {
1779                         pgrp->epgrppbc = epgrppbcREFAT;
1780                         pull->bRefAt   = TRUE;
1781                     }
1782                     else
1783                     {
1784                         if (pgrp->params.weight != NULL)
1785                         {
1786                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1787                         }
1788                         pgrp->epgrppbc = epgrppbcCOS;
1789                         if (pull->cosdim >= 0 && pull->cosdim != m)
1790                         {
1791                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
1792                         }
1793                         pull->cosdim = m;
1794                     }
1795                 }
1796             }
1797             /* Set the indices */
1798             init_pull_group_index(fplog, cr, g, pgrp,
1799                                   bConstraint, pulldim_con,
1800                                   mtop, ir, lambda);
1801         }
1802         else
1803         {
1804             /* Absolute reference, set the inverse mass to zero.
1805              * This is only relevant (and used) with constraint pulling.
1806              */
1807             pgrp->invtm  = 0;
1808             pgrp->wscale = 1;
1809         }
1810     }
1811
1812     /* If we use cylinder coordinates, do some initialising for them */
1813     if (pull->bCylinder)
1814     {
1815         snew(pull->dyna, pull->ncoord);
1816
1817         for (c = 0; c < pull->ncoord; c++)
1818         {
1819             const pull_coord_work_t *pcrd;
1820
1821             pcrd = &pull->coord[c];
1822
1823             if (pcrd->params.eGeom == epullgCYL)
1824             {
1825                 if (pull->group[pcrd->params.group[0]].params.nat == 0)
1826                 {
1827                     gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
1828                 }
1829             }
1830         }
1831     }
1832
1833     comm = &pull->comm;
1834
1835 #ifdef GMX_MPI
1836     /* Use a sub-communicator when we have more than 32 ranks */
1837     comm->bParticipateAll = (cr == NULL || !DOMAINDECOMP(cr) ||
1838                              cr->dd->nnodes <= 32 ||
1839                              getenv("GMX_PULL_PARTICIPATE_ALL") != NULL);
1840     /* This sub-commicator is not used with comm->bParticipateAll,
1841      * so we can always initialize it to NULL.
1842      */
1843     comm->mpi_comm_com    = MPI_COMM_NULL;
1844     comm->nparticipate    = 0;
1845 #else
1846     /* No MPI: 1 rank: all ranks pull */
1847     comm->bParticipateAll = TRUE;
1848 #endif
1849     comm->bParticipate    = comm->bParticipateAll;
1850     comm->setup_count     = 0;
1851     comm->must_count      = 0;
1852
1853     if (!comm->bParticipateAll && fplog != NULL)
1854     {
1855         fprintf(fplog, "Will use a sub-communicator for pull communication\n");
1856     }
1857
1858     comm->rbuf     = NULL;
1859     comm->dbuf     = NULL;
1860     comm->dbuf_cyl = NULL;
1861
1862     /* We still need to initialize the PBC reference coordinates */
1863     pull->bSetPBCatoms = TRUE;
1864
1865     /* Only do I/O when we are doing dynamics and if we are the MASTER */
1866     pull->out_x = NULL;
1867     pull->out_f = NULL;
1868     if (bOutFile)
1869     {
1870         if (pull->params.nstxout > 0)
1871         {
1872             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
1873                                         TRUE, Flags);
1874         }
1875         if (pull->params.nstfout > 0)
1876         {
1877             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1878                                         FALSE, Flags);
1879         }
1880     }
1881
1882     return pull;
1883 }
1884
1885 static void destroy_pull_group(pull_group_work_t *pgrp)
1886 {
1887     if (pgrp->ind_loc != pgrp->params.ind)
1888     {
1889         sfree(pgrp->ind_loc);
1890     }
1891     if (pgrp->weight_loc != pgrp->params.weight)
1892     {
1893         sfree(pgrp->weight_loc);
1894     }
1895     sfree(pgrp->mdw);
1896     sfree(pgrp->dv);
1897
1898     sfree(pgrp->params.ind);
1899     sfree(pgrp->params.weight);
1900 }
1901
1902 static void destroy_pull(struct pull_t *pull)
1903 {
1904     int i;
1905
1906     for (i = 0; i < pull->ngroup; i++)
1907     {
1908         destroy_pull_group(&pull->group[i]);
1909     }
1910     for (i = 0; i < pull->ncoord; i++)
1911     {
1912         if (pull->coord[i].params.eGeom == epullgCYL)
1913         {
1914             destroy_pull_group(&(pull->dyna[i]));
1915         }
1916     }
1917     sfree(pull->group);
1918     sfree(pull->dyna);
1919     sfree(pull->coord);
1920
1921 #ifdef GMX_MPI
1922     if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
1923     {
1924         MPI_Comm_free(&pull->comm.mpi_comm_com);
1925     }
1926 #endif
1927     sfree(pull->comm.rbuf);
1928     sfree(pull->comm.dbuf);
1929     sfree(pull->comm.dbuf_cyl);
1930
1931     sfree(pull);
1932 }
1933
1934 void finish_pull(struct pull_t *pull)
1935 {
1936     if (pull->out_x)
1937     {
1938         gmx_fio_fclose(pull->out_x);
1939     }
1940     if (pull->out_f)
1941     {
1942         gmx_fio_fclose(pull->out_f);
1943     }
1944
1945     destroy_pull(pull);
1946 }
1947
1948 gmx_bool pull_have_potential(const struct pull_t *pull)
1949 {
1950     return pull->bPotential;
1951 }
1952
1953 gmx_bool pull_have_constraint(const struct pull_t *pull)
1954 {
1955     return pull->bConstraint;
1956 }