src/mdlib/pull.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_pull_c = "$Id$";
  30
  31 #include "futil.h"
  32 #include "rdgroup.h"
  33 #include "statutil.h"
  34 #include "math.h"
  35 #include "stdio.h"
  36 #include "stdlib.h"
  37 #include "vec.h"
  38 #include "typedefs.h"
  39 #include "network.h"
  40 #include "filenm.h"
  41 #include "string.h"
  42 #include "smalloc.h"
  43 #include "pull.h"
  44
  45 /* check if we're close enough to the target position */
  46 static bool check_convergence(t_pull *pull) {
  47   bool bTest = TRUE;
  48   real tol;
  49   int i,m;
  50   rvec tmp1;
  51   rvec dr; /* distance between target and current position */
  52
  53   tol = pull->tolerance;
  54
  55   for (i=0;i<pull->pull.n;i++) {
  56     rvec_sub(pull->ref.x_unc[0],pull->pull.x_unc[i],tmp1);
  57     rvec_add(pull->pull.xtarget[i],tmp1,dr);
  58     /* multiply by pull->dims to pick the elements we want */
  59     for(m=0;m<DIM;m++)
  60       dr[m] = pull->dims[m]*dr[m];
  61     bTest = bTest && (norm(dr) < tol);
  62   }
  63
  64   return bTest;
  65 }
  66
  67 static void do_start(t_pull *pull, rvec *x, matrix box, t_mdatoms *md,
  68                      real dt, int step, t_topology *top)
  69 {
  70   int i,j,ii,m;
  71   rvec dr,dx,tmp;
  72   bool bThereYet,bDump;
  73   static int nout = 0;
  74   rvec ds;
  75
  76   /* pull.xtarget[i] is the desired position of group i with respect
  77      to the reference group */
  78
  79   /* bThereYet is true if all groups have reached their desired position */
  80   bThereYet = check_convergence(pull);
  81
  82   if (pull->bVerbose) {
  83     for (i=0;i<pull->pull.n;i++) {
  84       copy_rvec(pull->pull.xtarget[i],tmp);
  85       fprintf(stderr,"Group %d goal:%8.3f%8.3f%8.3f from reference\n",
  86               i,tmp[0],tmp[1],tmp[2]);
  87     }
  88   }
  89
  90   /* some groups are not there yet! */
  91   if (!bThereYet) {
  92     for (i=0;i<pull->pull.n;i++) {
  93       /* move this group closer to its target position */
  94       rvec_sub(pull->ref.x_unc[0],pull->pull.x_unc[i],tmp);
  95       rvec_add(pull->pull.xtarget[i],tmp,dr);
  96
  97       /* multiply by pull->dims to pick the elements we want */
  98       for(m=0;m<DIM;m++)
  99         dr[m] = pull->dims[m]*dr[m];
 100
 101       /* move over pull->xlt_rate in the direction of dr */
 102       svmul(pull->xlt_rate/norm(dr),dr,dx);
 103
 104       if (pull->bVerbose)
 105         fprintf(stderr,"To go:%10.2e %10.2e %10.2e\n",
 106                 dr[XX],dr[YY],dr[ZZ]);
 107
 108       /* update the actual positions */
 109       for (j=0;j<pull->pull.ngx[i];j++) {
 110         ii = pull->pull.idx[i][j];
 111         rvec_add(x[ii],dx,x[ii]);
 112       }
 113
 114       /* get the new center of mass in pull.x_unc[i] */
 115       (void)calc_com(x,pull->pull.ngx[i],pull->pull.idx[i],
 116                      md,pull->pull.x_unc[i],box);
 117     }
 118   }
 119
 120   bDump = check_convergence(pull);
 121
 122   if (bDump) {
 123     /* we have to write a coordinate file and reset the desired position */
 124     for (i=0;i<pull->pull.n;i++) {
 125       for (m=0;m<DIM;m++) {
 126         pull->pull.xtarget[i][m] += pull->pull.dir[i][m]*pull->xlt_incr/
 127           norm(pull->pull.dir[i]);
 128       }
 129
 130       if (pull->bVerbose)
 131         fprintf(stderr,"New target position: %8.3f%8.3f%8.3f\n",
 132                 pull->pull.xtarget[i][0],pull->pull.xtarget[i][1],
 133                 pull->pull.xtarget[i][2]);
 134     }
 135     dump_conf(pull,x,box,top,nout,step*dt/1000);
 136     nout++;
 137   }
 138 }
 139
 140 static void do_umbrella(t_pull *pull, rvec *x,rvec *f,matrix box,
 141                         t_mdatoms *md)
 142 {
 143   int i,j,m,g;
 144   real mi;
 145   rvec cm;    /* center of mass displacement of reference */
 146   rvec dr;    /* distance from com to potential center */
 147
 148   /* loop over the groups that are being umbrella sampled */
 149   for (i=0;i<pull->pull.n;i++) {
 150
 151     /* pull->dims selects which components (x,y,z) we want */
 152     for (m=0;m<DIM;m++) {
 153       dr[m]=pull->dims[m]*(pull->pull.x_ref[i][m]-pull->pull.x_unc[i][m]);
 154       if (dr[m] >  box[m][m]/2) dr[m]-=box[m][m];
 155       if (dr[m] < -box[m][m]/2) dr[m]+=box[m][m];
 156     }
 157
 158     /* f = um_width*x */
 159     svmul(pull->um_width,dr,pull->pull.f[i]);
 160
 161     /* distribute the force over all atoms in the group */
 162     for (j=0;j<pull->pull.ngx[i];j++)
 163       for (m=0;m<DIM;m++) {
 164         mi = md->massT[pull->pull.idx[i][j]];
 165         f[pull->pull.idx[i][j]][m] +=
 166           mi*(pull->pull.f[i][m])/(pull->pull.tmass[i]);
 167       }
 168   }
 169
 170   /* reset center of mass of the reference group to its starting value */
 171   rvec_sub(pull->ref.x_unc[0],pull->ref.x_con[0],cm);
 172   for (i=0;i<pull->ref.ngx[0];i++)
 173     rvec_sub(x[pull->ref.idx[0][i]],cm,x[pull->ref.idx[0][i]]);
 174 }
 175
 176 /* this implements a constraint run like SHAKE does. */
 177 static void do_constraint(t_pull *pull, rvec *x, matrix box, t_mdatoms *md,
 178                           real dt)
 179 {
 180
 181   rvec r_ij, /* x_con[i] com of i in prev. step. Obeys constr. -> r_ij   */
 182     unc_ij,  /* x_unc[i] com of i this step, before constr.    -> unc_ij */
 183     ref_ij;  /* x_ref[i] com of i at t0, not updated           -> ref_ij */
 184
 185   rvec *rinew;           /* current 'new' position of group i */
 186   rvec *rjnew;           /* current 'new' position of group j */
 187   real *direction;       /* direction of dr relative to r_ij */
 188   double lambda, rm, mass, tol = 0.00001;
 189   bool bConverged = FALSE;
 190   int n=0,i,ii,j,m,max_iter=1000;
 191   int ref;
 192   double x1,x2,q,a,b,c;  /* for solving the quadratic equation,
 193                             see Num. Recipes in C ed 2 p. 184 */
 194   rvec *dr;              /* correction for group i */
 195   rvec *ref_dr;          /* correction for group j */
 196   rvec tmp,tmp2,tmp3,sum;
 197
 198   if (pull->bCyl) {
 199     snew(ref_dr,pull->pull.n);
 200     snew(rjnew,pull->pull.n);
 201   } else {
 202     snew(ref_dr,1);
 203     snew(rjnew,1);
 204   }
 205   snew(dr,pull->pull.n);
 206   snew(rinew,pull->pull.n);
 207   snew(direction,pull->pull.n);
 208
 209   /* copy the current unconstraint positions for use in iterations. We
 210      iterate until rinew[i] and rjnew[j] obey the constraints. Then
 211      rinew - pull.x_unc[i] is the correction dr to group i */
 212   for (i=0;i<pull->pull.n;i++)
 213     copy_rvec(pull->pull.x_unc[i],rinew[i]);
 214   if (pull->bCyl)
 215     for (i=0;i<pull->pull.n;i++)
 216       copy_rvec(pull->dyna.x_unc[i],rjnew[i]);
 217   else
 218     copy_rvec(pull->ref.x_unc[0],rjnew[0]);
 219
 220   while (!bConverged && n<max_iter) {
 221     /* loop over all constraints */
 222     for (i=0;(i<pull->pull.n);i++) {
 223
 224       if (pull->bVerbose)
 225         fprintf(stderr,"group %d, iteration %d\n",i,n);
 226
 227       if (pull->bCyl) {
 228         rvec_sub(pull->dyna.x_con[i],pull->pull.x_con[i],r_ij);
 229         rvec_sub(rjnew[i],rinew[i],unc_ij);
 230         rvec_sub(pull->dyna.x_ref[i],pull->pull.x_ref[i],ref_ij);
 231       } else {
 232         rvec_sub(pull->ref.x_con[0],pull->pull.x_con[i],r_ij);
 233         rvec_sub(rjnew[0],rinew[i],unc_ij);
 234         rvec_sub(pull->ref.x_ref[0],pull->pull.x_ref[i],ref_ij);
 235       }
 236
 237       /* select components we want */
 238       for (m=0;m<DIM;m++) {
 239         r_ij[m]     *= pull->dims[m];
 240         unc_ij[m]   *= pull->dims[m];
 241         ref_ij[m]   *= pull->dims[m];
 242
 243         /* correct for PBC. Is this correct? */
 244         if (r_ij[m]   < -0.5*box[m][m]) r_ij[m]   += box[m][m];
 245         if (r_ij[m]   >  0.5*box[m][m]) r_ij[m]   -= box[m][m];
 246         if (unc_ij[m] < -0.5*box[m][m]) unc_ij[m] += box[m][m];
 247         if (unc_ij[m] >  0.5*box[m][m]) unc_ij[m] -= box[m][m];
 248         if (ref_ij[m] < -0.5*box[m][m]) ref_ij[m] += box[m][m];
 249         if (ref_ij[m] >  0.5*box[m][m]) ref_ij[m] -= box[m][m];
 250       }
 251
 252
 253       if (pull->bCyl)
 254         rm = 1/pull->pull.tmass[i] + 1/pull->dyna.tmass[i];
 255       else
 256         rm = 1/pull->pull.tmass[i] + 1/pull->ref.tmass[0];
 257
 258       a = iprod(r_ij,r_ij)*dt*dt*dt*dt*rm*rm;
 259       b = iprod(unc_ij,r_ij)*2*dt*dt*rm;
 260       c = iprod(unc_ij,unc_ij) - norm2(ref_ij);
 261
 262       if (b<0)
 263         q = -0.5*(b-sqrt(b*b-4*a*c));
 264       else
 265         q = -0.5*(b+sqrt(b*b-4*a*c));
 266       x1 = q/a; x2 = c/q;
 267       lambda = x1 > 0 ? x1 : x2;
 268
 269       if (pull->bVerbose)
 270         fprintf(stderr,"\nax^2+bx+c=0: a=%e b=%e c=%e\n"
 271                 "x1=%e x2=%e sum:%e,%e, lambda:%e\n",a,b,c,x1,x2,
 272                 a*x1*x1+b*x1+c,a*x2*x2+b*x2+c,lambda);
 273
 274       /* the position corrections dr due to the constraint are: */
 275       if (pull->bCyl) {
 276         svmul(-dt*dt*lambda/pull->pull.tmass[i],r_ij,dr[i]);
 277         svmul(dt*dt*lambda/pull->dyna.tmass[i],r_ij,ref_dr[i]);
 278       } else {
 279         svmul(-dt*dt*lambda/pull->pull.tmass[i],r_ij,dr[i]);
 280         svmul(dt*dt*lambda/pull->ref.tmass[0],r_ij,ref_dr[0]);
 281       }
 282
 283       /* and the direction of the constraint force: */
 284       direction[i] = cos_angle(r_ij,dr[i]);
 285
 286       /* DEBUG */
 287       if (pull->bVerbose) {
 288         fprintf(stderr,"Direction: %f\n",direction[i]);
 289         if (pull->bCyl) {
 290           rvec_sub(rinew[i],rjnew[i],tmp);
 291           rvec_sub(pull->pull.x_ref[i],pull->dyna.x_ref[i],tmp2);
 292         } else {
 293           rvec_sub(pull->pull.x_ref[i],pull->ref.x_ref[0],tmp2);
 294           rvec_sub(rinew[i],rjnew[0],tmp);
 295         }
 296         rvec_sub(dr[i],ref_dr[0],tmp3);
 297         for (m=0;m<DIM;m++) {
 298           tmp[m]  *= pull->dims[m];
 299           tmp2[m] *= pull->dims[m];
 300           tmp3[m] *= pull->dims[m];
 301           if (tmp[m]  < -0.5*box[m][m]) tmp[m]  += box[m][m];
 302           if (tmp[m]  >  0.5*box[m][m]) tmp[m]  -= box[m][m];
 303           if (tmp2[m] < -0.5*box[m][m]) tmp2[m] += box[m][m];
 304           if (tmp2[m] >  0.5*box[m][m]) tmp2[m] -= box[m][m];
 305           if (tmp3[m] < -0.5*box[m][m]) tmp3[m] += box[m][m];
 306           if (tmp3[m] >  0.5*box[m][m]) tmp3[m] -= box[m][m];
 307         }
 308
 309         if (pull->bCyl)
 310           fprintf(stderr,
 311                   "cur. i:%f %f %f j:%f %f %f d: %f\n"
 312                   "ref. i:%f %f %f j:%f %f %f d: %f\n"
 313                   "cor. i:%f %f %f j:%f %f %f d: %f\n",
 314                   rinew[i][0],rinew[i][1],rinew[i][2],
 315                   rjnew[i][0],rjnew[i][1],rjnew[i][2], norm(tmp),
 316                   pull->pull.x_ref[i][0],pull->pull.x_ref[i][1],
 317                   pull->pull.x_ref[i][2],pull->dyna.x_ref[i][0],
 318                   pull->dyna.x_ref[i][1],pull->dyna.x_ref[i][2],
 319                   norm(tmp2),
 320                   dr[i][0],dr[i][1],dr[i][2],
 321                   ref_dr[0][0],ref_dr[0][1],ref_dr[0][2],
 322                   norm(tmp3));
 323         else
 324           fprintf(stderr,
 325                   "cur. i:%f %f %f j:%f %f %f d: %f\n"
 326                   "ref. i:%f %f %f j:%f %f %f d: %f\n"
 327                   "cor. i:%f %f %f j:%f %f %f d: %f\n",
 328                   rinew[i][0],rinew[i][1],rinew[i][2],
 329                   rjnew[0][0],rjnew[0][1],rjnew[0][2], norm(tmp),
 330                   pull->pull.x_ref[i][0],pull->pull.x_ref[i][1],
 331                   pull->pull.x_ref[i][2],pull->ref.x_ref[0][0],
 332                   pull->ref.x_ref[0][1],pull->ref.x_ref[0][2],
 333                   norm(tmp2),
 334                   dr[i][0],dr[i][1],dr[i][2],
 335                   ref_dr[0][0],ref_dr[0][1],ref_dr[0][2],
 336                   norm(tmp3));
 337       } /* END DEBUG */
 338
 339       /* update the positions with dr */
 340       rvec_add(rinew[i],dr[i],rinew[i]);
 341
 342       if (pull->bCyl) {
 343         rvec_add(rjnew[i],ref_dr[i],rjnew[i]);
 344
 345         /* calculate new distance between the two groups */
 346         rvec_sub(rjnew[i],rinew[i],unc_ij);
 347
 348         /* select components and check PBC again */
 349         for (m=0;m<DIM;m++) {
 350           unc_ij[m] *= pull->dims[m];
 351           if (unc_ij[m] < -0.5*box[m][m]) unc_ij[m] += box[m][m];
 352           if (unc_ij[m] >  0.5*box[m][m]) unc_ij[m] -= box[m][m];
 353         }
 354       } else {
 355         rvec_add(rjnew[0],ref_dr[0],rjnew[0]);
 356
 357         /* calculate new distance between the two groups */
 358         rvec_sub(rjnew[0],rinew[i],unc_ij);
 359
 360         /* select components again and check PBC again */
 361         for (m=0;m<DIM;m++) {
 362           unc_ij[m] *= pull->dims[m];
 363           if (unc_ij[m] < -0.5*box[m][m]) unc_ij[m] += box[m][m];
 364           if (unc_ij[m] >  0.5*box[m][m]) unc_ij[m] -= box[m][m];
 365         }
 366       }
 367     }
 368
 369     /* check if all constraints are fullfilled now */
 370     bConverged = TRUE;
 371     for (i=0;i<pull->pull.n;i++) {
 372       if (pull->bCyl) {
 373         rvec_sub(rjnew[i],rinew[i],unc_ij);
 374         rvec_sub(pull->dyna.x_ref[i],pull->pull.x_ref[i],ref_ij);
 375       } else {
 376         rvec_sub(rjnew[0],rinew[i],unc_ij);
 377         rvec_sub(pull->ref.x_ref[0],pull->pull.x_ref[i],ref_ij);
 378       }
 379
 380       for (m=0;m<DIM;m++) {
 381         ref_ij[m] *= pull->dims[m];
 382         unc_ij[m] *= pull->dims[m];
 383         if (unc_ij[m] < -0.5*box[m][m]) unc_ij[m] += box[m][m];
 384         if (unc_ij[m] >  0.5*box[m][m]) unc_ij[m] -= box[m][m];
 385         if (ref_ij[m] < -0.5*box[m][m]) ref_ij[m] += box[m][m];
 386         if (ref_ij[m] >  0.5*box[m][m]) ref_ij[m] -= box[m][m];
 387       }
 388
 389       bConverged = bConverged && (fabs(norm(unc_ij)-norm(ref_ij)) < tol);
 390     }
 391
 392
 393     /* DEBUG */
 394     if (pull->bVerbose) {
 395       if (!bConverged)
 396         fprintf(stderr,"NOT CONVERGED YET: Group %d (%s):"
 397                 "d_ref = %f, current d = %f\n",
 398                 i,pull->pull.grps[i], norm(ref_ij),norm(unc_ij));
 399     } /* END DEBUG */
 400
 401
 402     n++;
 403     /* if after all constraints are dealt with and bConverged is still TRUE
 404        we're finished, if not we do another iteration */
 405   }
 406   if (n>max_iter)
 407     fatal_error(0,"Too many iterations for constraint run");
 408
 409   /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
 410
 411   /* update the normal groups */
 412   for (i=0;i<pull->pull.n;i++) {
 413     /* get the final dr and constraint force for group i */
 414     rvec_sub(rinew[i],pull->pull.x_unc[i],dr[i]);
 415     /* select components of dr */
 416     for (m=0;m<DIM;m++) dr[i][m] *= pull->dims[m];
 417     svmul(pull->pull.tmass[i]/(dt*dt),dr[i],tmp);
 418     /* get the direction of dr */
 419     pull->pull.f[i][ZZ] = -norm(tmp)*direction[i];
 420
 421     /* copy the new x_unc to x_con */
 422     copy_rvec(rinew[i],pull->pull.x_con[i]);
 423
 424     /* update the atom positions */
 425     clear_rvec(sum);
 426     for (j=0;j<pull->pull.ngx[i];j++) {
 427       ii = pull->pull.idx[i][j];
 428       rvec_add(x[ii],dr[i],x[ii]);
 429       svmul(md->massT[ii],dr[i],tmp);
 430       rvec_add(tmp,sum,sum);
 431     }
 432     if (pull->bVerbose)
 433       fprintf(stderr,"Group %i: correction %e %e %e\n",
 434               i,sum[0],sum[1],sum[2]);
 435   }
 436
 437   /* update the reference groups */
 438   if (pull->bCyl) {
 439     /* update the dynamic reference groups */
 440     for (i=0;i<pull->pull.n;i++) {
 441       rvec_sub(rjnew[i],pull->dyna.x_unc[i],ref_dr[i]);
 442       /* copy the new x_unc to x_con */
 443       copy_rvec(rjnew[i],pull->dyna.x_con[i]);
 444       /* select components of ref_dr */
 445       for (m=0;m<DIM;m++) ref_dr[i][m] *= pull->dims[m];
 446
 447       clear_rvec(sum);
 448       for (j=0;j<pull->dyna.ngx[i];j++) {
 449         /* reset the atoms with dr, weighted by w_i */
 450         svmul(pull->dyna.weights[i][j],ref_dr[i],tmp);
 451         ii = pull->dyna.idx[i][j];
 452         rvec_add(x[ii],tmp,x[ii]);
 453         svmul(md->massT[ii],tmp,tmp2);
 454         rvec_add(tmp2,sum,sum);
 455       }
 456       if (pull->bVerbose)
 457         fprintf(stderr,"Dyna grp %i: correction %e %e %e\n",
 458                 i,sum[0],sum[1],sum[2]);
 459     }
 460   } else {
 461     /* update the reference group */
 462     rvec_sub(rjnew[0],pull->ref.x_unc[0], ref_dr[0]);
 463     /* copy the new x_unc to x_con */
 464     copy_rvec(rjnew[0],pull->ref.x_con[0]);
 465     /* select components of ref_dr */
 466     for (m=0;m<DIM;m++) ref_dr[0][m] *= pull->dims[m];
 467
 468     clear_rvec(sum);
 469     for (j=0;j<pull->ref.ngx[0];j++) {
 470       ii = pull->ref.idx[0][j];
 471       rvec_add(x[ii],ref_dr[0],x[ii]);
 472       svmul(md->massT[ii],ref_dr[0],tmp);
 473       rvec_add(tmp,sum,sum);
 474     }
 475     if (pull->bVerbose)
 476       fprintf(stderr,"Reference: correction %e %e %e\n",
 477               sum[0],sum[1],sum[2]);
 478
 479   }
 480
 481   /* finished! I hope. Give back some memory */
 482   sfree(ref_dr);
 483   sfree(rinew);
 484   sfree(rjnew);
 485   sfree(dr);
 486   sfree(direction);
 487 }
 488
 489 /* mimicks an AFM experiment, groups are pulled via a spring */
 490 static void do_afm(t_pull *pull,rvec *f,matrix box,t_mdatoms *md)
 491 {
 492   int i,ii,j,m,g;
 493   real mi;
 494   rvec cm;     /* center of mass displacement of reference */
 495   rvec dr;     /* extension of the springs */
 496
 497   /* loop over the groups that are being pulled */
 498   for (i=0;i<pull->pull.n;i++) {
 499     /* compute how far the springs are stretched */
 500     for (m=0;m<DIM;m++) {
 501       dr[m]=pull->dims[m]*(pull->pull.spring[i][m]-pull->pull.x_unc[i][m]);
 502       while (dr[m] >  box[m][m]/2) dr[m]-=box[m][m];
 503       while (dr[m] < -box[m][m]/2) dr[m]+=box[m][m];
 504     }
 505
 506     /* calculate force from the spring on the pull group: f = - k*dr */
 507     for (m=0;m<DIM;m++)
 508       pull->pull.f[i][m] = pull->k*dr[m];
 509
 510     /* distribute force on com over individual atoms in the group */
 511     for (j=0;j<pull->pull.ngx[i];j++) {
 512       ii = pull->pull.idx[i][j];
 513       for (m=0;m<DIM;m++) {
 514         mi = md->massT[ii];
 515         f[ii][m] +=  mi*(pull->pull.f[i][m])/(pull->pull.tmass[i]);
 516       }
 517     }
 518
 519     /* move pulling spring along dir, over pull->rate  */
 520     for (m=0;m<DIM;m++)
 521       pull->pull.spring[i][m]+=pull->pull.dir[i][m]*pull->rate;
 522   }
 523   /* done */
 524 }
 525
 526 void pull(t_pull *pull,rvec *x,rvec *f,matrix box, t_topology *top,
 527           real dt, int step, int natoms, t_mdatoms *md)
 528 {
 529   int i;
 530   static rvec *x_s = NULL;
 531   bool bShakeFirst;
 532
 533   bShakeFirst = (f == NULL);
 534
 535   if (!x_s)
 536     snew(x_s,md->nr); /* can't rely on natoms */
 537
 538   /* copy x to temp array x_s. We assume all molecules are whole already */
 539   for (i=0;i<md->nr;i++)
 540     copy_rvec(x[i],x_s[i]);
 541
 542   switch (pull->runtype) {
 543   case eAfm:
 544     if (!bShakeFirst) {
 545       /* calculate center of mass of the pull groups */
 546       for (i=0;i<pull->pull.n;i++)
 547         (void)calc_com(x_s,pull->pull.ngx[i],pull->pull.idx[i],md,
 548                        pull->pull.x_unc[i],box);
 549       do_afm(pull,f,box,md);
 550       print_afm(pull,step);
 551     }
 552     break;
 553
 554   case eStart:
 555     if (!bShakeFirst) {
 556       for (i=0;i<pull->pull.n;i++)
 557         (void)calc_com(x_s,pull->pull.ngx[i],pull->pull.idx[i],md,
 558                        pull->pull.x_unc[i],box);
 559       if (pull->bCyl)
 560         make_refgrps(pull,x_s,md);
 561       else
 562         (void)calc_com(x_s,pull->ref.ngx[0],pull->ref.idx[0],md,
 563                        pull->ref.x_unc[0],box);
 564       do_start(pull,x,box,md,dt,step,top);
 565       print_start(pull,step);
 566     }
 567     break;
 568
 569   case eConstraint:
 570     /* if necessary, correct for particles jumping across the box
 571        this makes sure pull->ref.x0 has the pbc-corrected coordinates
 572        Else just copy the normal coordinates to ref.x0
 573      */
 574     if (pull->reftype == eComT0 || pull->reftype == eDynT0) {
 575       if (pull->bVerbose)
 576         fprintf(stderr,"\nCalling correct_t0_pbc\n");
 577       correct_t0_pbc(pull,x_s,md,box);
 578     } else {
 579       for (i=0;i<pull->ref.ngx[0];i++)
 580         copy_rvec(x_s[pull->ref.idx[0][i]],pull->ref.x0[0][i]);
 581     }
 582
 583     /* get centers of mass for the pull groups. Does this work correctly
 584        with pbc? */
 585     for (i=0;i<pull->pull.n;i++) {
 586       (void)calc_com(x_s,pull->pull.ngx[i],pull->pull.idx[i],md,
 587                      pull->pull.x_unc[i],box);
 588     }
 589
 590     /* get new centers of mass for reference groups, from the, possibly
 591        corrected, pull->ref.x0 */
 592     /* dynamic case */
 593     if (pull->bCyl) {
 594       if (step % pull->update == 0)    /* make new ones? */
 595         make_refgrps(pull,x_s,md);
 596       else {
 597         for (i=0;i<pull->pull.n;i++) {
 598           (void)calc_com2(pull->ref.x0[0],pull->dyna.ngx[i],pull->dyna.idx[i],
 599                           md,pull->dyna.x_unc[i],box);
 600           if (pull->bVerbose)
 601             fprintf(stderr,"dynacom: %8.3f%8.3f%8.3f\n",pull->dyna.x_unc[i][0],
 602                     pull->dyna.x_unc[i][1],pull->dyna.x_unc[i][2]);
 603         }
 604       }
 605     }
 606
 607     /* normal case */
 608     if (!pull->bCyl)
 609       (void)calc_com2(pull->ref.x0[0],pull->ref.ngx[0],pull->ref.idx[0],
 610                       md,pull->ref.x_unc[0],box);
 611
 612     /* if necessary, do a running average over the last reflag steps for com */
 613     if (pull->reflag > 1) {
 614       if (pull->bVerbose)
 615         fprintf(stderr,"Calling calc_running_com\n");
 616       calc_running_com(pull);
 617     }
 618
 619     /* print some debug info if necessary */
 620     if (pull->bVerbose) {
 621       if (pull->bCyl) {
 622         for (i=0;i<pull->pull.n;i++) {
 623           fprintf(stderr,"I      :%9.6f %9.6f %9.6f\n",pull->pull.x_unc[i][0],
 624                   pull->pull.x_unc[i][1],pull->pull.x_unc[i][2]);
 625           fprintf(stderr,"dyna xref: unconstr. com:%9.6f %9.6f %9.6f\n",
 626                   pull->dyna.x_unc[i][0],
 627                   pull->dyna.x_unc[i][1],pull->dyna.x_unc[i][2]);
 628         }
 629       } else {
 630         fprintf(stderr,"xref: unconstr. com:%9.6f %9.6f %9.6f\n",
 631                 pull->ref.x_unc[0][0], pull->ref.x_unc[0][1],
 632                 pull->ref.x_unc[0][2]);
 633       }
 634     }
 635
 636     /* do the actual constraint calculation */
 637     do_constraint(pull,x,box,md,dt);
 638     print_constraint(pull,f,step,box);
 639     break;
 640
 641   case eUmbrella:
 642     if (!bShakeFirst) {
 643       do_umbrella(pull,x,f,box,md);
 644       print_umbrella(pull,step);
 645     }
 646     break;
 647
 648   case eTest:
 649     if (!bShakeFirst) {
 650       /* code to test reference groups, without actually doing anything
 651          else
 652          */
 653       (void)calc_com(x,pull->ref.ngx[0],pull->ref.idx[0],
 654                      md,pull->ref.x_unc[0],box);
 655       fprintf(stderr,"ref: %8.3f %8.3f %8.3f\n",pull->ref.x_unc[0][XX],
 656               pull->ref.x_unc[0][YY],pull->ref.x_unc[0][ZZ]);
 657       correct_t0_pbc(pull,x,md,box);
 658       (void)calc_com2(pull->ref.x0[0],pull->ref.ngx[0],pull->ref.idx[0],
 659                       md,pull->ref.x_unc[0],box);
 660       fprintf(stderr,"ref_t0: %8.3f %8.3f %8.3f\n",pull->ref.x_unc[0][XX],
 661               pull->ref.x_unc[0][YY],pull->ref.x_unc[0][ZZ]);
 662     }
 663     break;
 664
 665   default:
 666     fatal_error(0,"undetermined runtype");
 667   }
 668 }
 669
 670
 671
 672
 673
 674
 675
 676
 677