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