src/mdlib/groupcoord.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  * $Id: gmx_matrix.c,v 1.4 2008/12/02 18:27:57 spoel Exp $
   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 4.5
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2008, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Groningen Machine for Chemical Simulation
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40
  41 #include "groupcoord.h"
  42 #include "mpelogging.h"
  43 #include "network.h"
  44 #include "pbc.h"
  45 #include "vec.h"
  46 #include "smalloc.h"
  47 #include "gmx_ga2la.h"
  48
  49 #define MIN(a,b) (((a)<(b))?(a):(b))
  50
  51
  52
  53 /* Select the indices of the group's atoms which are local and store them in
  54  * anrs_loc[0..nr_loc]. The indices are saved in coll_ind[] for later reduction
  55  * in communicate_group_positions()
  56  */
  57 extern void dd_make_local_group_indices(
  58         gmx_ga2la_t   ga2la,
  59         const int     nr,          /* IN:  Total number of atoms in the group */
  60         int           anrs[],      /* IN:  Global atom numbers of the groups atoms */
  61         int           *nr_loc,     /* OUT: Number of group atoms found locally */
  62         int           *anrs_loc[], /* OUT: Local atom numbers of the group  */
  63         int           *nalloc_loc, /* IN+OUT: Allocation size of anrs_loc */
  64         int           coll_ind[])  /* OUT (opt): Where is this position found in the collective array? */
  65 {
  66     int  i,ii;
  67     int  localnr;
  68
  69
  70     /* Loop over all the atom indices of the group to check
  71      * which ones are on the local node */
  72     localnr = 0;
  73     for(i=0; i<nr; i++)
  74     {
  75         if (ga2la_get_home(ga2la,anrs[i],&ii))
  76         {
  77             /* The atom with this index is a home atom */
  78             if (localnr >= *nalloc_loc) /* Check whether memory suffices */
  79             {
  80                 *nalloc_loc = over_alloc_dd(localnr+1);
  81                 /* We never need more memory than the number of atoms in the group */
  82                 *nalloc_loc = MIN(*nalloc_loc, nr);
  83                 srenew(*anrs_loc,*nalloc_loc);
  84             }
  85             /* Save the atoms index in the local atom numbers array */
  86             (*anrs_loc)[localnr] = ii;
  87
  88             if (coll_ind != NULL)
  89             {
  90                 /* Keep track of where this local atom belongs in the collective index array.
  91                  * This is needed when reducing the local arrays to a collective/global array
  92                  * in communicate_group_positions */
  93                 coll_ind[localnr] = i;
  94             }
  95
  96             /* add one to the local atom count */
  97             localnr++;
  98         }
  99     }
 100
 101     /* Return the number of local atoms that were found */
 102     *nr_loc = localnr;
 103 }
 104
 105
 106 static void get_shifts_group(
 107         int    npbcdim,
 108         matrix box,
 109         rvec   *xcoll,     /* IN:  Collective set of positions [0..nr] */
 110         int    nr,         /* IN:  Total number of atoms in the group */
 111         rvec   *xcoll_old, /* IN:  Positions from the last time step [0...nr] */
 112         ivec   *shifts)    /* OUT: Shifts for xcoll */
 113 {
 114     int  i,m,d;
 115     rvec dx;
 116
 117
 118     /* Get the shifts such that each atom is within closest
 119      * distance to its position at the last NS time step after shifting.
 120      * If we start with a whole group, and always keep track of
 121      * shift changes, the group will stay whole this way */
 122     for (i=0; i < nr; i++)
 123         clear_ivec(shifts[i]);
 124
 125     for (i=0; i<nr; i++)
 126     {
 127         /* The distance this atom moved since the last time step */
 128         /* If this is more than just a bit, it has changed its home pbc box */
 129         rvec_sub(xcoll[i],xcoll_old[i],dx);
 130
 131         for(m=npbcdim-1; m>=0; m--)
 132         {
 133             while (dx[m] < -0.5*box[m][m])
 134             {
 135                 for(d=0; d<DIM; d++)
 136                     dx[d] += box[m][d];
 137                 shifts[i][m]++;
 138             }
 139             while (dx[m] >= 0.5*box[m][m])
 140             {
 141                 for(d=0; d<DIM; d++)
 142                     dx[d] -= box[m][d];
 143                 shifts[i][m]--;
 144             }
 145         }
 146     }
 147 }
 148
 149
 150 static void shift_positions_group(
 151         matrix box,
 152         rvec   x[],      /* The positions [0..nr] */
 153         ivec   *is,      /* The shifts [0..nr] */
 154         int    nr)       /* The number of positions and shifts */
 155 {
 156     int      i,tx,ty,tz;
 157
 158
 159     GMX_MPE_LOG(ev_shift_start);
 160
 161     /* Loop over the group's atoms */
 162     if(TRICLINIC(box))
 163     {
 164         for (i=0; i < nr; i++)
 165         {
 166             tx=is[i][XX];
 167             ty=is[i][YY];
 168             tz=is[i][ZZ];
 169
 170             x[i][XX]=x[i][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
 171             x[i][YY]=x[i][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
 172             x[i][ZZ]=x[i][ZZ]+tz*box[ZZ][ZZ];
 173         }
 174     } else
 175     {
 176         for (i=0; i < nr; i++)
 177         {
 178             tx=is[i][XX];
 179             ty=is[i][YY];
 180             tz=is[i][ZZ];
 181
 182             x[i][XX]=x[i][XX]+tx*box[XX][XX];
 183             x[i][YY]=x[i][YY]+ty*box[YY][YY];
 184             x[i][ZZ]=x[i][ZZ]+tz*box[ZZ][ZZ];
 185         }
 186     }
 187     GMX_MPE_LOG(ev_shift_finish);
 188 }
 189
 190
 191 /* Assemble the positions of the group such that every node has all of them.
 192  * The atom indices are retrieved from anrs_loc[0..nr_loc]
 193  * Note that coll_ind[i] = i is needed in the serial case */
 194 extern void communicate_group_positions(
 195         t_commrec  *cr,
 196         rvec       *xcoll,        /* OUT: Collective array of positions */
 197         ivec       *shifts,       /* IN+OUT: Collective array of shifts for xcoll */
 198         ivec       *extra_shifts, /* BUF: Extra shifts since last time step */
 199         const gmx_bool bNS,           /* IN:  NS step, the shifts have changed */
 200         rvec       *x_loc,        /* IN:  Local positions on this node */
 201         const int  nr,            /* IN:  Total number of atoms in the group */
 202         const int  nr_loc,        /* IN:  Local number of atoms in the group */
 203         int        *anrs_loc,     /* IN:  Local atom numbers */
 204         int        *coll_ind,     /* IN:  Collective index */
 205         rvec       *xcoll_old,    /* IN+OUT: Positions from the last time step, used to make group whole */
 206         matrix     box)
 207 {
 208     int i;
 209
 210
 211     GMX_MPE_LOG(ev_get_group_x_start);
 212
 213     /* Zero out the groups' global position array */
 214     clear_rvecs(nr, xcoll);
 215
 216     /* Put the local positions that this node has into the right place of
 217      * the collective array. Note that in the serial case, coll_ind[i] = i */
 218     for (i=0; i<nr_loc; i++)
 219         copy_rvec(x_loc[anrs_loc[i]], xcoll[coll_ind[i]]);
 220
 221     if (PAR(cr))
 222     {
 223         /* Add the arrays from all nodes together */
 224         gmx_sum(nr*3, xcoll[0], cr);
 225
 226         /* To make the group whole, start with a whole group and each
 227          * step move the assembled positions at closest distance to the positions
 228          * from the last step. First shift the positions with the saved shift
 229          * vectors (these are 0 when this routine is called for the first time!) */
 230         shift_positions_group(box, xcoll, shifts, nr);
 231
 232         /* Now check if some shifts changed since the last step.
 233          * This only needs to be done when the shifts are expected to have changed,
 234          * i.e. after neighboursearching */
 235         if (bNS)
 236         {
 237             get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
 238
 239             /* Shift with the additional shifts such that we get a whole group now */
 240             shift_positions_group(box, xcoll, extra_shifts, nr);
 241
 242             /* Add the shift vectors together for the next time step */
 243             for (i=0; i<nr; i++)
 244             {
 245                 shifts[i][XX] += extra_shifts[i][XX];
 246                 shifts[i][YY] += extra_shifts[i][YY];
 247                 shifts[i][ZZ] += extra_shifts[i][ZZ];
 248             }
 249
 250             /* Store current correctly-shifted positions for comparison in the next NS time step */
 251             for (i=0; i<nr; i++)
 252                 copy_rvec(xcoll[i],xcoll_old[i]);
 253         }
 254     }
 255
 256     GMX_MPE_LOG(ev_get_group_x_finish);
 257 }
 258
 259
 260 /* Determine the (weighted) sum vector from positions x */
 261 extern double get_sum_of_positions(rvec x[], real weight[], const int nat, dvec dsumvec)
 262 {
 263     int i;
 264     rvec x_weighted;
 265     double weight_sum = 0.0;
 266
 267
 268     /* Zero out the center */
 269     clear_dvec(dsumvec);
 270
 271     /* Loop over all atoms and add their weighted position vectors */
 272     if (weight != NULL)
 273     {
 274         for (i=0; i<nat; i++)
 275         {
 276             weight_sum += weight[i];
 277             svmul(weight[i], x[i], x_weighted);
 278             dsumvec[XX] += x_weighted[XX];
 279             dsumvec[YY] += x_weighted[YY];
 280             dsumvec[ZZ] += x_weighted[ZZ];
 281         }
 282     }
 283     else
 284     {
 285         for (i=0; i<nat; i++)
 286         {
 287             dsumvec[XX] += x[i][XX];
 288             dsumvec[YY] += x[i][YY];
 289             dsumvec[ZZ] += x[i][ZZ];
 290         }
 291     }
 292     return weight_sum;
 293 }
 294
 295
 296 /* Determine center of structure from collective positions x */
 297 extern void get_center(rvec x[], real weight[], const int nr, rvec rcenter)
 298 {
 299     dvec   dcenter;
 300     double weight_sum, denom;
 301
 302
 303     weight_sum = get_sum_of_positions(x, weight, nr, dcenter);
 304
 305     if (weight != NULL)
 306         denom = weight_sum; /* Divide by the sum of weight */
 307     else
 308         denom = nr;        /* Divide by the number of atoms */
 309
 310     dsvmul(1.0/denom, dcenter, dcenter);
 311
 312     rcenter[XX] = dcenter[XX];
 313     rcenter[YY] = dcenter[YY];
 314     rcenter[ZZ] = dcenter[ZZ];
 315 }
 316
 317
 318 /* Get the center from local positions that already have the correct
 319  * PBC representation */
 320 extern void get_center_comm(
 321         t_commrec *cr,
 322         rvec x_loc[],       /* Local positions */
 323         real weight_loc[],  /* Local masses or other weights */
 324         int nr_loc,         /* Local number of atoms */
 325         int nr_group,       /* Total number of atoms of the group */
 326         rvec center)        /* Weighted center */
 327 {
 328     double weight_sum, denom;
 329     dvec   dsumvec;
 330     double buf[4];
 331
 332
 333     weight_sum = get_sum_of_positions(x_loc, weight_loc, nr_loc, dsumvec);
 334
 335     /* Add the local contributions from all nodes. Put the sum vector and the
 336      * weight in a buffer array so that we get along with a single communication
 337      * call. */
 338     if (PAR(cr))
 339     {
 340         buf[0] = dsumvec[XX];
 341         buf[1] = dsumvec[YY];
 342         buf[2] = dsumvec[ZZ];
 343         buf[3] = weight_sum;
 344
 345         /* Communicate buffer */
 346         gmx_sumd(4, buf, cr);
 347
 348         dsumvec[XX] = buf[0];
 349         dsumvec[YY] = buf[1];
 350         dsumvec[ZZ] = buf[2];
 351         weight_sum  = buf[3];
 352     }
 353
 354     if (weight_loc != NULL)
 355         denom = 1.0/weight_sum; /* Divide by the sum of weight to get center of mass e.g. */
 356     else
 357         denom = 1.0/nr_group;   /* Divide by the number of atoms to get the geometrical center */
 358
 359     center[XX] = dsumvec[XX]*denom;
 360     center[YY] = dsumvec[YY]*denom;
 361     center[ZZ] = dsumvec[ZZ]*denom;
 362 }
 363
 364
 365 /* Translate x with transvec */
 366 extern void translate_x(rvec x[], const int nr, const rvec transvec)
 367 {
 368     int i;
 369
 370
 371     for (i=0; i<nr; i++)
 372         rvec_inc(x[i], transvec);
 373 }
 374
 375
 376 extern void rotate_x(rvec x[], const int nr, matrix rmat)
 377 {
 378     int i,j,k;
 379     rvec x_old;
 380
 381
 382     /* Apply the rotation matrix */
 383     for (i=0; i<nr; i++)
 384     {
 385         for (j=0; j<3; j++)
 386             x_old[j] = x[i][j];
 387         for (j=0; j<3; j++)
 388         {
 389             x[i][j] = 0;
 390             for (k=0; k<3; k++)
 391                 x[i][j] += rmat[k][j]*x_old[k];
 392         }
 393     }
 394 }
 395