src/mdlib/pullutil.c

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