src/gromacs/pulling/pull.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, 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 <math.h>
  40 #include <stdio.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include "gromacs/utility/futil.h"
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/legacyheaders/types/commrec.h"
  47 #include "gromacs/legacyheaders/network.h"
  48 #include "pull.h"
  49 #include "gromacs/legacyheaders/names.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/topology/mtop_util.h"
  52 #include "gromacs/legacyheaders/mdrun.h"
  53 #include "gromacs/legacyheaders/gmx_ga2la.h"
  54 #include "gromacs/legacyheaders/copyrite.h"
  55 #include "gromacs/legacyheaders/macros.h"
  56
  57 #include "gromacs/fileio/filenm.h"
  58 #include "gromacs/fileio/gmxfio.h"
  59 #include "gromacs/fileio/xvgr.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
  64 {
  65     int m;
  66
  67     for (m = 0; m < DIM; m++)
  68     {
  69         if (dim[m])
  70         {
  71             fprintf(out, "\t%g", pgrp->x[m]);
  72         }
  73     }
  74 }
  75
  76 static void pull_print_coord_dr(FILE *out, ivec dim, const t_pull_coord *pcrd)
  77 {
  78     int m;
  79
  80     for (m = 0; m < DIM; m++)
  81     {
  82         if (dim[m])
  83         {
  84             fprintf(out, "\t%g", pcrd->dr[m]);
  85         }
  86     }
  87 }
  88
  89 static void pull_print_x(FILE *out, t_pull *pull, double t)
  90 {
  91     int                 c;
  92     const t_pull_coord *pcrd;
  93
  94     fprintf(out, "%.4f", t);
  95
  96     for (c = 0; c < pull->ncoord; c++)
  97     {
  98         pcrd = &pull->coord[c];
  99
 100         if (pull->bPrintRef)
 101         {
 102             if (PULL_CYL(pull))
 103             {
 104                 pull_print_group_x(out, pull->dim, &pull->dyna[c]);
 105             }
 106             else
 107             {
 108                 pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]);
 109             }
 110         }
 111         pull_print_coord_dr(out, pull->dim, pcrd);
 112     }
 113     fprintf(out, "\n");
 114 }
 115
 116 static void pull_print_f(FILE *out, t_pull *pull, double t)
 117 {
 118     int c, d;
 119
 120     fprintf(out, "%.4f", t);
 121
 122     for (c = 0; c < pull->ncoord; c++)
 123     {
 124         fprintf(out, "\t%g", pull->coord[c].f_scal);
 125     }
 126     fprintf(out, "\n");
 127 }
 128
 129 void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
 130 {
 131     if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
 132     {
 133         pull_print_x(pull->out_x, pull, time);
 134     }
 135
 136     if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
 137     {
 138         pull_print_f(pull->out_f, pull, time);
 139     }
 140 }
 141
 142 static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
 143                            gmx_bool bCoord, unsigned long Flags)
 144 {
 145     FILE  *fp;
 146     int    nsets, c, m;
 147     char **setname, buf[10];
 148
 149     if (Flags & MD_APPENDFILES)
 150     {
 151         fp = gmx_fio_fopen(fn, "a+");
 152     }
 153     else
 154     {
 155         fp = gmx_fio_fopen(fn, "w+");
 156         if (bCoord)
 157         {
 158             xvgr_header(fp, "Pull COM",  "Time (ps)", "Position (nm)",
 159                         exvggtXNY, oenv);
 160         }
 161         else
 162         {
 163             xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
 164                         exvggtXNY, oenv);
 165         }
 166
 167         snew(setname, 2*pull->ncoord*DIM);
 168         nsets = 0;
 169         for (c = 0; c < pull->ncoord; c++)
 170         {
 171             if (bCoord)
 172             {
 173                 if (pull->bPrintRef)
 174                 {
 175                     for (m = 0; m < DIM; m++)
 176                     {
 177                         if (pull->dim[m])
 178                         {
 179                             sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
 180                             setname[nsets] = gmx_strdup(buf);
 181                             nsets++;
 182                         }
 183                     }
 184                 }
 185                 for (m = 0; m < DIM; m++)
 186                 {
 187                     if (pull->dim[m])
 188                     {
 189                         sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
 190                         setname[nsets] = gmx_strdup(buf);
 191                         nsets++;
 192                     }
 193                 }
 194             }
 195             else
 196             {
 197                 sprintf(buf, "%d", c+1);
 198                 setname[nsets] = gmx_strdup(buf);
 199                 nsets++;
 200             }
 201         }
 202         if (nsets > 1)
 203         {
 204             xvgr_legend(fp, nsets, (const char**)setname, oenv);
 205         }
 206         for (c = 0; c < nsets; c++)
 207         {
 208             sfree(setname[c]);
 209         }
 210         sfree(setname);
 211     }
 212
 213     return fp;
 214 }
 215
 216 /* Apply forces in a mass weighted fashion */
 217 static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
 218                              const dvec f_pull, int sign, rvec *f)
 219 {
 220     int    i, ii, m;
 221     double wmass, inv_wm;
 222
 223     inv_wm = pgrp->wscale*pgrp->invtm;
 224
 225     for (i = 0; i < pgrp->nat_loc; i++)
 226     {
 227         ii    = pgrp->ind_loc[i];
 228         wmass = md->massT[ii];
 229         if (pgrp->weight_loc)
 230         {
 231             wmass *= pgrp->weight_loc[i];
 232         }
 233
 234         for (m = 0; m < DIM; m++)
 235         {
 236             f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
 237         }
 238     }
 239 }
 240
 241 /* Apply forces in a mass weighted fashion */
 242 static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
 243 {
 244     int                 c;
 245     const t_pull_coord *pcrd;
 246
 247     for (c = 0; c < pull->ncoord; c++)
 248     {
 249         pcrd = &pull->coord[c];
 250
 251         if (PULL_CYL(pull))
 252         {
 253             apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f);
 254         }
 255         else
 256         {
 257             if (pull->group[pcrd->group[0]].nat > 0)
 258             {
 259                 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
 260             }
 261         }
 262         apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
 263     }
 264 }
 265
 266 static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc)
 267 {
 268     double max_d2;
 269     int    m;
 270
 271     max_d2 = GMX_DOUBLE_MAX;
 272
 273     if (pull->eGeom != epullgDIRPBC)
 274     {
 275         for (m = 0; m < pbc->ndim_ePBC; m++)
 276         {
 277             if (pull->dim[m] != 0)
 278             {
 279                 max_d2 = min(max_d2, norm2(pbc->box[m]));
 280             }
 281         }
 282     }
 283
 284     return 0.25*max_d2;
 285 }
 286
 287 static void low_get_pull_coord_dr(const t_pull *pull,
 288                                   const t_pull_coord *pcrd,
 289                                   const t_pbc *pbc, double t,
 290                                   dvec xg, dvec xref, double max_dist2,
 291                                   dvec dr)
 292 {
 293     const t_pull_group *pgrp0, *pgrp1;
 294     int                 m;
 295     dvec                xrefr, dref = {0, 0, 0};
 296     double              dr2;
 297
 298     pgrp0 = &pull->group[pcrd->group[0]];
 299     pgrp1 = &pull->group[pcrd->group[1]];
 300
 301     /* Only the first group can be an absolute reference, in that case nat=0 */
 302     if (pgrp0->nat == 0)
 303     {
 304         for (m = 0; m < DIM; m++)
 305         {
 306             xref[m] = pcrd->origin[m];
 307         }
 308     }
 309
 310     copy_dvec(xref, xrefr);
 311
 312     if (pull->eGeom == epullgDIRPBC)
 313     {
 314         for (m = 0; m < DIM; m++)
 315         {
 316             dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
 317         }
 318         /* Add the reference position, so we use the correct periodic image */
 319         dvec_inc(xrefr, dref);
 320     }
 321
 322     pbc_dx_d(pbc, xg, xrefr, dr);
 323     dr2 = 0;
 324     for (m = 0; m < DIM; m++)
 325     {
 326         dr[m] *= pull->dim[m];
 327         dr2   += dr[m]*dr[m];
 328     }
 329     if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
 330     {
 331         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",
 332                   pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
 333     }
 334
 335     if (pull->eGeom == epullgDIRPBC)
 336     {
 337         dvec_inc(dr, dref);
 338     }
 339 }
 340
 341 static void get_pull_coord_dr(const t_pull *pull,
 342                               int coord_ind,
 343                               const t_pbc *pbc, double t,
 344                               dvec dr)
 345 {
 346     double              md2;
 347     const t_pull_coord *pcrd;
 348
 349     if (pull->eGeom == epullgDIRPBC)
 350     {
 351         md2 = -1;
 352     }
 353     else
 354     {
 355         md2 = max_pull_distance2(pull, pbc);
 356     }
 357
 358     pcrd = &pull->coord[coord_ind];
 359
 360     low_get_pull_coord_dr(pull, pcrd, pbc, t,
 361                           pull->group[pcrd->group[1]].x,
 362                           PULL_CYL(pull) ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
 363                           md2,
 364                           dr);
 365 }
 366
 367 void get_pull_coord_distance(const t_pull *pull,
 368                              int coord_ind,
 369                              const t_pbc *pbc, double t,
 370                              dvec dr, double *dev)
 371 {
 372     static gmx_bool     bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
 373                                             but is fairly benign */
 374     const t_pull_coord *pcrd;
 375     int                 m;
 376     double              ref, drs, inpr;
 377
 378     pcrd = &pull->coord[coord_ind];
 379
 380     get_pull_coord_dr(pull, coord_ind, pbc, t, dr);
 381
 382     ref = pcrd->init + pcrd->rate*t;
 383
 384     switch (pull->eGeom)
 385     {
 386         case epullgDIST:
 387             /* Pull along the vector between the com's */
 388             if (ref < 0 && !bWarned)
 389             {
 390                 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
 391                 bWarned = TRUE;
 392             }
 393             drs = dnorm(dr);
 394             if (drs == 0)
 395             {
 396                 /* With no vector we can not determine the direction for the force,
 397                  * so we set the force to zero.
 398                  */
 399                 *dev = 0;
 400             }
 401             else
 402             {
 403                 /* Determine the deviation */
 404                 *dev = drs - ref;
 405             }
 406             break;
 407         case epullgDIR:
 408         case epullgDIRPBC:
 409         case epullgCYL:
 410             /* Pull along vec */
 411             inpr = 0;
 412             for (m = 0; m < DIM; m++)
 413             {
 414                 inpr += pcrd->vec[m]*dr[m];
 415             }
 416             *dev = inpr - ref;
 417             break;
 418     }
 419 }
 420
 421 void clear_pull_forces(t_pull *pull)
 422 {
 423     int i;
 424
 425     /* Zeroing the forces is only required for constraint pulling.
 426      * It can happen that multiple constraint steps need to be applied
 427      * and therefore the constraint forces need to be accumulated.
 428      */
 429     for (i = 0; i < pull->ncoord; i++)
 430     {
 431         clear_dvec(pull->coord[i].f);
 432         pull->coord[i].f_scal = 0;
 433     }
 434 }
 435
 436 /* Apply constraint using SHAKE */
 437 static void do_constraint(t_pull *pull, t_pbc *pbc,
 438                           rvec *x, rvec *v,
 439                           gmx_bool bMaster, tensor vir,
 440                           double dt, double t)
 441 {
 442
 443     dvec         *r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
 444     dvec          unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
 445     dvec         *rnew;   /* current 'new' positions of the groups */
 446     double       *dr_tot; /* the total update of the coords */
 447     double        ref;
 448     dvec          vec;
 449     double        d0, inpr;
 450     double        lambda, rm, mass, invdt = 0;
 451     gmx_bool      bConverged_all, bConverged = FALSE;
 452     int           niter = 0, g, c, ii, j, m, max_iter = 100;
 453     double        a;
 454     dvec          f;       /* the pull force */
 455     dvec          tmp, tmp3;
 456     t_pull_group *pdyna, *pgrp0, *pgrp1;
 457     t_pull_coord *pcrd;
 458
 459     snew(r_ij,   pull->ncoord);
 460     snew(dr_tot, pull->ncoord);
 461
 462     snew(rnew, pull->ngroup);
 463
 464     /* copy the current unconstrained positions for use in iterations. We
 465        iterate until rinew[i] and rjnew[j] obey the constraints. Then
 466        rinew - pull.x_unc[i] is the correction dr to group i */
 467     for (g = 0; g < pull->ngroup; g++)
 468     {
 469         copy_dvec(pull->group[g].xp, rnew[g]);
 470     }
 471     if (PULL_CYL(pull))
 472     {
 473         /* There is only one pull coordinate and reference group */
 474         copy_dvec(pull->dyna[0].xp, rnew[pull->coord[0].group[0]]);
 475     }
 476
 477     /* Determine the constraint directions from the old positions */
 478     for (c = 0; c < pull->ncoord; c++)
 479     {
 480         get_pull_coord_dr(pull, c, pbc, t, r_ij[c]);
 481         /* Store the difference vector at time t for printing */
 482         copy_dvec(r_ij[c], pull->coord[c].dr);
 483         if (debug)
 484         {
 485             fprintf(debug, "Pull coord %d dr %f %f %f\n",
 486                     c, r_ij[c][XX], r_ij[c][YY], r_ij[c][ZZ]);
 487         }
 488
 489         if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
 490         {
 491             /* Select the component along vec */
 492             a = 0;
 493             for (m = 0; m < DIM; m++)
 494             {
 495                 a += pull->coord[c].vec[m]*r_ij[c][m];
 496             }
 497             for (m = 0; m < DIM; m++)
 498             {
 499                 r_ij[c][m] = a*pull->coord[c].vec[m];
 500             }
 501         }
 502
 503         if (dnorm2(r_ij[c]) == 0)
 504         {
 505             gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
 506         }
 507     }
 508
 509     bConverged_all = FALSE;
 510     while (!bConverged_all && niter < max_iter)
 511     {
 512         bConverged_all = TRUE;
 513
 514         /* loop over all constraints */
 515         for (c = 0; c < pull->ncoord; c++)
 516         {
 517             dvec dr0, dr1;
 518
 519             pcrd  = &pull->coord[c];
 520             pgrp0 = &pull->group[pcrd->group[0]];
 521             pgrp1 = &pull->group[pcrd->group[1]];
 522
 523             /* Get the current difference vector */
 524             low_get_pull_coord_dr(pull, pcrd, pbc, t,
 525                                   rnew[pcrd->group[1]],
 526                                   rnew[pcrd->group[0]],
 527                                   -1, unc_ij);
 528
 529             ref = pcrd->init + pcrd->rate*t;
 530
 531             if (debug)
 532             {
 533                 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
 534             }
 535
 536             rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
 537
 538             switch (pull->eGeom)
 539             {
 540                 case epullgDIST:
 541                     if (ref <= 0)
 542                     {
 543                         gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
 544                     }
 545
 546                     {
 547                         double q, c_a, c_b, c_c;
 548
 549                         c_a = diprod(r_ij[c], r_ij[c]);
 550                         c_b = diprod(unc_ij, r_ij[c])*2;
 551                         c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
 552
 553                         if (c_b < 0)
 554                         {
 555                             q      = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
 556                             lambda = -q/c_a;
 557                         }
 558                         else
 559                         {
 560                             q      = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
 561                             lambda = -c_c/q;
 562                         }
 563
 564                         if (debug)
 565                         {
 566                             fprintf(debug,
 567                                     "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
 568                                     c_a, c_b, c_c, lambda);
 569                         }
 570                     }
 571
 572                     /* The position corrections dr due to the constraints */
 573                     dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
 574                     dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
 575                     dr_tot[c] += -lambda*dnorm(r_ij[c]);
 576                     break;
 577                 case epullgDIR:
 578                 case epullgDIRPBC:
 579                 case epullgCYL:
 580                     /* A 1-dimensional constraint along a vector */
 581                     a = 0;
 582                     for (m = 0; m < DIM; m++)
 583                     {
 584                         vec[m] = pcrd->vec[m];
 585                         a     += unc_ij[m]*vec[m];
 586                     }
 587                     /* Select only the component along the vector */
 588                     dsvmul(a, vec, unc_ij);
 589                     lambda = a - ref;
 590                     if (debug)
 591                     {
 592                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
 593                     }
 594
 595                     /* The position corrections dr due to the constraints */
 596                     dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
 597                     dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
 598                     dr_tot[c] += -lambda;
 599                     break;
 600             }
 601
 602             /* DEBUG */
 603             if (debug)
 604             {
 605                 int g0, g1;
 606
 607                 g0 = pcrd->group[0];
 608                 g1 = pcrd->group[1];
 609                 low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
 610                 low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
 611                 fprintf(debug,
 612                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 613                         rnew[g0][0], rnew[g0][1], rnew[g0][2],
 614                         rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
 615                 fprintf(debug,
 616                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
 617                         "", "", "", "", "", "", ref);
 618                 fprintf(debug,
 619                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
 620                         dr0[0], dr0[1], dr0[2],
 621                         dr1[0], dr1[1], dr1[2],
 622                         dnorm(tmp3));
 623             } /* END DEBUG */
 624
 625             /* Update the COMs with dr */
 626             dvec_inc(rnew[pcrd->group[1]], dr1);
 627             dvec_inc(rnew[pcrd->group[0]], dr0);
 628         }
 629
 630         /* Check if all constraints are fullfilled now */
 631         for (c = 0; c < pull->ncoord; c++)
 632         {
 633             pcrd = &pull->coord[c];
 634             ref  = pcrd->init + pcrd->rate*t;
 635
 636             low_get_pull_coord_dr(pull, pcrd, pbc, t,
 637                                   rnew[pcrd->group[1]],
 638                                   rnew[pcrd->group[0]],
 639                                   -1, unc_ij);
 640
 641             switch (pull->eGeom)
 642             {
 643                 case epullgDIST:
 644                     bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
 645                     break;
 646                 case epullgDIR:
 647                 case epullgDIRPBC:
 648                 case epullgCYL:
 649                     for (m = 0; m < DIM; m++)
 650                     {
 651                         vec[m] = pcrd->vec[m];
 652                     }
 653                     inpr = diprod(unc_ij, vec);
 654                     dsvmul(inpr, vec, unc_ij);
 655                     bConverged =
 656                         fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
 657                     break;
 658             }
 659
 660             if (!bConverged)
 661             {
 662                 if (debug)
 663                 {
 664                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
 665                             "d_ref = %f, current d = %f\n",
 666                             g, ref, dnorm(unc_ij));
 667                 }
 668
 669                 bConverged_all = FALSE;
 670             }
 671         }
 672
 673         niter++;
 674         /* if after all constraints are dealt with and bConverged is still TRUE
 675            we're finished, if not we do another iteration */
 676     }
 677     if (niter > max_iter)
 678     {
 679         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
 680     }
 681
 682     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 683
 684     if (v)
 685     {
 686         invdt = 1/dt;
 687     }
 688
 689     /* update atoms in the groups */
 690     for (g = 0; g < pull->ngroup; g++)
 691     {
 692         const t_pull_group *pgrp;
 693         dvec                dr;
 694
 695         if (PULL_CYL(pull) && g == pull->coord[0].group[0])
 696         {
 697             pgrp = &pull->dyna[0];
 698         }
 699         else
 700         {
 701             pgrp = &pull->group[g];
 702         }
 703
 704         /* get the final constraint displacement dr for group g */
 705         dvec_sub(rnew[g], pgrp->xp, dr);
 706         /* select components of dr */
 707         for (m = 0; m < DIM; m++)
 708         {
 709             dr[m] *= pull->dim[m];
 710         }
 711
 712         /* update the atom positions */
 713         copy_dvec(dr, tmp);
 714         for (j = 0; j < pgrp->nat_loc; j++)
 715         {
 716             ii = pgrp->ind_loc[j];
 717             if (pgrp->weight_loc)
 718             {
 719                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
 720             }
 721             for (m = 0; m < DIM; m++)
 722             {
 723                 x[ii][m] += tmp[m];
 724             }
 725             if (v)
 726             {
 727                 for (m = 0; m < DIM; m++)
 728                 {
 729                     v[ii][m] += invdt*tmp[m];
 730                 }
 731             }
 732         }
 733     }
 734
 735     /* calculate the constraint forces, used for output and virial only */
 736     for (c = 0; c < pull->ncoord; c++)
 737     {
 738         pcrd         = &pull->coord[c];
 739         pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
 740
 741         if (vir && bMaster)
 742         {
 743             double f_invr;
 744
 745             /* Add the pull contribution to the virial */
 746             /* We have already checked above that r_ij[c] != 0 */
 747             f_invr = pcrd->f_scal/dnorm(r_ij[c]);
 748
 749             for (j = 0; j < DIM; j++)
 750             {
 751                 for (m = 0; m < DIM; m++)
 752                 {
 753                     vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
 754                 }
 755             }
 756         }
 757     }
 758
 759     /* finished! I hope. Give back some memory */
 760     sfree(r_ij);
 761     sfree(dr_tot);
 762     sfree(rnew);
 763 }
 764
 765 /* Pulling with a harmonic umbrella potential or constant force */
 766 static void do_pull_pot(int ePull,
 767                         t_pull *pull, t_pbc *pbc, double t, real lambda,
 768                         real *V, tensor vir, real *dVdl)
 769 {
 770     int           c, j, m;
 771     double        dev, ndr, invdr;
 772     real          k, dkdl;
 773     t_pull_coord *pcrd;
 774
 775     /* loop over the pull coordinates */
 776     *V    = 0;
 777     *dVdl = 0;
 778     for (c = 0; c < pull->ncoord; c++)
 779     {
 780         pcrd = &pull->coord[c];
 781
 782         get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev);
 783
 784         k    = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
 785         dkdl = pcrd->kB - pcrd->k;
 786
 787         switch (pull->eGeom)
 788         {
 789             case epullgDIST:
 790                 ndr   = dnorm(pcrd->dr);
 791                 if (ndr > 0)
 792                 {
 793                     invdr = 1/ndr;
 794                 }
 795                 else
 796                 {
 797                     /* With an harmonic umbrella, the force is 0 at r=0,
 798                      * so we can set invdr to any value.
 799                      * With a constant force, the force at r=0 is not defined,
 800                      * so we zero it (this is anyhow a very rare event).
 801                      */
 802                     invdr = 0;
 803                 }
 804                 if (ePull == epullUMBRELLA)
 805                 {
 806                     pcrd->f_scal  =       -k*dev;
 807                     *V           += 0.5*   k*dsqr(dev);
 808                     *dVdl        += 0.5*dkdl*dsqr(dev);
 809                 }
 810                 else
 811                 {
 812                     pcrd->f_scal  =   -k;
 813                     *V           +=    k*ndr;
 814                     *dVdl        += dkdl*ndr;
 815                 }
 816                 for (m = 0; m < DIM; m++)
 817                 {
 818                     pcrd->f[m]    = pcrd->f_scal*pcrd->dr[m]*invdr;
 819                 }
 820                 break;
 821             case epullgDIR:
 822             case epullgDIRPBC:
 823             case epullgCYL:
 824                 if (ePull == epullUMBRELLA)
 825                 {
 826                     pcrd->f_scal  =       -k*dev;
 827                     *V           += 0.5*   k*dsqr(dev);
 828                     *dVdl        += 0.5*dkdl*dsqr(dev);
 829                 }
 830                 else
 831                 {
 832                     ndr = 0;
 833                     for (m = 0; m < DIM; m++)
 834                     {
 835                         ndr += pcrd->vec[m]*pcrd->dr[m];
 836                     }
 837                     pcrd->f_scal  =   -k;
 838                     *V           +=    k*ndr;
 839                     *dVdl        += dkdl*ndr;
 840                 }
 841                 for (m = 0; m < DIM; m++)
 842                 {
 843                     pcrd->f[m]    = pcrd->f_scal*pcrd->vec[m];
 844                 }
 845                 break;
 846         }
 847
 848         if (vir)
 849         {
 850             /* Add the pull contribution to the virial */
 851             for (j = 0; j < DIM; j++)
 852             {
 853                 for (m = 0; m < DIM; m++)
 854                 {
 855                     vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
 856                 }
 857             }
 858         }
 859     }
 860 }
 861
 862 real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
 863                     t_commrec *cr, double t, real lambda,
 864                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
 865 {
 866     real V, dVdl;
 867
 868     pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
 869
 870     do_pull_pot(ePull, pull, pbc, t, lambda,
 871                 &V, pull->bVirial && MASTER(cr) ? vir : NULL, &dVdl);
 872
 873     /* Distribute forces over pulled groups */
 874     apply_forces(pull, md, f);
 875
 876     if (MASTER(cr))
 877     {
 878         *dvdlambda += dVdl;
 879     }
 880
 881     return (MASTER(cr) ? V : 0.0);
 882 }
 883
 884 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
 885                      t_commrec *cr, double dt, double t,
 886                      rvec *x, rvec *xp, rvec *v, tensor vir)
 887 {
 888     pull_calc_coms(cr, pull, md, pbc, t, x, xp);
 889
 890     do_constraint(pull, pbc, xp, v, pull->bVirial && MASTER(cr), vir, dt, t);
 891 }
 892
 893 static void make_local_pull_group(gmx_ga2la_t ga2la,
 894                                   t_pull_group *pg, int start, int end)
 895 {
 896     int i, ii;
 897
 898     pg->nat_loc = 0;
 899     for (i = 0; i < pg->nat; i++)
 900     {
 901         ii = pg->ind[i];
 902         if (ga2la)
 903         {
 904             if (!ga2la_get_home(ga2la, ii, &ii))
 905             {
 906                 ii = -1;
 907             }
 908         }
 909         if (ii >= start && ii < end)
 910         {
 911             /* This is a home atom, add it to the local pull group */
 912             if (pg->nat_loc >= pg->nalloc_loc)
 913             {
 914                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
 915                 srenew(pg->ind_loc, pg->nalloc_loc);
 916                 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
 917                 {
 918                     srenew(pg->weight_loc, pg->nalloc_loc);
 919                 }
 920             }
 921             pg->ind_loc[pg->nat_loc] = ii;
 922             if (pg->weight)
 923             {
 924                 pg->weight_loc[pg->nat_loc] = pg->weight[i];
 925             }
 926             pg->nat_loc++;
 927         }
 928     }
 929 }
 930
 931 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
 932 {
 933     gmx_ga2la_t ga2la;
 934     int         g;
 935
 936     if (dd)
 937     {
 938         ga2la = dd->ga2la;
 939     }
 940     else
 941     {
 942         ga2la = NULL;
 943     }
 944
 945     for (g = 0; g < pull->ngroup; g++)
 946     {
 947         make_local_pull_group(ga2la, &pull->group[g],
 948                               0, md->homenr);
 949     }
 950
 951     /* Since the PBC of atoms might have changed, we need to update the PBC */
 952     pull->bSetPBCatoms = TRUE;
 953 }
 954
 955 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
 956                                   int g, t_pull_group *pg, ivec pulldims,
 957                                   gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
 958 {
 959     int                   i, ii, d, nfrozen, ndim;
 960     real                  m, w, mbd;
 961     double                tmass, wmass, wwmass;
 962     gmx_groups_t         *groups;
 963     gmx_mtop_atomlookup_t alook;
 964     t_atom               *atom;
 965
 966     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
 967     {
 968         /* There are no masses in the integrator.
 969          * But we still want to have the correct mass-weighted COMs.
 970          * So we store the real masses in the weights.
 971          * We do not set nweight, so these weights do not end up in the tpx file.
 972          */
 973         if (pg->nweight == 0)
 974         {
 975             snew(pg->weight, pg->nat);
 976         }
 977     }
 978
 979     if (cr && PAR(cr))
 980     {
 981         pg->nat_loc    = 0;
 982         pg->nalloc_loc = 0;
 983         pg->ind_loc    = NULL;
 984         pg->weight_loc = NULL;
 985     }
 986     else
 987     {
 988         pg->nat_loc = pg->nat;
 989         pg->ind_loc = pg->ind;
 990         if (pg->epgrppbc == epgrppbcCOS)
 991         {
 992             snew(pg->weight_loc, pg->nat);
 993         }
 994         else
 995         {
 996             pg->weight_loc = pg->weight;
 997         }
 998     }
 999
1000     groups = &mtop->groups;
1001
1002     alook = gmx_mtop_atomlookup_init(mtop);
1003
1004     nfrozen = 0;
1005     tmass   = 0;
1006     wmass   = 0;
1007     wwmass  = 0;
1008     for (i = 0; i < pg->nat; i++)
1009     {
1010         ii = pg->ind[i];
1011         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1012         if (ir->opts.nFreeze)
1013         {
1014             for (d = 0; d < DIM; d++)
1015             {
1016                 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1017                 {
1018                     nfrozen++;
1019                 }
1020             }
1021         }
1022         if (ir->efep == efepNO)
1023         {
1024             m = atom->m;
1025         }
1026         else
1027         {
1028             m = (1 - lambda)*atom->m + lambda*atom->mB;
1029         }
1030         if (pg->nweight > 0)
1031         {
1032             w = pg->weight[i];
1033         }
1034         else
1035         {
1036             w = 1;
1037         }
1038         if (EI_ENERGY_MINIMIZATION(ir->eI))
1039         {
1040             /* Move the mass to the weight */
1041             w            *= m;
1042             m             = 1;
1043             pg->weight[i] = w;
1044         }
1045         else if (ir->eI == eiBD)
1046         {
1047             if (ir->bd_fric)
1048             {
1049                 mbd = ir->bd_fric*ir->delta_t;
1050             }
1051             else
1052             {
1053                 if (groups->grpnr[egcTC] == NULL)
1054                 {
1055                     mbd = ir->delta_t/ir->opts.tau_t[0];
1056                 }
1057                 else
1058                 {
1059                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1060                 }
1061             }
1062             w            *= m/mbd;
1063             m             = mbd;
1064             pg->weight[i] = w;
1065         }
1066         tmass  += m;
1067         wmass  += m*w;
1068         wwmass += m*w*w;
1069     }
1070
1071     gmx_mtop_atomlookup_destroy(alook);
1072
1073     if (wmass == 0)
1074     {
1075         gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1076                   pg->weight ? " weighted" : "", g);
1077     }
1078     if (fplog)
1079     {
1080         fprintf(fplog,
1081                 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1082         if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1083         {
1084             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1085         }
1086         if (pg->epgrppbc == epgrppbcCOS)
1087         {
1088             fprintf(fplog, ", cosine weighting will be used");
1089         }
1090         fprintf(fplog, "\n");
1091     }
1092
1093     if (nfrozen == 0)
1094     {
1095         /* A value > 0 signals not frozen, it is updated later */
1096         pg->invtm  = 1.0;
1097     }
1098     else
1099     {
1100         ndim = 0;
1101         for (d = 0; d < DIM; d++)
1102         {
1103             ndim += pulldims[d]*pg->nat;
1104         }
1105         if (fplog && nfrozen > 0 && nfrozen < ndim)
1106         {
1107             fprintf(fplog,
1108                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1109                     "         that are subject to pulling are frozen.\n"
1110                     "         For pulling the whole group will be frozen.\n\n",
1111                     g);
1112         }
1113         pg->invtm  = 0.0;
1114         pg->wscale = 1.0;
1115     }
1116 }
1117
1118 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1119                gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1120                gmx_bool bOutFile, unsigned long Flags)
1121 {
1122     t_pull       *pull;
1123     t_pull_group *pgrp;
1124     int           c, g, start = 0, end = 0, m;
1125
1126     pull = ir->pull;
1127
1128     pull->ePBC = ir->ePBC;
1129     switch (pull->ePBC)
1130     {
1131         case epbcNONE: pull->npbcdim = 0; break;
1132         case epbcXY:   pull->npbcdim = 2; break;
1133         default:       pull->npbcdim = 3; break;
1134     }
1135
1136     if (fplog)
1137     {
1138         gmx_bool bAbs, bCos;
1139
1140         bAbs = FALSE;
1141         for (c = 0; c < pull->ncoord; c++)
1142         {
1143             if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1144                 pull->group[pull->coord[c].group[1]].nat == 0)
1145             {
1146                 bAbs = TRUE;
1147             }
1148         }
1149
1150         fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
1151                 EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
1152         fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1153                 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1154                 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1155         if (bAbs)
1156         {
1157             fprintf(fplog, "with an absolute reference\n");
1158         }
1159         bCos = FALSE;
1160         for (g = 0; g < pull->ngroup; g++)
1161         {
1162             if (pull->group[g].nat > 1 &&
1163                 pull->group[g].pbcatom < 0)
1164             {
1165                 /* We are using cosine weighting */
1166                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1167                 bCos = TRUE;
1168             }
1169         }
1170         if (bCos)
1171         {
1172             please_cite(fplog, "Engin2010");
1173         }
1174     }
1175
1176     /* We always add the virial contribution,
1177      * except for geometry = direction_periodic where this is impossible.
1178      */
1179     pull->bVirial = (pull->eGeom != epullgDIRPBC);
1180     if (getenv("GMX_NO_PULLVIR") != NULL)
1181     {
1182         if (fplog)
1183         {
1184             fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
1185         }
1186         pull->bVirial = FALSE;
1187     }
1188
1189     pull->rbuf     = NULL;
1190     pull->dbuf     = NULL;
1191     pull->dbuf_cyl = NULL;
1192     pull->bRefAt   = FALSE;
1193     pull->cosdim   = -1;
1194     for (g = 0; g < pull->ngroup; g++)
1195     {
1196         pgrp           = &pull->group[g];
1197         pgrp->epgrppbc = epgrppbcNONE;
1198         if (pgrp->nat > 0)
1199         {
1200             /* Determine if we need to take PBC into account for calculating
1201              * the COM's of the pull groups.
1202              */
1203             for (m = 0; m < pull->npbcdim; m++)
1204             {
1205                 if (pull->dim[m] && pgrp->nat > 1)
1206                 {
1207                     if (pgrp->pbcatom >= 0)
1208                     {
1209                         pgrp->epgrppbc = epgrppbcREFAT;
1210                         pull->bRefAt   = TRUE;
1211                     }
1212                     else
1213                     {
1214                         if (pgrp->weight)
1215                         {
1216                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1217                         }
1218                         pgrp->epgrppbc = epgrppbcCOS;
1219                         if (pull->cosdim >= 0 && pull->cosdim != m)
1220                         {
1221                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1222                         }
1223                         pull->cosdim = m;
1224                     }
1225                 }
1226             }
1227             /* Set the indices */
1228             init_pull_group_index(fplog, cr, g, pgrp, pull->dim, mtop, ir, lambda);
1229             if (PULL_CYL(pull) && pgrp->invtm == 0)
1230             {
1231                 gmx_fatal(FARGS, "Can not have frozen atoms in a cylinder pull group");
1232             }
1233         }
1234         else
1235         {
1236             /* Absolute reference, set the inverse mass to zero */
1237             pgrp->invtm  = 0;
1238             pgrp->wscale = 1;
1239         }
1240     }
1241
1242     /* if we use dynamic reference groups, do some initialising for them */
1243     if (PULL_CYL(pull))
1244     {
1245         if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1)
1246         {
1247             /* We can't easily update the single reference group with multiple
1248              * constraints. This would require recalculating COMs.
1249              */
1250             gmx_fatal(FARGS, "Constraint COM pulling supports only one coordinate with geometry=cylinder, you can use umbrella pulling with multiple coordinates");
1251         }
1252
1253         for (c = 0; c < pull->ncoord; c++)
1254         {
1255             if (pull->group[pull->coord[c].group[0]].nat == 0)
1256             {
1257                 gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
1258             }
1259         }
1260
1261         snew(pull->dyna, pull->ncoord);
1262     }
1263
1264     /* We still need to initialize the PBC reference coordinates */
1265     pull->bSetPBCatoms = TRUE;
1266
1267     /* Only do I/O when we are doing dynamics and if we are the MASTER */
1268     pull->out_x = NULL;
1269     pull->out_f = NULL;
1270     if (bOutFile)
1271     {
1272         if (pull->nstxout > 0)
1273         {
1274             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags);
1275         }
1276         if (pull->nstfout > 0)
1277         {
1278             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1279                                         FALSE, Flags);
1280         }
1281     }
1282 }
1283
1284 void finish_pull(t_pull *pull)
1285 {
1286     if (pull->out_x)
1287     {
1288         gmx_fio_fclose(pull->out_x);
1289     }
1290     if (pull->out_f)
1291     {
1292         gmx_fio_fclose(pull->out_f);
1293     }
1294 }