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 "rdgroup.h"
  52 #include "symtab.h"
  53 #include "index.h"
  54 #include "confio.h"
  55 #include "network.h"
  56 #include "pbc.h"
  57 #include "pull.h"
  58 #include "gmx_ga2la.h"
  59
  60 void pull_d_pbc_dx(int npbcdim,matrix box,
  61                    const dvec x1, const dvec x2, dvec dx)
  62 {
  63   int m,d;
  64
  65   /* Only correct for atom pairs with a distance within
  66    * half of the smallest diagonal element of box.
  67    */
  68   dvec_sub(x1,x2,dx);
  69   for(m=npbcdim-1; m>=0; m--) {
  70     while (dx[m] < -0.5*box[m][m]) {
  71       for(d=0; d<DIM; d++)
  72         dx[d] += box[m][d];
  73     }
  74     while (dx[m] >=  0.5*box[m][m]) {
  75       for(d=0; d<DIM; d++)
  76         dx[d] -= box[m][d];
  77     }
  78   }
  79 }
  80
  81 static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg,
  82                              t_mdatoms *md, rvec *x,
  83                              rvec x_pbc)
  84 {
  85   int a,m;
  86
  87   if (cr && PAR(cr)) {
  88     if (DOMAINDECOMP(cr)) {
  89       if (!ga2la_get_home(cr->dd->ga2la,pg->pbcatom,&a)) {
  90         a = -1;
  91       }
  92     } else {
  93       a = pg->pbcatom;
  94     }
  95
  96     if (a >= md->start && a < md->start+md->homenr) {
  97       copy_rvec(x[a],x_pbc);
  98     } else {
  99       clear_rvec(x_pbc);
 100     }
 101   } else {
 102     copy_rvec(x[pg->pbcatom],x_pbc);
 103   }
 104 }
 105
 106 static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
 107                               t_mdatoms *md, rvec *x,
 108                               rvec *x_pbc)
 109 {
 110   int g,n,m;
 111
 112   n = 0;
 113   for(g=0; g<1+pull->ngrp; g++) {
 114     if ((g==0 && PULL_CYL(pull)) || pull->grp[g].pbcatom == -1) {
 115       clear_rvec(x_pbc[g]);
 116     } else {
 117       pull_set_pbcatom(cr,&pull->grp[g],md,x,x_pbc[g]);
 118       for(m=0; m<DIM; m++) {
 119         if (pull->dim[m] == 0) {
 120           x_pbc[g][m] = 0.0;
 121         }
 122       }
 123       n++;
 124     }
 125   }
 126
 127   if (cr && PAR(cr) && n > 0) {
 128     /* Sum over the nodes to get x_pbc from the home node of pbcatom */
 129     gmx_sum((1+pull->ngrp)*DIM,x_pbc[0],cr);
 130   }
 131 }
 132
 133 /* switch function, x between r and w */
 134 static real get_weight(real x, real r1, real r0)
 135 {
 136   real weight;
 137
 138   if (x >= r0)
 139     weight = 0;
 140   else if (x <= r1)
 141     weight = 1;
 142   else
 143     weight = (r0 - x)/(r0 - r1);
 144
 145   return weight;
 146 }
 147
 148 static void make_cyl_refgrps(t_commrec *cr,t_pull *pull,t_mdatoms *md,
 149                              t_pbc *pbc,double t,rvec *x,rvec *xp)
 150 {
 151   int g,i,ii,m,start,end;
 152   rvec g_x,dx,dir;
 153   double r0_2,sum_a,sum_ap,dr2,mass,weight,wmass,wwmass,inp;
 154   t_pullgrp *pref,*pgrp,*pdyna;
 155   gmx_ga2la_t ga2la=NULL;
 156
 157   if (pull->dbuf_cyl == NULL) {
 158     snew(pull->dbuf_cyl,pull->ngrp*4);
 159   }
 160
 161   if (cr && DOMAINDECOMP(cr))
 162     ga2la = cr->dd->ga2la;
 163
 164   start = md->start;
 165   end   = md->homenr + start;
 166
 167   r0_2 = dsqr(pull->cyl_r0);
 168
 169   /* loop over all groups to make a reference group for each*/
 170   pref = &pull->grp[0];
 171   for(g=1; g<1+pull->ngrp; g++) {
 172     pgrp  = &pull->grp[g];
 173     pdyna = &pull->dyna[g];
 174     copy_rvec(pgrp->vec,dir);
 175     sum_a = 0;
 176     sum_ap = 0;
 177     wmass = 0;
 178     wwmass = 0;
 179     pdyna->nat_loc = 0;
 180
 181     for(m=0; m<DIM; m++) {
 182       g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
 183     }
 184
 185     /* loop over all atoms in the main ref group */
 186     for(i=0; i<pref->nat; i++) {
 187       ii = pull->grp[0].ind[i];
 188       if (ga2la) {
 189         if (!ga2la_get_home(ga2la,pref->ind[i],&ii)) {
 190           ii = -1;
 191         }
 192       }
 193       if (ii >= start && ii < end) {
 194         pbc_dx_aiuc(pbc,x[ii],g_x,dx);
 195         inp = iprod(dir,dx);
 196         dr2 = 0;
 197         for(m=0; m<DIM; m++) {
 198           dr2 += dsqr(dx[m] - inp*dir[m]);
 199         }
 200
 201         if (dr2 < r0_2) {
 202           /* add to index, to sum of COM, to weight array */
 203           if (pdyna->nat_loc >= pdyna->nalloc_loc) {
 204             pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
 205             srenew(pdyna->ind_loc,pdyna->nalloc_loc);
 206             srenew(pdyna->weight_loc,pdyna->nalloc_loc);
 207           }
 208           pdyna->ind_loc[pdyna->nat_loc] = ii;
 209           mass = md->massT[ii];
 210           weight = get_weight(sqrt(dr2),pull->cyl_r1,pull->cyl_r0);
 211           pdyna->weight_loc[pdyna->nat_loc] = weight;
 212           sum_a += mass*weight*inp;
 213           if (xp) {
 214             pbc_dx_aiuc(pbc,xp[ii],g_x,dx);
 215             inp = iprod(dir,dx);
 216             sum_ap += mass*weight*inp;
 217           }
 218           wmass += mass*weight;
 219           wwmass += mass*sqr(weight);
 220           pdyna->nat_loc++;
 221         }
 222       }
 223     }
 224     pull->dbuf_cyl[(g-1)*4+0] = wmass;
 225     pull->dbuf_cyl[(g-1)*4+1] = wwmass;
 226     pull->dbuf_cyl[(g-1)*4+2] = sum_a;
 227     pull->dbuf_cyl[(g-1)*4+3] = sum_ap;
 228   }
 229
 230   if (cr && PAR(cr)) {
 231     /* Sum the contributions over the nodes */
 232     gmx_sumd(pull->ngrp*4,pull->dbuf_cyl,cr);
 233   }
 234
 235   for(g=1; g<1+pull->ngrp; g++) {
 236     pgrp  = &pull->grp[g];
 237     pdyna = &pull->dyna[g];
 238
 239     wmass        = pull->dbuf_cyl[(g-1)*4+0];
 240     wwmass       = pull->dbuf_cyl[(g-1)*4+1];
 241     pdyna->wscale = wmass/wwmass;
 242     pdyna->invtm = 1.0/(pdyna->wscale*wmass);
 243
 244     for(m=0; m<DIM; m++) {
 245       g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
 246       pdyna->x[m]    = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+2]/wmass;
 247       if (xp) {
 248         pdyna->xp[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+3]/wmass;
 249       }
 250     }
 251
 252     if (debug) {
 253       fprintf(debug,"Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
 254               g,pdyna->x[0],pdyna->x[1],
 255               pdyna->x[2],1.0/pdyna->invtm);
 256     }
 257   }
 258 }
 259
 260 static double atan2_0_2pi(double y,double x)
 261 {
 262   double a;
 263
 264   a = atan2(y,x);
 265   if (a < 0) {
 266     a += 2.0*M_PI;
 267   }
 268   return a;
 269 }
 270
 271 /* calculates center of mass of selection index from all coordinates x */
 272 void pull_calc_coms(t_commrec *cr,
 273                     t_pull *pull, t_mdatoms *md, t_pbc *pbc,double t,
 274                     rvec x[], rvec *xp)
 275 {
 276   int  g,i,ii,m;
 277   real mass,w,wm,twopi_box=0;
 278   double wmass,wwmass,invwmass;
 279   dvec com,comp;
 280   double cm,sm,cmp,smp,ccm,csm,ssm,csw,snw;
 281   rvec *xx[2],x_pbc={0,0,0},dx;
 282   t_pullgrp *pgrp;
 283
 284   if (pull->rbuf == NULL) {
 285     snew(pull->rbuf,1+pull->ngrp);
 286   }
 287   if (pull->dbuf == NULL) {
 288     snew(pull->dbuf,3*(1+pull->ngrp));
 289   }
 290
 291   if (pull->bRefAt) {
 292     pull_set_pbcatoms(cr,pull,md,x,pull->rbuf);
 293   }
 294
 295   if (pull->cosdim >= 0) {
 296     for(m=pull->cosdim+1; m<pull->npbcdim; m++) {
 297       if (pbc->box[m][pull->cosdim] != 0) {
 298         gmx_fatal(FARGS,"Can not do cosine weighting for trilinic dimensions");
 299       }
 300     }
 301     twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
 302   }
 303
 304   for (g=0; g<1+pull->ngrp; g++) {
 305     pgrp = &pull->grp[g];
 306     clear_dvec(com);
 307     clear_dvec(comp);
 308     wmass  = 0;
 309     wwmass = 0;
 310     cm  = 0;
 311     sm  = 0;
 312     cmp = 0;
 313     smp = 0;
 314     ccm = 0;
 315     csm = 0;
 316     ssm = 0;
 317     if (!(g==0 && PULL_CYL(pull))) {
 318       if (pgrp->epgrppbc == epgrppbcREFAT) {
 319         /* Set the pbc atom */
 320         copy_rvec(pull->rbuf[g],x_pbc);
 321       }
 322       w = 1;
 323       for(i=0; i<pgrp->nat_loc; i++) {
 324         ii = pgrp->ind_loc[i];
 325         mass = md->massT[ii];
 326         if (pgrp->epgrppbc != epgrppbcCOS) {
 327           if (pgrp->weight_loc) {
 328             w = pgrp->weight_loc[i];
 329           }
 330           wm = w*mass;
 331           wmass  += wm;
 332           wwmass += wm*w;
 333           if (pgrp->epgrppbc == epgrppbcNONE) {
 334             /* Plain COM: sum the coordinates */
 335             for(m=0; m<DIM; m++)
 336               com[m]    += wm*x[ii][m];
 337             if (xp) {
 338               for(m=0; m<DIM; m++)
 339                 comp[m] += wm*xp[ii][m];
 340             }
 341           } else {
 342             /* Sum the difference with the reference atom */
 343             pbc_dx(pbc,x[ii],x_pbc,dx);
 344             for(m=0; m<DIM; m++) {
 345               com[m]    += wm*dx[m];
 346             }
 347             if (xp) {
 348               /* For xp add the difference between xp and x to dx,
 349                * such that we use the same periodic image,
 350                * also when xp has a large displacement.
 351                */
 352               for(m=0; m<DIM; m++) {
 353                 comp[m] += wm*(dx[m] + xp[ii][m] - x[ii][m]);
 354               }
 355             }
 356           }
 357         } else {
 358           /* Determine cos and sin sums */
 359           csw = cos(x[ii][pull->cosdim]*twopi_box);
 360           snw = sin(x[ii][pull->cosdim]*twopi_box);
 361           cm  += csw*mass;
 362           sm  += snw*mass;
 363           ccm += csw*csw*mass;
 364           csm += csw*snw*mass;
 365           ssm += snw*snw*mass;
 366
 367           if (xp) {
 368             csw = cos(xp[ii][pull->cosdim]*twopi_box);
 369             snw = sin(xp[ii][pull->cosdim]*twopi_box);
 370             cmp += csw*mass;
 371             smp += snw*mass;
 372           }
 373         }
 374       }
 375     }
 376
 377     /* Copy local sums to a buffer for global summing */
 378     switch (pgrp->epgrppbc) {
 379     case epgrppbcNONE:
 380     case epgrppbcREFAT:
 381       copy_dvec(com,pull->dbuf[g*3]);
 382       copy_dvec(comp,pull->dbuf[g*3+1]);
 383       pull->dbuf[g*3+2][0] = wmass;
 384       pull->dbuf[g*3+2][1] = wwmass;
 385       pull->dbuf[g*3+2][2] = 0;
 386       break;
 387     case epgrppbcCOS:
 388       pull->dbuf[g*3  ][0] = cm;
 389       pull->dbuf[g*3  ][1] = sm;
 390       pull->dbuf[g*3  ][2] = 0;
 391       pull->dbuf[g*3+1][0] = ccm;
 392       pull->dbuf[g*3+1][1] = csm;
 393       pull->dbuf[g*3+1][2] = ssm;
 394       pull->dbuf[g*3+2][0] = cmp;
 395       pull->dbuf[g*3+2][1] = smp;
 396       pull->dbuf[g*3+2][2] = 0;
 397       break;
 398     }
 399   }
 400
 401   if (cr && PAR(cr)) {
 402     /* Sum the contributions over the nodes */
 403     gmx_sumd((1+pull->ngrp)*3*DIM,pull->dbuf[0],cr);
 404   }
 405
 406   for (g=0; g<1+pull->ngrp; g++) {
 407     pgrp = &pull->grp[g];
 408     if (pgrp->nat > 0 && !(g==0 && PULL_CYL(pull))) {
 409       if (pgrp->epgrppbc != epgrppbcCOS) {
 410         /* Determine the inverse mass */
 411         wmass  = pull->dbuf[g*3+2][0];
 412         wwmass = pull->dbuf[g*3+2][1];
 413         invwmass = 1/wmass;
 414         /* invtm==0 signals a frozen group, so then we should keep it zero */
 415         if (pgrp->invtm > 0) {
 416           pgrp->wscale = wmass/wwmass;
 417           pgrp->invtm  = 1.0/(pgrp->wscale*wmass);
 418         }
 419         /* Divide by the total mass */
 420         for(m=0; m<DIM; m++) {
 421           pgrp->x[m]    = pull->dbuf[g*3  ][m]*invwmass;
 422           if (xp) {
 423             pgrp->xp[m] = pull->dbuf[g*3+1][m]*invwmass;
 424           }
 425           if (pgrp->epgrppbc == epgrppbcREFAT) {
 426             pgrp->x[m]    += pull->rbuf[g][m];
 427             if (xp) {
 428               pgrp->xp[m] += pull->rbuf[g][m];
 429             }
 430           }
 431         }
 432       } else {
 433         /* Determine the optimal location of the cosine weight */
 434         csw = pull->dbuf[g*3][0];
 435         snw = pull->dbuf[g*3][1];
 436         pgrp->x[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
 437         /* Set the weights for the local atoms */
 438         wmass = sqrt(csw*csw + snw*snw);
 439         wwmass = (pull->dbuf[g*3+1][0]*csw*csw +
 440                   pull->dbuf[g*3+1][1]*csw*snw +
 441                   pull->dbuf[g*3+1][2]*snw*snw)/(wmass*wmass);
 442         pgrp->wscale = wmass/wwmass;
 443         pgrp->invtm  = 1.0/(pgrp->wscale*wmass);
 444         /* Set the weights for the local atoms */
 445         csw *= pgrp->invtm;
 446         snw *= pgrp->invtm;
 447         for(i=0; i<pgrp->nat_loc; i++) {
 448           ii = pgrp->ind_loc[i];
 449           pgrp->weight_loc[i] = csw*cos(twopi_box*x[ii][pull->cosdim]) +
 450                                 snw*sin(twopi_box*x[ii][pull->cosdim]);
 451         }
 452         if (xp) {
 453           csw = pull->dbuf[g*3+2][0];
 454           snw = pull->dbuf[g*3+2][1];
 455           pgrp->xp[pull->cosdim] = atan2_0_2pi(snw,csw)/twopi_box;
 456         }
 457       }
 458       if (debug) {
 459         fprintf(debug,"Pull group %d wmass %f wwmass %f invtm %f\n",
 460                 g,wmass,wwmass,pgrp->invtm);
 461       }
 462     }
 463   }
 464
 465   if (PULL_CYL(pull)) {
 466     /* Calculate the COMs for the cyclinder reference groups */
 467     make_cyl_refgrps(cr,pull,md,pbc,t,x,xp);
 468   }
 469 }