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