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     if (pgrp->nat == 1 && pgrp->nat_loc == 1)
 306     {
 307         /* Only one atom and our rank has this atom: we can skip
 308          * the mass weighting, which means that this code also works
 309          * for mass=0, e.g. with a virtual site.
 310          */
 311         for (m = 0; m < DIM; m++)
 312         {
 313             f[pgrp->ind_loc[0]][m] += sign*f_pull[m];
 314         }
 315     }
 316     else
 317     {
 318         for (i = 0; i < pgrp->nat_loc; i++)
 319         {
 320             ii    = pgrp->ind_loc[i];
 321             wmass = md->massT[ii];
 322             if (pgrp->weight_loc)
 323             {
 324                 wmass *= pgrp->weight_loc[i];
 325             }
 326
 327             for (m = 0; m < DIM; m++)
 328             {
 329                 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
 330             }
 331         }
 332     }
 333 }
 334
 335 /* Apply forces in a mass weighted fashion to a cylinder group */
 336 static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
 337                                  const t_mdatoms *md,
 338                                  const dvec f_pull, double f_scal,
 339                                  int sign, rvec *f)
 340 {
 341     int    i, ii, m;
 342     double mass, weight, inv_wm, dv_com;
 343
 344     inv_wm = pgrp->mwscale;
 345
 346     for (i = 0; i < pgrp->nat_loc; i++)
 347     {
 348         ii     = pgrp->ind_loc[i];
 349         mass   = md->massT[ii];
 350         weight = pgrp->weight_loc[i];
 351         /* The stored axial distance from the cylinder center (dv) needs
 352          * to be corrected for an offset (dv_corr), which was unknown when
 353          * we calculated dv.
 354          */
 355         dv_com = pgrp->dv[i] + dv_corr;
 356
 357         /* Here we not only add the pull force working along vec (f_pull),
 358          * but also a radial component, due to the dependence of the weights
 359          * on the radial distance.
 360          */
 361         for (m = 0; m < DIM; m++)
 362         {
 363             f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
 364                                      pgrp->mdw[i][m]*dv_com*f_scal);
 365         }
 366     }
 367 }
 368
 369 /* Apply torque forces in a mass weighted fashion to the groups that define
 370  * the pull vector direction for pull coordinate pcrd.
 371  */
 372 static void apply_forces_vec_torque(const t_pull       *pull,
 373                                     const t_pull_coord *pcrd,
 374                                     const t_mdatoms    *md,
 375                                     rvec               *f)
 376 {
 377     double inpr;
 378     int    m;
 379     dvec   f_perp;
 380
 381     /* The component inpr along the pull vector is accounted for in the usual
 382      * way. Here we account for the component perpendicular to vec.
 383      */
 384     inpr = 0;
 385     for (m = 0; m < DIM; m++)
 386     {
 387         inpr += pcrd->dr[m]*pcrd->vec[m];
 388     }
 389     /* The torque force works along the component of the distance vector
 390      * of between the two "usual" pull groups that is perpendicular to
 391      * the pull vector. The magnitude of this force is the "usual" scale force
 392      * multiplied with the ratio of the distance between the two "usual" pull
 393      * groups and the distance between the two groups that define the vector.
 394      */
 395     for (m = 0; m < DIM; m++)
 396     {
 397         f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal;
 398     }
 399
 400     /* Apply the force to the groups defining the vector using opposite signs */
 401     apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f);
 402     apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp,  1, f);
 403 }
 404
 405 /* Apply forces in a mass weighted fashion */
 406 static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
 407 {
 408     int                 c;
 409     const t_pull_coord *pcrd;
 410
 411     for (c = 0; c < pull->ncoord; c++)
 412     {
 413         pcrd = &pull->coord[c];
 414
 415         if (pcrd->eGeom == epullgCYL)
 416         {
 417             dvec f_tot;
 418             int  m;
 419
 420             apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
 421                                  pcrd->f, pcrd->f_scal, -1, f);
 422
 423             /* Sum the force along the vector and the radial force */
 424             for (m = 0; m < DIM; m++)
 425             {
 426                 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
 427             }
 428             apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f);
 429         }
 430         else
 431         {
 432             if (pcrd->eGeom == epullgDIRRELATIVE)
 433             {
 434                 /* We need to apply the torque forces to the pull groups
 435                  * that define the pull vector.
 436                  */
 437                 apply_forces_vec_torque(pull, pcrd, md, f);
 438             }
 439
 440             if (pull->group[pcrd->group[0]].nat > 0)
 441             {
 442                 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
 443             }
 444             apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
 445         }
 446     }
 447 }
 448
 449 static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
 450 {
 451     double max_d2;
 452     int    m;
 453
 454     max_d2 = GMX_DOUBLE_MAX;
 455
 456     for (m = 0; m < pbc->ndim_ePBC; m++)
 457     {
 458         if (pcrd->dim[m] != 0)
 459         {
 460             max_d2 = min(max_d2, norm2(pbc->box[m]));
 461         }
 462     }
 463
 464     return 0.25*max_d2;
 465 }
 466
 467 /* This function returns the distance based on coordinates xg and xref.
 468  * Note that the pull coordinate struct pcrd is not modified.
 469  */
 470 static void low_get_pull_coord_dr(const t_pull *pull,
 471                                   const t_pull_coord *pcrd,
 472                                   const t_pbc *pbc, double t,
 473                                   dvec xg, dvec xref, double max_dist2,
 474                                   dvec dr)
 475 {
 476     const t_pull_group *pgrp0, *pgrp1;
 477     int                 m;
 478     dvec                xrefr, dref = {0, 0, 0};
 479     double              dr2;
 480
 481     pgrp0 = &pull->group[pcrd->group[0]];
 482     pgrp1 = &pull->group[pcrd->group[1]];
 483
 484     /* Only the first group can be an absolute reference, in that case nat=0 */
 485     if (pgrp0->nat == 0)
 486     {
 487         for (m = 0; m < DIM; m++)
 488         {
 489             xref[m] = pcrd->origin[m];
 490         }
 491     }
 492
 493     copy_dvec(xref, xrefr);
 494
 495     if (pcrd->eGeom == epullgDIRPBC)
 496     {
 497         for (m = 0; m < DIM; m++)
 498         {
 499             dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
 500         }
 501         /* Add the reference position, so we use the correct periodic image */
 502         dvec_inc(xrefr, dref);
 503     }
 504
 505     pbc_dx_d(pbc, xg, xrefr, dr);
 506     dr2 = 0;
 507     for (m = 0; m < DIM; m++)
 508     {
 509         dr[m] *= pcrd->dim[m];
 510         dr2   += dr[m]*dr[m];
 511     }
 512     if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
 513     {
 514         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",
 515                   pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
 516     }
 517
 518     if (pcrd->eGeom == epullgDIRPBC)
 519     {
 520         dvec_inc(dr, dref);
 521     }
 522 }
 523
 524 /* This function returns the distance based on the contents of the pull struct.
 525  * pull->coord[coord_ind].dr, and potentially vector, are updated.
 526  */
 527 static void get_pull_coord_dr(const t_pull *pull,
 528                               int coord_ind,
 529                               const t_pbc *pbc, double t)
 530 {
 531     double        md2;
 532     t_pull_coord *pcrd;
 533
 534     pcrd = &pull->coord[coord_ind];
 535
 536     if (pcrd->eGeom == epullgDIRPBC)
 537     {
 538         md2 = -1;
 539     }
 540     else
 541     {
 542         md2 = max_pull_distance2(pcrd, pbc);
 543     }
 544
 545     if (pcrd->eGeom == epullgDIRRELATIVE)
 546     {
 547         /* We need to determine the pull vector */
 548         const t_pull_group *pgrp2, *pgrp3;
 549         dvec                vec;
 550         int                 m;
 551         double              invlen;
 552
 553         pgrp2 = &pull->group[pcrd->group[2]];
 554         pgrp3 = &pull->group[pcrd->group[3]];
 555
 556         pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
 557
 558         for (m = 0; m < DIM; m++)
 559         {
 560             vec[m] *= pcrd->dim[m];
 561         }
 562         pcrd->vec_len = dnorm(vec);
 563         for (m = 0; m < DIM; m++)
 564         {
 565             pcrd->vec[m] = vec[m]/pcrd->vec_len;
 566         }
 567         if (debug)
 568         {
 569             fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
 570                     coord_ind,
 571                     vec[XX], vec[YY], vec[ZZ],
 572                     pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]);
 573         }
 574     }
 575
 576     low_get_pull_coord_dr(pull, pcrd, pbc, t,
 577                           pull->group[pcrd->group[1]].x,
 578                           pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
 579                           md2,
 580                           pcrd->dr);
 581 }
 582
 583 void get_pull_coord_distance(const t_pull *pull,
 584                              int coord_ind,
 585                              const t_pbc *pbc, double t,
 586                              double *dev)
 587 {
 588     static gmx_bool     bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
 589                                             but is fairly benign */
 590     const t_pull_coord *pcrd;
 591     int                 m;
 592     double              ref, drs, inpr;
 593
 594     pcrd = &pull->coord[coord_ind];
 595
 596     get_pull_coord_dr(pull, coord_ind, pbc, t);
 597
 598     ref = pcrd->init + pcrd->rate*t;
 599
 600     switch (pcrd->eGeom)
 601     {
 602         case epullgDIST:
 603             /* Pull along the vector between the com's */
 604             if (ref < 0 && !bWarned)
 605             {
 606                 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
 607                 bWarned = TRUE;
 608             }
 609             drs = dnorm(pcrd->dr);
 610             if (drs == 0)
 611             {
 612                 /* With no vector we can not determine the direction for the force,
 613                  * so we set the force to zero.
 614                  */
 615                 *dev = 0;
 616             }
 617             else
 618             {
 619                 /* Determine the deviation */
 620                 *dev = drs - ref;
 621             }
 622             break;
 623         case epullgDIR:
 624         case epullgDIRPBC:
 625         case epullgDIRRELATIVE:
 626         case epullgCYL:
 627             /* Pull along vec */
 628             inpr = 0;
 629             for (m = 0; m < DIM; m++)
 630             {
 631                 inpr += pcrd->vec[m]*pcrd->dr[m];
 632             }
 633             *dev = inpr - ref;
 634             break;
 635     }
 636 }
 637
 638 void clear_pull_forces(t_pull *pull)
 639 {
 640     int i;
 641
 642     /* Zeroing the forces is only required for constraint pulling.
 643      * It can happen that multiple constraint steps need to be applied
 644      * and therefore the constraint forces need to be accumulated.
 645      */
 646     for (i = 0; i < pull->ncoord; i++)
 647     {
 648         clear_dvec(pull->coord[i].f);
 649         pull->coord[i].f_scal = 0;
 650     }
 651 }
 652
 653 /* Apply constraint using SHAKE */
 654 static void do_constraint(t_pull *pull, t_pbc *pbc,
 655                           rvec *x, rvec *v,
 656                           gmx_bool bMaster, tensor vir,
 657                           double dt, double t)
 658 {
 659
 660     dvec         *r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
 661     dvec          unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
 662     dvec         *rnew;   /* current 'new' positions of the groups */
 663     double       *dr_tot; /* the total update of the coords */
 664     double        ref;
 665     dvec          vec;
 666     double        inpr;
 667     double        lambda, rm, invdt = 0;
 668     gmx_bool      bConverged_all, bConverged = FALSE;
 669     int           niter = 0, g, c, ii, j, m, max_iter = 100;
 670     double        a;
 671     dvec          f;       /* the pull force */
 672     dvec          tmp, tmp3;
 673     t_pull_group *pgrp0, *pgrp1;
 674     t_pull_coord *pcrd;
 675
 676     snew(r_ij,   pull->ncoord);
 677     snew(dr_tot, pull->ncoord);
 678
 679     snew(rnew, pull->ngroup);
 680
 681     /* copy the current unconstrained positions for use in iterations. We
 682        iterate until rinew[i] and rjnew[j] obey the constraints. Then
 683        rinew - pull.x_unc[i] is the correction dr to group i */
 684     for (g = 0; g < pull->ngroup; g++)
 685     {
 686         copy_dvec(pull->group[g].xp, rnew[g]);
 687     }
 688
 689     /* Determine the constraint directions from the old positions */
 690     for (c = 0; c < pull->ncoord; c++)
 691     {
 692         pcrd = &pull->coord[c];
 693
 694         if (pcrd->eType != epullCONSTRAINT)
 695         {
 696             continue;
 697         }
 698
 699         get_pull_coord_dr(pull, c, pbc, t);
 700         /* Store the difference vector at time t for printing */
 701         if (debug)
 702         {
 703             fprintf(debug, "Pull coord %d dr %f %f %f\n",
 704                     c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
 705         }
 706
 707         if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC)
 708         {
 709             /* Select the component along vec */
 710             a = 0;
 711             for (m = 0; m < DIM; m++)
 712             {
 713                 a += pcrd->vec[m]*pcrd->dr[m];
 714             }
 715             for (m = 0; m < DIM; m++)
 716             {
 717                 r_ij[c][m] = a*pcrd->vec[m];
 718             }
 719         }
 720         else
 721         {
 722             copy_dvec(pcrd->dr, r_ij[c]);
 723         }
 724
 725         if (dnorm2(r_ij[c]) == 0)
 726         {
 727             gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
 728         }
 729     }
 730
 731     bConverged_all = FALSE;
 732     while (!bConverged_all && niter < max_iter)
 733     {
 734         bConverged_all = TRUE;
 735
 736         /* loop over all constraints */
 737         for (c = 0; c < pull->ncoord; c++)
 738         {
 739             dvec dr0, dr1;
 740
 741             pcrd  = &pull->coord[c];
 742
 743             if (pcrd->eType != epullCONSTRAINT)
 744             {
 745                 continue;
 746             }
 747
 748             pgrp0 = &pull->group[pcrd->group[0]];
 749             pgrp1 = &pull->group[pcrd->group[1]];
 750
 751             /* Get the current difference vector */
 752             low_get_pull_coord_dr(pull, pcrd, pbc, t,
 753                                   rnew[pcrd->group[1]],
 754                                   rnew[pcrd->group[0]],
 755                                   -1, unc_ij);
 756
 757             ref = pcrd->init + pcrd->rate*t;
 758
 759             if (debug)
 760             {
 761                 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
 762             }
 763
 764             rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
 765
 766             switch (pcrd->eGeom)
 767             {
 768                 case epullgDIST:
 769                     if (ref <= 0)
 770                     {
 771                         gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
 772                     }
 773
 774                     {
 775                         double q, c_a, c_b, c_c;
 776
 777                         c_a = diprod(r_ij[c], r_ij[c]);
 778                         c_b = diprod(unc_ij, r_ij[c])*2;
 779                         c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
 780
 781                         if (c_b < 0)
 782                         {
 783                             q      = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
 784                             lambda = -q/c_a;
 785                         }
 786                         else
 787                         {
 788                             q      = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
 789                             lambda = -c_c/q;
 790                         }
 791
 792                         if (debug)
 793                         {
 794                             fprintf(debug,
 795                                     "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
 796                                     c_a, c_b, c_c, lambda);
 797                         }
 798                     }
 799
 800                     /* The position corrections dr due to the constraints */
 801                     dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
 802                     dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
 803                     dr_tot[c] += -lambda*dnorm(r_ij[c]);
 804                     break;
 805                 case epullgDIR:
 806                 case epullgDIRPBC:
 807                 case epullgCYL:
 808                     /* A 1-dimensional constraint along a vector */
 809                     a = 0;
 810                     for (m = 0; m < DIM; m++)
 811                     {
 812                         vec[m] = pcrd->vec[m];
 813                         a     += unc_ij[m]*vec[m];
 814                     }
 815                     /* Select only the component along the vector */
 816                     dsvmul(a, vec, unc_ij);
 817                     lambda = a - ref;
 818                     if (debug)
 819                     {
 820                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
 821                     }
 822
 823                     /* The position corrections dr due to the constraints */
 824                     dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
 825                     dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
 826                     dr_tot[c] += -lambda;
 827                     break;
 828             }
 829
 830             /* DEBUG */
 831             if (debug)
 832             {
 833                 int g0, g1;
 834
 835                 g0 = pcrd->group[0];
 836                 g1 = pcrd->group[1];
 837                 low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
 838                 low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
 839                 fprintf(debug,
 840                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 841                         rnew[g0][0], rnew[g0][1], rnew[g0][2],
 842                         rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
 843                 fprintf(debug,
 844                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
 845                         "", "", "", "", "", "", ref);
 846                 fprintf(debug,
 847                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 848                         dr0[0], dr0[1], dr0[2],
 849                         dr1[0], dr1[1], dr1[2],
 850                         dnorm(tmp3));
 851             } /* END DEBUG */
 852
 853             /* Update the COMs with dr */
 854             dvec_inc(rnew[pcrd->group[1]], dr1);
 855             dvec_inc(rnew[pcrd->group[0]], dr0);
 856         }
 857
 858         /* Check if all constraints are fullfilled now */
 859         for (c = 0; c < pull->ncoord; c++)
 860         {
 861             pcrd = &pull->coord[c];
 862
 863             if (pcrd->eType != epullCONSTRAINT)
 864             {
 865                 continue;
 866             }
 867
 868             ref  = pcrd->init + pcrd->rate*t;
 869
 870             low_get_pull_coord_dr(pull, pcrd, pbc, t,
 871                                   rnew[pcrd->group[1]],
 872                                   rnew[pcrd->group[0]],
 873                                   -1, unc_ij);
 874
 875             switch (pcrd->eGeom)
 876             {
 877                 case epullgDIST:
 878                     bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
 879                     break;
 880                 case epullgDIR:
 881                 case epullgDIRPBC:
 882                 case epullgCYL:
 883                     for (m = 0; m < DIM; m++)
 884                     {
 885                         vec[m] = pcrd->vec[m];
 886                     }
 887                     inpr = diprod(unc_ij, vec);
 888                     dsvmul(inpr, vec, unc_ij);
 889                     bConverged =
 890                         fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
 891                     break;
 892             }
 893
 894             if (!bConverged)
 895             {
 896                 if (debug)
 897                 {
 898                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
 899                             "d_ref = %f, current d = %f\n",
 900                             g, ref, dnorm(unc_ij));
 901                 }
 902
 903                 bConverged_all = FALSE;
 904             }
 905         }
 906
 907         niter++;
 908         /* if after all constraints are dealt with and bConverged is still TRUE
 909            we're finished, if not we do another iteration */
 910     }
 911     if (niter > max_iter)
 912     {
 913         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
 914     }
 915
 916     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 917
 918     if (v)
 919     {
 920         invdt = 1/dt;
 921     }
 922
 923     /* update atoms in the groups */
 924     for (g = 0; g < pull->ngroup; g++)
 925     {
 926         const t_pull_group *pgrp;
 927         dvec                dr;
 928
 929         pgrp = &pull->group[g];
 930
 931         /* get the final constraint displacement dr for group g */
 932         dvec_sub(rnew[g], pgrp->xp, dr);
 933
 934         if (dnorm2(dr) == 0)
 935         {
 936             /* No displacement, no update necessary */
 937             continue;
 938         }
 939
 940         /* update the atom positions */
 941         copy_dvec(dr, tmp);
 942         for (j = 0; j < pgrp->nat_loc; j++)
 943         {
 944             ii = pgrp->ind_loc[j];
 945             if (pgrp->weight_loc)
 946             {
 947                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
 948             }
 949             for (m = 0; m < DIM; m++)
 950             {
 951                 x[ii][m] += tmp[m];
 952             }
 953             if (v)
 954             {
 955                 for (m = 0; m < DIM; m++)
 956                 {
 957                     v[ii][m] += invdt*tmp[m];
 958                 }
 959             }
 960         }
 961     }
 962
 963     /* calculate the constraint forces, used for output and virial only */
 964     for (c = 0; c < pull->ncoord; c++)
 965     {
 966         pcrd = &pull->coord[c];
 967
 968         if (pcrd->eType != epullCONSTRAINT)
 969         {
 970             continue;
 971         }
 972
 973         pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
 974
 975         if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster)
 976         {
 977             double f_invr;
 978
 979             /* Add the pull contribution to the virial */
 980             /* We have already checked above that r_ij[c] != 0 */
 981             f_invr = pcrd->f_scal/dnorm(r_ij[c]);
 982
 983             for (j = 0; j < DIM; j++)
 984             {
 985                 for (m = 0; m < DIM; m++)
 986                 {
 987                     vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
 988                 }
 989             }
 990         }
 991     }
 992
 993     /* finished! I hope. Give back some memory */
 994     sfree(r_ij);
 995     sfree(dr_tot);
 996     sfree(rnew);
 997 }
 998
 999 /* Pulling with a harmonic umbrella potential or constant force */
