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