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