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