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