src/mdlib/pull.c

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