src/gromacs/mdlib/domdec_box.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 #include "gmxpre.h"
  37
  38 #include "gromacs/legacyheaders/domdec.h"
  39 #include "gromacs/legacyheaders/domdec_network.h"
  40 #include "gromacs/legacyheaders/network.h"
  41 #include "gromacs/legacyheaders/nsgrid.h"
  42 #include "gromacs/legacyheaders/typedefs.h"
  43 #include "gromacs/legacyheaders/types/commrec.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/pbcutil/pbc.h"
  46 #include "gromacs/utility/fatalerror.h"
  47
  48 static void calc_cgcm_av_stddev(t_block *cgs, int n, rvec *x, rvec av, rvec stddev,
  49                                 t_commrec *cr_sum)
  50 {
  51     int   *cgindex;
  52     dvec   s1, s2;
  53     double buf[7];
  54     int    cg, d, k0, k1, k, nrcg;
  55     real   inv_ncg;
  56     rvec   cg_cm;
  57
  58     clear_dvec(s1);
  59     clear_dvec(s2);
  60
  61     cgindex = cgs->index;
  62     for (cg = 0; cg < n; cg++)
  63     {
  64         k0      = cgindex[cg];
  65         k1      = cgindex[cg+1];
  66         nrcg    = k1 - k0;
  67         if (nrcg == 1)
  68         {
  69             copy_rvec(x[k0], cg_cm);
  70         }
  71         else
  72         {
  73             inv_ncg = 1.0/nrcg;
  74
  75             clear_rvec(cg_cm);
  76             for (k = k0; (k < k1); k++)
  77             {
  78                 rvec_inc(cg_cm, x[k]);
  79             }
  80             for (d = 0; (d < DIM); d++)
  81             {
  82                 cg_cm[d] *= inv_ncg;
  83             }
  84         }
  85         for (d = 0; d < DIM; d++)
  86         {
  87             s1[d] += cg_cm[d];
  88             s2[d] += cg_cm[d]*cg_cm[d];
  89         }
  90     }
  91
  92     if (cr_sum != NULL)
  93     {
  94         for (d = 0; d < DIM; d++)
  95         {
  96             buf[d]     = s1[d];
  97             buf[DIM+d] = s2[d];
  98         }
  99         buf[6] = n;
 100         gmx_sumd(7, buf, cr_sum);
 101         for (d = 0; d < DIM; d++)
 102         {
 103             s1[d] = buf[d];
 104             s2[d] = buf[DIM+d];
 105         }
 106         n = (int)(buf[6] + 0.5);
 107     }
 108
 109     dsvmul(1.0/n, s1, s1);
 110     dsvmul(1.0/n, s2, s2);
 111
 112     for (d = 0; d < DIM; d++)
 113     {
 114         av[d]     = s1[d];
 115         stddev[d] = sqrt(s2[d] - s1[d]*s1[d]);
 116     }
 117 }
 118
 119 static void set_tric_dir(ivec *dd_nc, gmx_ddbox_t *ddbox, matrix box)
 120 {
 121     int   npbcdim, d, i, j;
 122     rvec *v, *normal;
 123     real  dep, inv_skew_fac2;
 124
 125     npbcdim = ddbox->npbcdim;
 126     normal  = ddbox->normal;
 127     for (d = 0; d < DIM; d++)
 128     {
 129         ddbox->tric_dir[d] = 0;
 130         for (j = d+1; j < npbcdim; j++)
 131         {
 132             if (box[j][d] != 0)
 133             {
 134                 ddbox->tric_dir[d] = 1;
 135                 if (dd_nc != NULL && (*dd_nc)[j] > 1 && (*dd_nc)[d] == 1)
 136                 {
 137                     gmx_fatal(FARGS, "Domain decomposition has not been implemented for box vectors that have non-zero components in directions that do not use domain decomposition: ncells = %d %d %d, box vector[%d] = %f %f %f",
 138                               dd_nc[XX], dd_nc[YY], dd_nc[ZZ],
 139                               j+1, box[j][XX], box[j][YY], box[j][ZZ]);
 140                 }
 141             }
 142         }
 143
 144         /* Convert box vectors to orthogonal vectors for this dimension,
 145          * for use in distance calculations.
 146          * Set the trilinic skewing factor that translates
 147          * the thickness of a slab perpendicular to this dimension
 148          * into the real thickness of the slab.
 149          */
 150         if (ddbox->tric_dir[d])
 151         {
 152             inv_skew_fac2 = 1;
 153             v             = ddbox->v[d];
 154             if (d == XX || d == YY)
 155             {
 156                 /* Normalize such that the "diagonal" is 1 */
 157                 svmul(1/box[d+1][d+1], box[d+1], v[d+1]);
 158                 for (i = 0; i < d; i++)
 159                 {
 160                     v[d+1][i] = 0;
 161                 }
 162                 inv_skew_fac2 += sqr(v[d+1][d]);
 163                 if (d == XX)
 164                 {
 165                     /* Normalize such that the "diagonal" is 1 */
 166                     svmul(1/box[d+2][d+2], box[d+2], v[d+2]);
 167                     for (i = 0; i < d; i++)
 168                     {
 169                         v[d+2][i] = 0;
 170                     }
 171                     /* Make vector [d+2] perpendicular to vector [d+1],
 172                      * this does not affect the normalization.
 173                      */
 174                     dep = iprod(v[d+1], v[d+2])/norm2(v[d+1]);
 175                     for (i = 0; i < DIM; i++)
 176                     {
 177                         v[d+2][i] -= dep*v[d+1][i];
 178                     }
 179                     inv_skew_fac2 += sqr(v[d+2][d]);
 180
 181                     cprod(v[d+1], v[d+2], normal[d]);
 182                 }
 183                 else
 184                 {
 185                     /* cross product with (1,0,0) */
 186                     normal[d][XX] =  0;
 187                     normal[d][YY] =  v[d+1][ZZ];
 188                     normal[d][ZZ] = -v[d+1][YY];
 189                 }
 190                 if (debug)
 191                 {
 192                     fprintf(debug, "box[%d]  %.3f %.3f %.3f\n",
 193                             d, box[d][XX], box[d][YY], box[d][ZZ]);
 194                     for (i = d+1; i < DIM; i++)
 195                     {
 196                         fprintf(debug, "  v[%d]  %.3f %.3f %.3f\n",
 197                                 i, v[i][XX], v[i][YY], v[i][ZZ]);
 198                     }
 199                 }
 200             }
 201             ddbox->skew_fac[d] = 1.0/sqrt(inv_skew_fac2);
 202             /* Set the normal vector length to skew_fac */
 203             dep = ddbox->skew_fac[d]/norm(normal[d]);
 204             svmul(dep, normal[d], normal[d]);
 205
 206             if (debug)
 207             {
 208                 fprintf(debug, "skew_fac[%d] = %f\n", d, ddbox->skew_fac[d]);
 209                 fprintf(debug, "normal[%d]  %.3f %.3f %.3f\n",
 210                         d, normal[d][XX], normal[d][YY], normal[d][ZZ]);
 211             }
 212         }
 213         else
 214         {
 215             ddbox->skew_fac[d] = 1;
 216
 217             for (i = 0; i < DIM; i++)
 218             {
 219                 clear_rvec(ddbox->v[d][i]);
 220                 ddbox->v[d][i][i] = 1;
 221             }
 222             clear_rvec(normal[d]);
 223             normal[d][d] = 1;
 224         }
 225     }
 226 }
 227
 228 static void low_set_ddbox(t_inputrec *ir, ivec *dd_nc, matrix box,
 229                           gmx_bool bCalcUnboundedSize, int ncg, t_block *cgs, rvec *x,
 230                           t_commrec *cr_sum,
 231                           gmx_ddbox_t *ddbox)
 232 {
 233     rvec av, stddev;
 234     real b0, b1;
 235     int  d;
 236
 237     ddbox->npbcdim     = ePBC2npbcdim(ir->ePBC);
 238     ddbox->nboundeddim = inputrec2nboundeddim(ir);
 239
 240     for (d = 0; d < ddbox->nboundeddim; d++)
 241     {
 242         ddbox->box0[d]     = 0;
 243         ddbox->box_size[d] = box[d][d];
 244     }
 245
 246     if (ddbox->nboundeddim < DIM && bCalcUnboundedSize)
 247     {
 248         calc_cgcm_av_stddev(cgs, ncg, x, av, stddev, cr_sum);
 249
 250         /* GRID_STDDEV_FAC * stddev
 251          * gives a uniform load for a rectangular block of cg's.
 252          * For a sphere it is not a bad approximation for 4x1x1 up to 4x2x2.
 253          */
 254         for (d = ddbox->nboundeddim; d < DIM; d++)
 255         {
 256             b0 = av[d] - GRID_STDDEV_FAC*stddev[d];
 257             b1 = av[d] + GRID_STDDEV_FAC*stddev[d];
 258             if (debug)
 259             {
 260                 fprintf(debug, "Setting global DD grid boundaries to %f - %f\n",
 261                         b0, b1);
 262             }
 263             ddbox->box0[d]     = b0;
 264             ddbox->box_size[d] = b1 - b0;
 265         }
 266     }
 267
 268     set_tric_dir(dd_nc, ddbox, box);
 269 }
 270
 271 void set_ddbox(gmx_domdec_t *dd, gmx_bool bMasterState, t_commrec *cr_sum,
 272                t_inputrec *ir, matrix box,
 273                gmx_bool bCalcUnboundedSize, t_block *cgs, rvec *x,
 274                gmx_ddbox_t *ddbox)
 275 {
 276     if (!bMasterState || DDMASTER(dd))
 277     {
 278         low_set_ddbox(ir, &dd->nc, box, bCalcUnboundedSize,
 279                       bMasterState ? cgs->nr : dd->ncg_home, cgs, x,
 280                       bMasterState ? NULL : cr_sum,
 281                       ddbox);
 282     }
 283
 284     if (bMasterState)
 285     {
 286         dd_bcast(dd, sizeof(gmx_ddbox_t), ddbox);
 287     }
 288 }
 289
 290 void set_ddbox_cr(t_commrec *cr, ivec *dd_nc,
 291                   t_inputrec *ir, matrix box, t_block *cgs, rvec *x,
 292                   gmx_ddbox_t *ddbox)
 293 {
 294     if (MASTER(cr))
 295     {
 296         low_set_ddbox(ir, dd_nc, box, TRUE, cgs->nr, cgs, x, NULL, ddbox);
 297     }
 298
 299     gmx_bcast(sizeof(gmx_ddbox_t), ddbox, cr);
 300 }