src/mdlib/domdec_box.c

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