src/mdlib/domdec_box.c

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