1000 static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
1001                         real *V, tensor vir, real *dVdl)
1002 {
1003     int           c, j, m;
1004     double        dev, ndr, invdr = 0;
1005     real          k, dkdl;
1006     t_pull_coord *pcrd;
1007
1008     /* loop over the pull coordinates */
1009     *V    = 0;
1010     *dVdl = 0;
1011     for (c = 0; c < pull->ncoord; c++)
1012     {
1013         pcrd = &pull->coord[c];
1014
1015         if (pcrd->eType == epullCONSTRAINT)
1016         {
1017             continue;
1018         }
1019
1020         get_pull_coord_distance(pull, c, pbc, t, &dev);
1021
1022         k    = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
1023         dkdl = pcrd->kB - pcrd->k;
1024
1025         if (pcrd->eGeom == epullgDIST)
1026         {
1027             ndr   = dnorm(pcrd->dr);
1028             if (ndr > 0)
1029             {
1030                 invdr = 1/ndr;
1031             }
1032             else
1033             {
1034                 /* With an harmonic umbrella, the force is 0 at r=0,
1035                  * so we can set invdr to any value.
1036                  * With a constant force, the force at r=0 is not defined,
1037                  * so we zero it (this is anyhow a very rare event).
1038                  */
1039                 invdr = 0;
1040             }
1041         }
1042         else
1043         {
1044             ndr = 0;
1045             for (m = 0; m < DIM; m++)
1046             {
1047                 ndr += pcrd->vec[m]*pcrd->dr[m];
1048             }
1049         }
1050
1051         switch (pcrd->eType)
1052         {
1053             case epullUMBRELLA:
1054             case epullFLATBOTTOM:
1055                 /* The only difference between an umbrella and a flat-bottom
1056                  * potential is that a flat-bottom is zero below zero.
1057                  */
1058                 if (pcrd->eType == epullFLATBOTTOM && dev < 0)
1059                 {
1060                     dev = 0;
1061                 }
1062
1063                 pcrd->f_scal  =       -k*dev;
1064                 *V           += 0.5*   k*dsqr(dev);
1065                 *dVdl        += 0.5*dkdl*dsqr(dev);
1066                 break;
1067             case epullCONST_F:
1068                 pcrd->f_scal  =   -k;
1069                 *V           +=    k*ndr;
1070                 *dVdl        += dkdl*ndr;
1071                 break;
1072             default:
1073                 gmx_incons("Unsupported pull type in do_pull_pot");
1074         }
1075
1076         if (pcrd->eGeom == epullgDIST)
1077         {
1078             for (m = 0; m < DIM; m++)
1079             {
1080                 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
1081             }
1082         }
1083         else
1084         {
1085             for (m = 0; m < DIM; m++)
1086             {
1087                 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
1088             }
1089         }
1090
1091         if (vir != NULL && pcrd->eGeom != epullgDIRPBC)
1092         {
1093             /* Add the pull contribution to the virial */
1094             for (j = 0; j < DIM; j++)
1095             {
1096                 for (m = 0; m < DIM; m++)
1097                 {
1098                     vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
1099                 }
1100             }
1101         }
1102     }
1103 }
1104
1105 real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1106                     t_commrec *cr, double t, real lambda,
1107                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
1108 {
1109     real V, dVdl;
1110
1111     pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1112
1113     do_pull_pot(pull, pbc, t, lambda,
1114                 &V, MASTER(cr) ? vir : NULL, &dVdl);
1115
1116     /* Distribute forces over pulled groups */
1117     apply_forces(pull, md, f);
1118
1119     if (MASTER(cr))
1120     {
1121         *dvdlambda += dVdl;
1122     }
1123
1124     return (MASTER(cr) ? V : 0.0);
1125 }
1126
1127 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1128                      t_commrec *cr, double dt, double t,
1129                      rvec *x, rvec *xp, rvec *v, tensor vir)
1130 {
1131     pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1132
1133     do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1134 }
1135
1136 static void make_local_pull_group(gmx_ga2la_t ga2la,
1137                                   t_pull_group *pg, int start, int end)
1138 {
1139     int i, ii;
1140
1141     pg->nat_loc = 0;
1142     for (i = 0; i < pg->nat; i++)
1143     {
1144         ii = pg->ind[i];
1145         if (ga2la)
1146         {
1147             if (!ga2la_get_home(ga2la, ii, &ii))
1148             {
1149                 ii = -1;
1150             }
1151         }
1152         if (ii >= start && ii < end)
1153         {
1154             /* This is a home atom, add it to the local pull group */
1155             if (pg->nat_loc >= pg->nalloc_loc)
1156             {
1157                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1158                 srenew(pg->ind_loc, pg->nalloc_loc);
1159                 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
1160                 {
1161                     srenew(pg->weight_loc, pg->nalloc_loc);
1162                 }
1163             }
1164             pg->ind_loc[pg->nat_loc] = ii;
1165             if (pg->weight)
1166             {
1167                 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1168             }
1169             pg->nat_loc++;
1170         }
1171     }
1172 }
1173
1174 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
1175 {
1176     gmx_ga2la_t ga2la;
1177     int         g;
1178
1179     if (dd)
1180     {
1181         ga2la = dd->ga2la;
1182     }
1183     else
1184     {
1185         ga2la = NULL;
1186     }
1187
1188     for (g = 0; g < pull->ngroup; g++)
1189     {
1190         make_local_pull_group(ga2la, &pull->group[g],
1191                               0, md->homenr);
1192     }
1193
1194     /* Since the PBC of atoms might have changed, we need to update the PBC */
1195     pull->bSetPBCatoms = TRUE;
1196 }
1197
1198 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1199                                   int g, t_pull_group *pg,
1200                                   gmx_bool bConstraint, ivec pulldim_con,
1201                                   gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
1202 {
1203     int                   i, ii, d, nfrozen, ndim;
1204     real                  m, w, mbd;
1205     double                tmass, wmass, wwmass;
1206     gmx_groups_t         *groups;
1207     gmx_mtop_atomlookup_t alook;
1208     t_atom               *atom;
1209
1210     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1211     {
1212         /* There are no masses in the integrator.
1213          * But we still want to have the correct mass-weighted COMs.
1214          * So we store the real masses in the weights.
1215          * We do not set nweight, so these weights do not end up in the tpx file.
1216          */
1217         if (pg->nweight == 0)
1218         {
1219             snew(pg->weight, pg->nat);
1220         }
1221     }
1222
1223     if (cr && PAR(cr))
1224     {
1225         pg->nat_loc    = 0;
1226         pg->nalloc_loc = 0;
1227         pg->ind_loc    = NULL;
1228         pg->weight_loc = NULL;
1229     }
1230     else
1231     {
1232         pg->nat_loc = pg->nat;
1233         pg->ind_loc = pg->ind;
1234         if (pg->epgrppbc == epgrppbcCOS)
1235         {
1236             snew(pg->weight_loc, pg->nat);
1237         }
1238         else
1239         {
1240             pg->weight_loc = pg->weight;
1241         }
1242     }
1243
1244     groups = &mtop->groups;
1245
1246     alook = gmx_mtop_atomlookup_init(mtop);
1247
1248     nfrozen = 0;
1249     tmass   = 0;
1250     wmass   = 0;
1251     wwmass  = 0;
1252     for (i = 0; i < pg->nat; i++)
1253     {
1254         ii = pg->ind[i];
1255         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1256         if (bConstraint && ir->opts.nFreeze)
1257         {
1258             for (d = 0; d < DIM; d++)
1259             {
1260                 if (pulldim_con[d] == 1 &&
1261                     ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1262                 {
1263                     nfrozen++;
1264                 }
1265             }
1266         }
1267         if (ir->efep == efepNO)
1268         {
1269             m = atom->m;
1270         }
1271         else
1272         {
1273             m = (1 - lambda)*atom->m + lambda*atom->mB;
1274         }
1275         if (pg->nweight > 0)
1276         {
1277             w = pg->weight[i];
1278         }
1279         else
1280         {
1281             w = 1;
1282         }
1283         if (EI_ENERGY_MINIMIZATION(ir->eI))
1284         {
1285             /* Move the mass to the weight */
1286             w            *= m;
1287             m             = 1;
1288             pg->weight[i] = w;
1289         }
1290         else if (ir->eI == eiBD)
1291         {
1292             if (ir->bd_fric)
1293             {
1294                 mbd = ir->bd_fric*ir->delta_t;
1295             }
1296             else
1297             {
1298                 if (groups->grpnr[egcTC] == NULL)
1299                 {
1300                     mbd = ir->delta_t/ir->opts.tau_t[0];
1301                 }
1302                 else
1303                 {
1304                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1305                 }
1306             }
1307             w            *= m/mbd;
1308             m             = mbd;
1309             pg->weight[i] = w;
1310         }
1311         tmass  += m;
1312         wmass  += m*w;
1313         wwmass += m*w*w;
1314     }
1315
1316     gmx_mtop_atomlookup_destroy(alook);
1317
1318     if (wmass == 0)
1319     {
1320         /* We can have single atom groups with zero mass with potential pulling
1321          * without cosine weighting.
1322          */
1323         if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1324         {
1325             /* With one atom the mass doesn't matter */
1326             wwmass = 1;
1327         }
1328         else
1329         {
1330             gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1331                       pg->weight ? " weighted" : "", g);
1332         }
1333     }
1334     if (fplog)
1335     {
1336         fprintf(fplog,
1337                 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1338         if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1339         {
1340             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1341         }
1342         if (pg->epgrppbc == epgrppbcCOS)
1343         {
1344             fprintf(fplog, ", cosine weighting will be used");
1345         }
1346         fprintf(fplog, "\n");
1347     }
1348
1349     if (nfrozen == 0)
1350     {
1351         /* A value != 0 signals not frozen, it is updated later */
1352         pg->invtm  = -1.0;
1353     }
1354     else
1355     {
1356         ndim = 0;
1357         for (d = 0; d < DIM; d++)
1358         {
1359             ndim += pulldim_con[d]*pg->nat;
1360         }
1361         if (fplog && nfrozen > 0 && nfrozen < ndim)
1362         {
1363             fprintf(fplog,
1364                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1365                     "         that are subject to pulling are frozen.\n"
1366                     "         For constraint pulling the whole group will be frozen.\n\n",
1367                     g);
1368         }
1369         pg->invtm  = 0.0;
1370         pg->wscale = 1.0;
1371     }
1372 }
1373
1374 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1375                gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1376                gmx_bool bOutFile, unsigned long Flags)
1377 {
1378     t_pull       *pull;
1379     t_pull_group *pgrp;
1380     int           c, g, m;
1381
1382     pull = ir->pull;
1383
1384     pull->bPotential  = FALSE;
1385     pull->bConstraint = FALSE;
1386     pull->bCylinder   = FALSE;
1387     for (c = 0; c < pull->ncoord; c++)
1388     {
1389         t_pull_coord *pcrd;
1390         int           calc_com_start, calc_com_end, g;
1391
1392         pcrd = &pull->coord[c];
1393
1394         if (pcrd->eType == epullCONSTRAINT)
1395         {
1396             /* Check restrictions of the constraint pull code */
1397             if (pcrd->eGeom == epullgCYL ||
1398                 pcrd->eGeom == epullgDIRRELATIVE)
1399             {
1400                 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
1401                           epull_names[pcrd->eType],
1402                           epullg_names[pcrd->eGeom],
1403                           epull_names[epullUMBRELLA]);
1404             }
1405
1406             pull->bConstraint = TRUE;
1407         }
1408         else
1409         {
1410             pull->bPotential = TRUE;
1411         }
1412
1413         if (pcrd->eGeom == epullgCYL)
1414         {
1415             pull->bCylinder = TRUE;
1416         }
1417         /* We only need to calculate the plain COM of a group
1418          * when it is not only used as a cylinder group.
1419          */
1420         calc_com_start = (pcrd->eGeom == epullgCYL         ? 1 : 0);
1421         calc_com_end   = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2);
1422
1423         for (g = calc_com_start; g <= calc_com_end; g++)
1424         {
1425             pull->group[pcrd->group[g]].bCalcCOM = TRUE;
1426         }
1427     }
1428
1429     pull->ePBC = ir->ePBC;
1430     switch (pull->ePBC)
1431     {
1432         case epbcNONE: pull->npbcdim = 0; break;
1433         case epbcXY:   pull->npbcdim = 2; break;
1434         default:       pull->npbcdim = 3; break;
1435     }
1436
1437     if (fplog)
1438     {
1439         gmx_bool bAbs, bCos;
1440
1441         bAbs = FALSE;
1442         for (c = 0; c < pull->ncoord; c++)
1443         {
1444             if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1445                 pull->group[pull->coord[c].group[1]].nat == 0)
1446             {
1447                 bAbs = TRUE;
1448             }
1449         }
1450
1451         fprintf(fplog, "\n");
1452         if (pull->bPotential)
1453         {
1454             fprintf(fplog, "Will apply potential COM pulling\n");
1455         }
1456         if (pull->bConstraint)
1457         {
1458             fprintf(fplog, "Will apply constraint COM pulling\n");
1459         }
1460         fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1461                 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1462                 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1463         if (bAbs)
1464         {
1465             fprintf(fplog, "with an absolute reference\n");
1466         }
1467         bCos = FALSE;
1468         for (g = 0; g < pull->ngroup; g++)
1469         {
1470             if (pull->group[g].nat > 1 &&
1471                 pull->group[g].pbcatom < 0)
1472             {
1473                 /* We are using cosine weighting */
1474                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1475                 bCos = TRUE;
1476             }
1477         }
1478         if (bCos)
1479         {
1480             please_cite(fplog, "Engin2010");
1481         }
1482     }
1483
1484     pull->rbuf     = NULL;
1485     pull->dbuf     = NULL;
1486     pull->dbuf_cyl = NULL;
1487     pull->bRefAt   = FALSE;
1488     pull->cosdim   = -1;
1489     for (g = 0; g < pull->ngroup; g++)
1490     {
1491         pgrp           = &pull->group[g];
1492         pgrp->epgrppbc = epgrppbcNONE;
1493         if (pgrp->nat > 0)
1494         {
1495             /* There is an issue when a group is used in multiple coordinates
1496              * and constraints are applied in different dimensions with atoms
1497              * frozen in some, but not all dimensions.
1498              * Since there is only one mass array per group, we can't have
1499              * frozen/non-frozen atoms for different coords at the same time.
1500              * But since this is a very exotic case, we don't check for this.
1501              * A warning is printed in init_pull_group_index.
1502              */
1503
1504             gmx_bool bConstraint;
1505             ivec     pulldim, pulldim_con;
1506
1507             /* Loop over all pull coordinates to see along which dimensions
1508              * this group is pulled and if it is involved in constraints.
1509              */
1510             bConstraint = FALSE;
1511             clear_ivec(pulldim);
1512             clear_ivec(pulldim_con);
1513             for (c = 0; c < pull->ncoord; c++)
1514             {
1515                 if (pull->coord[c].group[0] == g ||
1516                     pull->coord[c].group[1] == g ||
1517                     (pull->coord[c].eGeom == epullgDIRRELATIVE &&
1518                      (pull->coord[c].group[2] == g ||
1519                       pull->coord[c].group[3] == g)))
1520                 {
1521                     for (m = 0; m < DIM; m++)
1522                     {
1523                         if (pull->coord[c].dim[m] == 1)
1524                         {
1525                             pulldim[m] = 1;
1526
1527                             if (pull->coord[c].eType == epullCONSTRAINT)
1528                             {
1529                                 bConstraint    = TRUE;
1530                                 pulldim_con[m] = 1;
1531                             }
1532                         }
1533                     }
1534                 }
1535             }
1536
1537             /* Determine if we need to take PBC into account for calculating
1538              * the COM's of the pull groups.
1539              */
1540             for (m = 0; m < pull->npbcdim; m++)
1541             {
1542                 if (pulldim[m] == 1 && pgrp->nat > 1)
1543                 {
1544                     if (pgrp->pbcatom >= 0)
1545                     {
1546                         pgrp->epgrppbc = epgrppbcREFAT;
1547                         pull->bRefAt   = TRUE;
1548                     }
1549                     else
1550                     {
1551                         if (pgrp->weight)
1552                         {
1553                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1554                         }
1555                         pgrp->epgrppbc = epgrppbcCOS;
1556                         if (pull->cosdim >= 0 && pull->cosdim != m)
1557                         {
1558                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
1559                         }
1560                         pull->cosdim = m;
1561                     }
1562                 }
1563             }
1564             /* Set the indices */
1565             init_pull_group_index(fplog, cr, g, pgrp,
1566                                   bConstraint, pulldim_con,
1567                                   mtop, ir, lambda);
1568         }
1569         else
1570         {
1571             /* Absolute reference, set the inverse mass to zero.
1572              * This is only relevant (and used) with constraint pulling.
1573              */
1574             pgrp->invtm  = 0;
1575             pgrp->wscale = 1;
1576         }
1577     }
1578
1579     /* If we use cylinder coordinates, do some initialising for them */
1580     if (pull->bCylinder)
1581     {
1582         snew(pull->dyna, pull->ncoord);
1583
1584         for (c = 0; c < pull->ncoord; c++)
1585         {
1586             const t_pull_coord *pcrd;
1587
1588             pcrd = &pull->coord[c];
1589
1590             if (pcrd->eGeom == epullgCYL)
1591             {
1592                 if (pull->group[pcrd->group[0]].nat == 0)
1593                 {
1594                     gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
1595                 }
1596             }
1597         }
1598     }
1599
1600     /* We still need to initialize the PBC reference coordinates */
1601     pull->bSetPBCatoms = TRUE;
1602
1603     /* Only do I/O when we are doing dynamics and if we are the MASTER */
1604     pull->out_x = NULL;
1605     pull->out_f = NULL;
1606     if (bOutFile)
1607     {
1608         if (pull->nstxout > 0)
1609         {
1610             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
1611                                         TRUE, Flags);
1612         }
1613         if (pull->nstfout > 0)
1614         {
1615             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1616                                         FALSE, Flags);
1617         }
1618     }
1619 }
1620
1621 void finish_pull(t_pull *pull)
1622 {
1623     if (pull->out_x)
1624     {
1625         gmx_fio_fclose(pull->out_x);
1626     }
1627     if (pull->out_f)
1628     {
1629         gmx_fio_fclose(pull->out_f);
1630     }
1631 }