src/mdlib/pullutil.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40 #include <stdlib.h>
  41 #include "sysstuff.h"
  42 #include "princ.h"
  43 #include "futil.h"
  44 #include "statutil.h"
  45 #include "vec.h"
  46 #include "smalloc.h"
  47 #include "typedefs.h"
  48 #include "names.h"
  49 #include "gmx_fatal.h"
  50 #include "macros.h"
  51 #include "symtab.h"
  52 #include "index.h"
  53 #include "confio.h"
  54 #include "network.h"
  55 #include "pbc.h"
  56 #include "pull.h"
  57 #include "gmx_ga2la.h"
  58
  59 static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg,
  60                              t_mdatoms *md, rvec *x,
  61                              rvec x_pbc)
  62 {
  63   int a,m;
  64
  65   if (cr && PAR(cr)) {
  66     if (DOMAINDECOMP(cr)) {
  67       if (!ga2la_get_home(cr->dd->ga2la,pg->pbcatom,&a)) {
  68         a = -1;
  69       }
  70     } else {
  71       a = pg->pbcatom;
  72     }
  73
  74     if (a >= md->start && a < md->start+md->homenr) {
  75       copy_rvec(x[a],x_pbc);
  76     } else {
  77       clear_rvec(x_pbc);
  78     }
  79   } else {
  80     copy_rvec(x[pg->pbcatom],x_pbc);
  81   }
  82 }
  83
  84 static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
  85                               t_mdatoms *md, rvec *x,
  86                               rvec *x_pbc)
  87 {
  88   int g,n,m;
  89
  90   n = 0;
  91   for(g=0; g<1+pull->ngrp; g++) {
  92     if ((g==0 && PULL_CYL(pull)) || pull->grp[g].pbcatom == -1) {
  93       clear_rvec(x_pbc[g]);
  94     } else {
  95       pull_set_pbcatom(cr,&pull->grp[g],md,x,x_pbc[g]);
  96       for(m=0; m<DIM; m++) {
  97         if (pull->dim[m] == 0) {
  98           x_pbc[g][m] = 0.0;
  99         }
 100       }
 101       n++;
 102     }
 103   }
 104
 105   if (cr && PAR(cr) && n > 0) {
 106     /* Sum over the nodes to get x_pbc from the home node of pbcatom */
 107     gmx_sum((1+pull->ngrp)*DIM,x_pbc[0],cr);
 108   }
 109 }
 110
 111 /* switch function, x between r and w */
 112 static real get_weight(real x, real r1, real r0)
 113 {
 114   real weight;
 115
 116   if (x >= r0)
 117     weight = 0;
 118   else if (x <= r1)
 119     weight = 1;
 120   else
 121     weight = (r0 - x)/(r0 - r1);
 122
 123   return weight;
 124 }
 125
 126 static void make_cyl_refgrps(t_commrec *cr,t_pull *pull,t_mdatoms *md,
 127                              t_pbc *pbc,double t,rvec *x,rvec *xp)
 128 {
 129   int g,i,ii,m,start,end;
 130   rvec g_x,dx,dir;
 131   double r0_2,sum_a,sum_ap,dr2,mass,weight,wmass,wwmass,inp;
 132   t_pullgrp *pref,*pgrp,*pdyna;
 133   gmx_ga2la_t ga2la=NULL;
 134
 135   if (pull->dbuf_cyl == NULL) {
 136     snew(pull->dbuf_cyl,pull->ngrp*4);
 137   }
 138
 139   if (cr && DOMAINDECOMP(cr))
 140     ga2la = cr->dd->ga2la;
 141
 142   start = md->start;
 143   end   = md->homenr + start;
 144
 145   r0_2 = dsqr(pull->cyl_r0);
 146
 147   /* loop over all groups to make a reference group for each*/
 148   pref = &pull->grp[0];
 149   for(g=1; g<1+pull->ngrp; g++) {
 150     pgrp  = &pull->grp[g];
 151     pdyna = &pull->dyna[g];
 152     copy_rvec(pgrp->vec,dir);
 153     sum_a = 0;
 154     sum_ap = 0;
 155     wmass = 0;
 156     wwmass = 0;
 157     pdyna->nat_loc = 0;
 158
 159     for(m=0; m<DIM; m++) {
 160       g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
 161     }
 162
 163     /* loop over all atoms in the main ref group */
 164     for(i=0; i<pref->nat; i++) {
 165       ii = pull->grp[0].ind[i];
 166       if (ga2la) {
 167         if (!ga2la_get_home(ga2la,pref->ind[i],&ii)) {
 168           ii = -1;
 169         }
 170       }
 171       if (ii >= start && ii < end) {
 172         pbc_dx_aiuc(pbc,x[ii],g_x,dx);
 173         inp = iprod(dir,dx);
 174         dr2 = 0;
 175         for(m=0; m<DIM; m++) {
 176           dr2 += dsqr(dx[m] - inp*dir[m]);
 177         }
 178
 179         if (dr2 < r0_2) {
 180           /* add to index, to sum of COM, to weight array */
 181           if (pdyna->nat_loc >= pdyna->nalloc_loc) {
 182             pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
 183             srenew(pdyna->ind_loc,pdyna->nalloc_loc);
 184             srenew(pdyna->weight_loc,pdyna->nalloc_loc);
 185           }
 186           pdyna->ind_loc[pdyna->nat_loc] = ii;
 187           mass = md->massT[ii];
 188           weight = get_weight(sqrt(dr2),pull->cyl_r1,pull->cyl_r0);
 189           pdyna->weight_loc[pdyna->nat_loc] = weight;
 190           sum_a += mass*weight*inp;
 191           if (xp) {
 192             pbc_dx_aiuc(pbc,xp[ii],g_x,dx);
 193             inp = iprod(dir,dx);
 194             sum_ap += mass*weight*inp;
 195           }
 196           wmass += mass*weight;
 197           wwmass += mass*sqr(weight);
 198           pdyna->nat_loc++;
 199         }
 200       }
 201     }
 202     pull->dbuf_cyl[(g-1)*4+0] = wmass;
 203     pull->dbuf_cyl[(g-1)*4+1] = wwmass;
 204     pull->dbuf_cyl[(g-1)*4+2] = sum_a;
 205     pull->dbuf_cyl[(g-1)*4+3] = sum_ap;
 206   }
 207
 208   if (cr && PAR(cr)) {
 209     /* Sum the contributions over the nodes */
 210     gmx_sumd(pull->ngrp*4,pull->dbuf_cyl,cr);
 211   }
 212
 213   for(g=1; g<1+pull->ngrp; g++) {
 214     pgrp  = &pull->grp[g];
 215     pdyna = &pull->dyna[g];
 216
 217     wmass        = pull->dbuf_cyl[(g-1)*4+0];
 218     wwmass       = pull->dbuf_cyl[(g-1)*4+1];
 219     pdyna->wscale = wmass/wwmass;
 220     pdyna->invtm = 1.0/(pdyna->wscale*wmass);
 221
 222     for(m=0; m<DIM; m++) {
 223       g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
 224       pdyna->x[m]    = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+2]/wmass;
 225       if (xp) {
 226         pdyna->xp[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+3]/wmass;
 227       }
 228     }
 229
 230     if (debug) {
 231       fprintf(debug,"Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
 232               g,pdyna->x[0],pdyna->x[1],
 233               pdyna->x[2],1.0/pdyna->invtm);
 234     }
 235   }
 236 }
 237
 238 static double atan2_0_2pi(double y,double x)
 239 {
 240   double a;
 241
 242   a = atan2(y,x);
 243   if (a < 0) {
 244     a += 2.0*M_PI;
 245   }
 246   return a;
 247 }
 248
 249 /* calculates center of mass of selection index from all coordinates x */
 250 void pull_calc_coms(t_commrec *cr,
 251                     t_pull *pull, t_mdatoms *md, t_pbc *pbc,double t,
 252                     rvec x[], rvec *xp)
 253 {
 254   int  g,i,ii,m;
 255   real mass,w,wm,twopi_box=0;
 256   double wmass,wwmass,invwmass;
 257   dvec com,comp;
 258   double cm,sm,cmp,smp,ccm,csm,ssm,csw,snw;
 259   rvec *xx[2],x_pbc={0,0,0},dx;
 260   t_pullgrp *pgrp;
 261
 262   if (pull->rbuf == NULL) {
 263     snew(pull->rbuf,1+pull->ngrp);
 264   }
 265   if (pull->dbuf == NULL) {
 266     snew(pull->dbuf,3*(1+pull->ngrp));
 267   }
 268
 269   if (pull->bRefAt) {
 270     pull_set_pbcatoms(cr,pull,md,x,pull->rbuf);
 271   }
 272
 273   if (pull->cosdim >= 0) {
 274     for(m=pull->cosdim+1; m<pull->npbcdim; m++) {
 275       if (pbc->box[m][pull->cosdim] != 0) {
 276         gmx_fatal(FARGS,"Can not do cosine weighting for trilinic dimensions");
 277       }
 278     }
 279     twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
 280   }
 281
 282   for (g=0; g<1+pull->ngrp; g++) {
 283     pgrp = &pull->grp[g];
 284     clear_dvec(com);
 285     clear_dvec(comp);
 286     wmass  = 0;
 287     wwmass = 0;
 288     cm  = 0;
 289     sm  = 0;
 290     cmp = 0;
 291     smp = 0;
 292     ccm = 0;
 293     csm = 0;
 294     ssm = 0;
 295     if (!(g==0 && PULL_CYL(pull))) {
 296       if (pgrp->epgrppbc == epgrppbcREFAT) {
 297         /* Set the pbc atom */
 298         copy_rvec(pull->rbuf[g],x_pbc);
 299       }
 300       w = 1;
 301       for(i=0; i<pgrp->nat_loc; i++) {
 302         ii = pgrp->ind_loc[i];
 303         mass = md->massT[ii];
 304         if (pgrp->epgrppbc != epgrppbcCOS) {
 305           if (pgrp->weight_loc) {
 306             w = pgrp->weight_loc[i];
 307           }
 308           wm = w*mass;
 309           wmass  += wm;
 310           wwmass += wm*w;
 311           if (pgrp->epgrppbc == epgrppbcNONE) {
 312             /* Plain COM: sum the coordinates */
 313             for(m=0; m<DIM; m++)
 314               com[m]    += wm*x[ii][m];
 315             if (xp) {
 316               for(m=0; m<DIM; m++)
 317                 comp[m] += wm*xp[ii][m];
 318             }
 319           } else {
 320             /* Sum the difference with the reference atom */
 321             pbc_dx(pbc,x[ii],x_pbc,dx);
 322             for(m=0; m<DIM; m++) {
 323               com[m]    += wm*dx[m];
 324             }
 325             if (xp) {
 326               /* For xp add the difference between xp and x to dx,
 327                * such that we use the same periodic image,
 328                * also when xp has a large displacement.
 329                */
 330               for(m=0; m<DIM; m++) {
 331                 comp[m] += wm*(dx[m] + xp[ii][m] - x[ii][m]);
 332               }
 333             }
 334           }
 335         } else {
 336           /* Determine cos and sin sums */
 337           csw = cos(x[ii][pull->cosdim]*twopi_box);
 338           snw = sin(x[ii][pull->cosdim]*twopi_box);
 339           cm  += csw*mass;
 340           sm  += snw*mass;
 341           ccm += csw*csw*mass;
 342           csm += csw*snw*mass;
 343           ssm += snw*snw*mass;
 344
 345           if (xp) {
 346             csw = cos(xp[ii][pull->cosdim]*twopi_box);
 347             snw = sin(xp[ii][pull->cosdim]*twopi_box);
 348             cmp += csw*mass;
 349             smp += snw*mass;
 350           }
 351         }
 352       }
 353     }
 354
 355     /* Copy local sums to a buffer for global summing */
 356     switch (pgrp->epgrppbc) {
 357     case epgrppbcNONE:
 358     case epgrppbcREFAT:
 359       copy_dvec(com,pull->dbuf[g*3]);
 360       copy_dvec(comp,pull->dbuf[g*3+1]);
 361       pull->dbuf[g*3+2][0] = wmass;
 362       pull->dbuf[g*3+2][1] = wwmass;
 363       pull->dbuf[g*3+2][2] = 0;
 364       break;
 365     case epgrppbcCOS:
 366       pull->dbuf[g*3  ][0] = cm;
 367       pull->dbuf[g*3  ][1] = sm;
 368       pull->dbuf[g*3  ][2] = 0;
 369       pull->dbuf[g*3+1][0] = ccm;
 370       pull->dbuf[g*3+1][1] = csm;
 371       pull->dbuf[g*3+1][2] = ssm;
 372       pull->dbuf[g*3+2][0] = cmp;
 373       pull->dbuf[g*3+2][1] = smp;
 374       pull->dbuf[g*3+2][2] = 0;
 375       break;
 376     }
 377   }
 378
 379   if (cr && PAR(cr)) {
 380     /* Sum the contributions over the nodes */
 381     gmx_sumd((1+pull->ngrp)*3*DIM,pull->dbuf[0],cr);
 382   }
 383
 384   for (g=0; g<1+pull->ngrp; g++) {
 385     pgrp = &pull->grp[g];
 386     if (pgrp->nat > 0 && !(g==0 && PULL_CYL(pull))) {
 387       if (pgrp->epgrppbc != epgrppbcCOS) {
 388         /* Determine the inverse mass */
 389         wmass  = pull->dbuf[g*3+2][0];
 390         wwmass = pull->dbuf[g*3+2][1];
 391         invwmass = 1/wmass;
 392         /* invtm==0 signals a frozen group, so then we should keep it zero */
 393         if (pgrp->invtm > 0) {
 394           pgrp->wscale = wmass/wwmass;
 395           pgrp->invtm  = 1.0/(pgrp->wscale*wmass);
 396         }
 397         /* Divide by the total mass */
 398         for(m=0; m<DIM; m++) {
 399           pgrp->x[m]    = pull->dbuf[g*3  ][m]*invwmass;
 400           if (xp) {
 401             pgrp->xp[m] = pull->dbuf[g*3+1][m]*invwmass;
 402           }
 403           if (pgrp->epgrppbc == epgrppbcREFAT) {
 404             pgrp->x[m]    += pull->rbuf[g][m];
 405             if (xp) {
 406               pgrp->xp[m] += pull->rbuf[g][m];
 407             }
 408           }
 409         }
 410       } else {
 411         /* Determine the optimal location of the cosine weight */
 412         csw = pull->dbuf[g*3][0];
 413         snw = pull->dbuf[g*3][1];
 414         pgrp->x[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
 415         /* Set the weights for the local atoms */
 416         wmass = sqrt(csw*csw + snw*snw);
 417         wwmass = (pull->dbuf[g*3+1][0]*csw*csw +
 418                   pull->dbuf[g*3+1][1]*csw*snw +
 419                   pull->dbuf[g*3+1][2]*snw*snw)/(wmass*wmass);
 420         pgrp->wscale = wmass/wwmass;
 421         pgrp->invtm  = 1.0/(pgrp->wscale*wmass);
 422         /* Set the weights for the local atoms */
 423         csw *= pgrp->invtm;
 424         snw *= pgrp->invtm;
 425         for(i=0; i<pgrp->nat_loc; i++) {
 426           ii = pgrp->ind_loc[i];
 427           pgrp->weight_loc[i] = csw*cos(twopi_box*x[ii][pull->cosdim]) +
 428                                 snw*sin(twopi_box*x[ii][pull->cosdim]);
 429         }
 430         if (xp) {
 431           csw = pull->dbuf[g*3+2][0];
 432           snw = pull->dbuf[g*3+2][1];
 433           pgrp->xp[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
 434         }
 435       }
 436       if (debug) {
 437         fprintf(debug,"Pull group %d wmass %f wwmass %f invtm %f\n",
 438                 g,wmass,wwmass,pgrp->invtm);
 439       }
 440     }
 441   }
 442
 443   if (PULL_CYL(pull)) {
 444     /* Calculate the COMs for the cyclinder reference groups */
 445     make_cyl_refgrps(cr,pull,md,pbc,t,x,xp);
 446   }
 447 }