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