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