src/gromacs/trajectoryanalysis/modules/surfacearea.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2007, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "surfacearea.h"
  40
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  44
  45 #include <cmath>
  46
  47 #include <algorithm>
  48 #include <vector>
  49
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/utilities.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/selection/nbsearch.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/gmxassert.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 using namespace gmx;
  60
  61 #define UNSP_ICO_DOD      9
  62 #define UNSP_ICO_ARC     10
  63
  64 #define FOURPI (4.*M_PI)
  65 #define TORAD(A)     ((A)*0.017453293)
  66 #define DP_TOL     0.001
  67
  68 static real safe_asin(real f)
  69 {
  70     if ( (fabs(f) < 1.00) )
  71     {
  72         return( asin(f) );
  73     }
  74     GMX_ASSERT(fabs(f) - 1.0 > DP_TOL, "Invalid argument");
  75     return(M_PI_2);
  76 }
  77
  78 /* routines for dot distributions on the surface of the unit sphere */
  79 static real icosaeder_vertices(real *xus)
  80 {
  81     const real rh = std::sqrt(1.-2.*cos(TORAD(72.)))/(1.-cos(TORAD(72.)));
  82     const real rg = cos(TORAD(72.))/(1.-cos(TORAD(72.)));
  83     /* icosaeder vertices */
  84     xus[ 0] = 0.;                  xus[ 1] = 0.;                  xus[ 2] = 1.;
  85     xus[ 3] = rh*cos(TORAD(72.));  xus[ 4] = rh*sin(TORAD(72.));  xus[ 5] = rg;
  86     xus[ 6] = rh*cos(TORAD(144.)); xus[ 7] = rh*sin(TORAD(144.)); xus[ 8] = rg;
  87     xus[ 9] = rh*cos(TORAD(216.)); xus[10] = rh*sin(TORAD(216.)); xus[11] = rg;
  88     xus[12] = rh*cos(TORAD(288.)); xus[13] = rh*sin(TORAD(288.)); xus[14] = rg;
  89     xus[15] = rh;                  xus[16] = 0;                   xus[17] = rg;
  90     xus[18] = rh*cos(TORAD(36.));  xus[19] = rh*sin(TORAD(36.));  xus[20] = -rg;
  91     xus[21] = rh*cos(TORAD(108.)); xus[22] = rh*sin(TORAD(108.)); xus[23] = -rg;
  92     xus[24] = -rh;                 xus[25] = 0;                   xus[26] = -rg;
  93     xus[27] = rh*cos(TORAD(252.)); xus[28] = rh*sin(TORAD(252.)); xus[29] = -rg;
  94     xus[30] = rh*cos(TORAD(324.)); xus[31] = rh*sin(TORAD(324.)); xus[32] = -rg;
  95     xus[33] = 0.;                  xus[34] = 0.;                  xus[35] = -1.;
  96     return rh;
  97 }
  98
  99
 100 static void divarc(real x1, real y1, real z1,
 101                    real x2, real y2, real z2,
 102                    int div1, int div2, real *xr, real *yr, real *zr)
 103 {
 104
 105     real xd, yd, zd, dd, d1, d2, s, x, y, z;
 106     real phi, sphi, cphi;
 107
 108     xd = y1*z2-y2*z1;
 109     yd = z1*x2-z2*x1;
 110     zd = x1*y2-x2*y1;
 111     dd = std::sqrt(xd*xd+yd*yd+zd*zd);
 112     GMX_ASSERT(dd >= DP_TOL, "Rotation axis vector too short");
 113
 114     d1 = x1*x1+y1*y1+z1*z1;
 115     d2 = x2*x2+y2*y2+z2*z2;
 116     GMX_ASSERT(d1 >= 0.5, "Vector 1 too short");
 117     GMX_ASSERT(d2 >= 0.5, "Vector 2 too short");
 118
 119     phi  = safe_asin(dd/std::sqrt(d1*d2));
 120     phi  = phi*(static_cast<real>(div1))/(static_cast<real>(div2));
 121     sphi = sin(phi); cphi = cos(phi);
 122     s    = (x1*xd+y1*yd+z1*zd)/dd;
 123
 124     x   = xd*s*(1.-cphi)/dd + x1 * cphi + (yd*z1-y1*zd)*sphi/dd;
 125     y   = yd*s*(1.-cphi)/dd + y1 * cphi + (zd*x1-z1*xd)*sphi/dd;
 126     z   = zd*s*(1.-cphi)/dd + z1 * cphi + (xd*y1-x1*yd)*sphi/dd;
 127     dd  = std::sqrt(x*x+y*y+z*z);
 128     *xr = x/dd; *yr = y/dd; *zr = z/dd;
 129 }
 130
 131 /* densit...required dots per unit sphere */
 132 static std::vector<real> ico_dot_arc(int densit)
 133 {
 134     /* dot distribution on a unit sphere based on an icosaeder *
 135      * great circle average refining of icosahedral face       */
 136
 137     int   i, j, k, tl, tl2, tn;
 138     real  a, d, x, y, z, x2, y2, z2, x3, y3, z3;
 139     real  xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
 140           xjk, yjk, zjk, xkj, ykj, zkj;
 141
 142     /* calculate tessalation level */
 143     a = std::sqrt(((static_cast<real>(densit))-2.)/10.);
 144     const int  tess = static_cast<int>(ceil(a));
 145     const int  ndot = 10*tess*tess+2;
 146     GMX_RELEASE_ASSERT(ndot >= densit, "Inconsistent surface dot formula");
 147
 148     std::vector<real> xus(3*ndot);
 149     const real        rh = icosaeder_vertices(xus.data());
 150
 151     if (tess > 1)
 152     {
 153         tn = 12;
 154         a  = rh*rh*2.*(1.-cos(TORAD(72.)));
 155         /* calculate tessalation of icosaeder edges */
 156         for (i = 0; i < 11; i++)
 157         {
 158             for (j = i+1; j < 12; j++)
 159             {
 160                 x = xus[3*i]-xus[3*j];
 161                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 162                 d = x*x+y*y+z*z;
 163                 if (std::fabs(a-d) > DP_TOL)
 164                 {
 165                     continue;
 166                 }
 167                 for (tl = 1; tl < tess; tl++)
 168                 {
 169                     GMX_ASSERT(tn < ndot, "Inconsistent precomputed surface dot count");
 170                     divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 171                            xus[3*j], xus[1+3*j], xus[2+3*j],
 172                            tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
 173                     tn++;
 174                 }
 175             }
 176         }
 177         /* calculate tessalation of icosaeder faces */
 178         for (i = 0; i < 10; i++)
 179         {
 180             for (j = i+1; j < 11; j++)
 181             {
 182                 x = xus[3*i]-xus[3*j];
 183                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 184                 d = x*x+y*y+z*z;
 185                 if (std::fabs(a-d) > DP_TOL)
 186                 {
 187                     continue;
 188                 }
 189
 190                 for (k = j+1; k < 12; k++)
 191                 {
 192                     x = xus[3*i]-xus[3*k];
 193                     y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 194                     d = x*x+y*y+z*z;
 195                     if (std::fabs(a-d) > DP_TOL)
 196                     {
 197                         continue;
 198                     }
 199                     x = xus[3*j]-xus[3*k];
 200                     y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 201                     d = x*x+y*y+z*z;
 202                     if (std::fabs(a-d) > DP_TOL)
 203                     {
 204                         continue;
 205                     }
 206                     for (tl = 1; tl < tess-1; tl++)
 207                     {
 208                         divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 209                                xus[3*i], xus[1+3*i], xus[2+3*i],
 210                                tl, tess, &xji, &yji, &zji);
 211                         divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 212                                xus[3*i], xus[1+3*i], xus[2+3*i],
 213                                tl, tess, &xki, &yki, &zki);
 214
 215                         for (tl2 = 1; tl2 < tess-tl; tl2++)
 216                         {
 217                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 218                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 219                                    tl2, tess, &xij, &yij, &zij);
 220                             divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 221                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 222                                    tl2, tess, &xkj, &ykj, &zkj);
 223                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 224                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 225                                    tess-tl-tl2, tess, &xik, &yik, &zik);
 226                             divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 227                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 228                                    tess-tl-tl2, tess, &xjk, &yjk, &zjk);
 229                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
 230                                    &x, &y, &z);
 231                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
 232                                    &x2, &y2, &z2);
 233                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
 234                                    &x3, &y3, &z3);
 235                             GMX_ASSERT(tn < ndot,
 236                                        "Inconsistent precomputed surface dot count");
 237                             x           = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
 238                             d           = std::sqrt(x*x+y*y+z*z);
 239                             xus[3*tn]   = x/d;
 240                             xus[1+3*tn] = y/d;
 241                             xus[2+3*tn] = z/d;
 242                             tn++;
 243                         } /* cycle tl2 */
 244                     }     /* cycle tl */
 245                 }         /* cycle k */
 246             }             /* cycle j */
 247         }                 /* cycle i */
 248         GMX_ASSERT(tn == ndot, "Inconsistent precomputed surface dot count");
 249     }                     /* end of if (tess > 1) */
 250
 251     return xus;
 252 }                                  /* end of routine ico_dot_arc */
 253
 254 /* densit...required dots per unit sphere */
 255 static std::vector<real> ico_dot_dod(int densit)
 256 {
 257     /* dot distribution on a unit sphere based on an icosaeder *
 258      * great circle average refining of icosahedral face       */
 259
 260     int   i, j, k, tl, tl2, tn, tess, j1, j2;
 261     real  a, d, x, y, z, x2, y2, z2, x3, y3, z3, ai_d, adod;
 262     real  xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki,
 263           xjk, yjk, zjk, xkj, ykj, zkj;
 264
 265     /* calculate tesselation level */
 266     a     = std::sqrt(((static_cast<real>(densit))-2.)/30.);
 267     tess  = std::max(static_cast<int>(ceil(a)), 1);
 268     const int ndot = 30*tess*tess+2;
 269     GMX_RELEASE_ASSERT(ndot >= densit, "Inconsistent surface dot formula");
 270
 271     std::vector<real> xus(3*ndot);
 272     const real        rh = icosaeder_vertices(xus.data());
 273
 274     tn = 12;
 275     /* square of the edge of an icosaeder */
 276     a = rh*rh*2.*(1.-cos(TORAD(72.)));
 277     /* dodecaeder vertices */
 278     for (i = 0; i < 10; i++)
 279     {
 280         for (j = i+1; j < 11; j++)
 281         {
 282             x = xus[3*i]-xus[3*j];
 283             y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 284             d = x*x+y*y+z*z;
 285             if (std::fabs(a-d) > DP_TOL)
 286             {
 287                 continue;
 288             }
 289             for (k = j+1; k < 12; k++)
 290             {
 291                 x = xus[3*i]-xus[3*k];
 292                 y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 293                 d = x*x+y*y+z*z;
 294                 if (std::fabs(a-d) > DP_TOL)
 295                 {
 296                     continue;
 297                 }
 298                 x = xus[3*j]-xus[3*k];
 299                 y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 300                 d = x*x+y*y+z*z;
 301                 if (std::fabs(a-d) > DP_TOL)
 302                 {
 303                     continue;
 304                 }
 305                 x         = xus[  3*i]+xus[  3*j]+xus[  3*k];
 306                 y         = xus[1+3*i]+xus[1+3*j]+xus[1+3*k];
 307                 z         = xus[2+3*i]+xus[2+3*j]+xus[2+3*k];
 308                 d         = std::sqrt(x*x+y*y+z*z);
 309                 xus[3*tn] = x/d; xus[1+3*tn] = y/d; xus[2+3*tn] = z/d;
 310                 tn++;
 311             }
 312         }
 313     }
 314
 315     if (tess > 1)
 316     {
 317         tn = 32;
 318         /* square of the edge of an dodecaeder */
 319         adod = 4.*(cos(TORAD(108.))-cos(TORAD(120.)))/(1.-cos(TORAD(120.)));
 320         /* square of the distance of two adjacent vertices of ico- and dodecaeder */
 321         ai_d = 2.*(1.-std::sqrt(1.-a/3.));
 322
 323         /* calculate tessalation of mixed edges */
 324         for (i = 0; i < 31; i++)
 325         {
 326             j1 = 12; j2 = 32; a = ai_d;
 327             if (i >= 12)
 328             {
 329                 j1 = i+1; a = adod;
 330             }
 331             for (j = j1; j < j2; j++)
 332             {
 333                 x = xus[3*i]-xus[3*j];
 334                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 335                 d = x*x+y*y+z*z;
 336                 if (std::fabs(a-d) > DP_TOL)
 337                 {
 338                     continue;
 339                 }
 340                 for (tl = 1; tl < tess; tl++)
 341                 {
 342                     GMX_ASSERT(tn < ndot,
 343                                "Inconsistent precomputed surface dot count");
 344                     divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 345                            xus[3*j], xus[1+3*j], xus[2+3*j],
 346                            tl, tess, &xus[3*tn], &xus[1+3*tn], &xus[2+3*tn]);
 347                     tn++;
 348                 }
 349             }
 350         }
 351         /* calculate tessalation of pentakisdodecahedron faces */
 352         for (i = 0; i < 12; i++)
 353         {
 354             for (j = 12; j < 31; j++)
 355             {
 356                 x = xus[3*i]-xus[3*j];
 357                 y = xus[1+3*i]-xus[1+3*j]; z = xus[2+3*i]-xus[2+3*j];
 358                 d = x*x+y*y+z*z;
 359                 if (std::fabs(ai_d-d) > DP_TOL)
 360                 {
 361                     continue;
 362                 }
 363
 364                 for (k = j+1; k < 32; k++)
 365                 {
 366                     x = xus[3*i]-xus[3*k];
 367                     y = xus[1+3*i]-xus[1+3*k]; z = xus[2+3*i]-xus[2+3*k];
 368                     d = x*x+y*y+z*z;
 369                     if (std::fabs(ai_d-d) > DP_TOL)
 370                     {
 371                         continue;
 372                     }
 373                     x = xus[3*j]-xus[3*k];
 374                     y = xus[1+3*j]-xus[1+3*k]; z = xus[2+3*j]-xus[2+3*k];
 375                     d = x*x+y*y+z*z;
 376                     if (std::fabs(adod-d) > DP_TOL)
 377                     {
 378                         continue;
 379                     }
 380                     for (tl = 1; tl < tess-1; tl++)
 381                     {
 382                         divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 383                                xus[3*i], xus[1+3*i], xus[2+3*i],
 384                                tl, tess, &xji, &yji, &zji);
 385                         divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 386                                xus[3*i], xus[1+3*i], xus[2+3*i],
 387                                tl, tess, &xki, &yki, &zki);
 388
 389                         for (tl2 = 1; tl2 < tess-tl; tl2++)
 390                         {
 391                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 392                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 393                                    tl2, tess, &xij, &yij, &zij);
 394                             divarc(xus[3*k], xus[1+3*k], xus[2+3*k],
 395                                    xus[3*j], xus[1+3*j], xus[2+3*j],
 396                                    tl2, tess, &xkj, &ykj, &zkj);
 397                             divarc(xus[3*i], xus[1+3*i], xus[2+3*i],
 398                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 399                                    tess-tl-tl2, tess, &xik, &yik, &zik);
 400                             divarc(xus[3*j], xus[1+3*j], xus[2+3*j],
 401                                    xus[3*k], xus[1+3*k], xus[2+3*k],
 402                                    tess-tl-tl2, tess, &xjk, &yjk, &zjk);
 403                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess-tl,
 404                                    &x, &y, &z);
 405                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess-tl2,
 406                                    &x2, &y2, &z2);
 407                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl+tl2,
 408                                    &x3, &y3, &z3);
 409                             GMX_ASSERT(tn < ndot,
 410                                        "Inconsistent precomputed surface dot count");
 411                             x           = x+x2+x3; y = y+y2+y3; z = z+z2+z3;
 412                             d           = std::sqrt(x*x+y*y+z*z);
 413                             xus[3*tn]   = x/d;
 414                             xus[1+3*tn] = y/d;
 415                             xus[2+3*tn] = z/d;
 416                             tn++;
 417                         } /* cycle tl2 */
 418                     }     /* cycle tl */
 419                 }         /* cycle k */
 420             }             /* cycle j */
 421         }                 /* cycle i */
 422         GMX_ASSERT(tn == ndot, "Inconsistent precomputed surface dot count");
 423     }                     /* end of if (tess > 1) */
 424
 425     return xus;
 426 }       /* end of routine ico_dot_dod */
 427
 428 static int unsp_type(int densit)
 429 {
 430     int i1, i2;
 431     i1 = 1;
 432     while (10*i1*i1+2 < densit)
 433     {
 434         i1++;
 435     }
 436     i2 = 1;
 437     while (30*i2*i2+2 < densit)
 438     {
 439         i2++;
 440     }
 441     if (10*i1*i1-2 < 30*i2*i2-2)
 442     {
 443         return UNSP_ICO_ARC;
 444     }
 445     else
 446     {
 447         return UNSP_ICO_DOD;
 448     }
 449 }
 450
 451 static std::vector<real> make_unsp(int densit, int cubus)
 452 {
 453     int              *ico_wk, *ico_pt;
 454     int               ico_cube, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
 455     int              *work;
 456     real              x, y, z;
 457
 458     int               mode = unsp_type(densit);
 459     std::vector<real> xus;
 460     if (mode == UNSP_ICO_ARC)
 461     {
 462         xus = ico_dot_arc(densit);
 463     }
 464     else if (mode == UNSP_ICO_DOD)
 465     {
 466         xus = ico_dot_dod(densit);
 467     }
 468     else
 469     {
 470         GMX_RELEASE_ASSERT(false, "Invalid unit sphere mode");
 471     }
 472
 473     const int ndot = static_cast<int>(xus.size())/3;
 474
 475     /* determine distribution of points in elementary cubes */
 476     if (cubus)
 477     {
 478         ico_cube = cubus;
 479     }
 480     else
 481     {
 482         i          = 1;
 483         while (i*i*i*2 < ndot)
 484         {
 485             i++;
 486         }
 487         ico_cube = std::max(i-1, 0);
 488     }
 489     ico_cube_cb = ico_cube*ico_cube*ico_cube;
 490     const real del_cube = 2./(static_cast<real>(ico_cube));
 491     snew(work, ndot);
 492     for (l = 0; l < ndot; l++)
 493     {
 494         i = std::max(static_cast<int>(floor((1.+xus[3*l])/del_cube)), 0);
 495         if (i >= ico_cube)
 496         {
 497             i = ico_cube-1;
 498         }
 499         j = std::max(static_cast<int>(floor((1.+xus[1+3*l])/del_cube)), 0);
 500         if (j >= ico_cube)
 501         {
 502             j = ico_cube-1;
 503         }
 504         k = std::max(static_cast<int>(floor((1.+xus[2+3*l])/del_cube)), 0);
 505         if (k >= ico_cube)
 506         {
 507             k = ico_cube-1;
 508         }
 509         ijk     = i+j*ico_cube+k*ico_cube*ico_cube;
 510         work[l] = ijk;
 511     }
 512
 513     snew(ico_wk, 2*ico_cube_cb+1);
 514
 515     ico_pt = ico_wk+ico_cube_cb;
 516     for (l = 0; l < ndot; l++)
 517     {
 518         ico_wk[work[l]]++; /* dots per elementary cube */
 519     }
 520
 521     /* reordering of the coordinate array in accordance with box number */
 522     tn = 0;
 523     for (i = 0; i < ico_cube; i++)
 524     {
 525         for (j = 0; j < ico_cube; j++)
 526         {
 527             for (k = 0; k < ico_cube; k++)
 528             {
 529                 tl            = 0;
 530                 tl2           = tn;
 531                 ijk           = i+ico_cube*j+ico_cube*ico_cube*k;
 532                 *(ico_pt+ijk) = tn;
 533                 for (l = tl2; l < ndot; l++)
 534                 {
 535                     if (ijk == work[l])
 536                     {
 537                         x          = xus[3*l]; y = xus[1+3*l]; z = xus[2+3*l];
 538                         xus[3*l]   = xus[3*tn];
 539                         xus[1+3*l] = xus[1+3*tn]; xus[2+3*l] = xus[2+3*tn];
 540                         xus[3*tn]  = x; xus[1+3*tn] = y; xus[2+3*tn] = z;
 541                         ijk        = work[l]; work[l] = work[tn]; work[tn] = ijk;
 542                         tn++; tl++;
 543                     }
 544                 }
 545                 *(ico_wk+ijk) = tl;
 546             } /* cycle k */
 547         }     /* cycle j */
 548     }         /* cycle i */
 549     sfree(ico_wk);
 550     sfree(work);
 551
 552     return xus;
 553 }
 554
 555 static void
 556 nsc_dclm_pbc(const rvec *coords, const ArrayRef<const real> &radius, int nat,
 557              const real *xus, int n_dot, int mode,
 558              real *value_of_area, real **at_area,
 559              real *value_of_vol,
 560              real **lidots, int *nu_dots,
 561              int index[], AnalysisNeighborhood *nb,
 562              const t_pbc *pbc)
 563 {
 564     const real dotarea = FOURPI/static_cast<real>(n_dot);
 565
 566     if (debug)
 567     {
 568         fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
 569     }
 570
 571     /* start with neighbour list */
 572     /* calculate neighbour list with the box algorithm */
 573     if (nat == 0)
 574     {
 575         return;
 576     }
 577     real        area = 0.0, vol = 0.0;
 578     real       *dots = nullptr, *atom_area = nullptr;
 579     int         lfnr = 0, maxdots = 0;
 580     if (mode & FLAG_VOLUME)
 581     {
 582         vol = 0.;
 583     }
 584     if (mode & FLAG_DOTS)
 585     {
 586         maxdots = (3*n_dot*nat)/10;
 587         snew(dots, maxdots);
 588         lfnr = 0;
 589     }
 590     if (mode & FLAG_ATOM_AREA)
 591     {
 592         snew(atom_area, nat);
 593     }
 594
 595     // Compute the center of the molecule for volume calculation.
 596     // In principle, the center should not influence the results, but that is
 597     // only true at the limit of infinite dot density, so this makes the
 598     // results translation-invariant.
 599     // With PBC, if the molecule is broken across the boundary, the computation
 600     // is broken in other ways as well, so it does not need to be considered
 601     // here.
 602     real xs = 0.0, ys = 0.0, zs = 0.0;
 603     for (int i = 0; i < nat; ++i)
 604     {
 605         const int iat = index[i];
 606         xs += coords[iat][XX];
 607         ys += coords[iat][YY];
 608         zs += coords[iat][ZZ];
 609     }
 610     xs /= nat;
 611     ys /= nat;
 612     zs /= nat;
 613
 614     AnalysisNeighborhoodPositions pos(coords, radius.size());
 615     pos.indexed(constArrayRefFromArray(index, nat));
 616     AnalysisNeighborhoodSearch    nbsearch(nb->initSearch(pbc, pos));
 617
 618     std::vector<int>              wkdot(n_dot);
 619
 620     for (int i = 0; i < nat; ++i)
 621     {
 622         const int                      iat  = index[i];
 623         const real                     ai   = radius[iat];
 624         const real                     aisq = ai*ai;
 625         AnalysisNeighborhoodPairSearch pairSearch(
 626                 nbsearch.startPairSearch(coords[iat]));
 627         AnalysisNeighborhoodPair       pair;
 628         std::fill(wkdot.begin(), wkdot.end(), 1);
 629         int currDotCount = n_dot;
 630         while (currDotCount > 0 && pairSearch.findNextPair(&pair))
 631         {
 632             const int  jat = index[pair.refIndex()];
 633             const real aj  = radius[jat];
 634             const real d2  = pair.distance2();
 635             if (iat == jat || d2 > gmx::square(ai+aj))
 636             {
 637                 continue;
 638             }
 639             const rvec &dx     = pair.dx();
 640             const real  refdot = (d2 + aisq - aj*aj)/(2*ai);
 641             // TODO: Consider whether micro-optimizations from the old
 642             // implementation would be useful, compared to the complexity that
 643             // they bring: instead of this direct loop, the neighbors were
 644             // stored into a temporary array, the loop order was
 645             // reversed (first over dots, then over neighbors), and for each
 646             // dot, it was first checked whether the same neighbor that
 647             // resulted in marking the previous dot covered would also cover
 648             // this dot. This presumably plays together with sorting of the
 649             // surface dots (done in make_unsp) to avoid some of the looping.
 650             // Alternatively, we could keep a skip list here to avoid
 651             // repeatedly looping over dots that have already marked as
 652             // covered.
 653             for (int j = 0; j < n_dot; ++j)
 654             {
 655                 if (wkdot[j] && iprod(&xus[3*j], dx) > refdot)
 656                 {
 657                     --currDotCount;
 658                     wkdot[j] = 0;
 659                 }
 660             }
 661         }
 662
 663         const real a = aisq * dotarea * currDotCount;
 664         area = area + a;
 665         if (mode & FLAG_ATOM_AREA)
 666         {
 667             atom_area[i] = a;
 668         }
 669         const real xi = coords[iat][XX];
 670         const real yi = coords[iat][YY];
 671         const real zi = coords[iat][ZZ];
 672         if (mode & FLAG_DOTS)
 673         {
 674             for (int l = 0; l < n_dot; l++)
 675             {
 676                 if (wkdot[l])
 677                 {
 678                     lfnr++;
 679                     if (maxdots <= 3*lfnr+1)
 680                     {
 681                         maxdots = maxdots+n_dot*3;
 682                         srenew(dots, maxdots);
 683                     }
 684                     dots[3*lfnr-3] = ai*xus[3*l]+xi;
 685                     dots[3*lfnr-2] = ai*xus[1+3*l]+yi;
 686                     dots[3*lfnr-1] = ai*xus[2+3*l]+zi;
 687                 }
 688             }
 689         }
 690         if (mode & FLAG_VOLUME)
 691         {
 692             real dx = 0.0, dy = 0.0, dz = 0.0;
 693             for (int l = 0; l < n_dot; l++)
 694             {
 695                 if (wkdot[l])
 696                 {
 697                     dx = dx+xus[3*l];
 698                     dy = dy+xus[1+3*l];
 699                     dz = dz+xus[2+3*l];
 700                 }
 701             }
 702             vol = vol+aisq*(dx*(xi-xs)+dy*(yi-ys)+dz*(zi-zs) + ai*currDotCount);
 703         }
 704     }
 705
 706     if (mode & FLAG_VOLUME)
 707     {
 708         *value_of_vol = vol*FOURPI/(3.*n_dot);
 709     }
 710     if (mode & FLAG_DOTS)
 711     {
 712         *nu_dots = lfnr;
 713         *lidots  = dots;
 714     }
 715     if (mode & FLAG_ATOM_AREA)
 716     {
 717         *at_area = atom_area;
 718     }
 719     *value_of_area = area;
 720
 721     if (debug)
 722     {
 723         fprintf(debug, "area=%8.3f\n", area);
 724     }
 725 }
 726
 727 namespace gmx
 728 {
 729
 730 class SurfaceAreaCalculator::Impl
 731 {
 732     public:
 733         Impl() : flags_(0)
 734         {
 735         }
 736
 737         std::vector<real>             unitSphereDots_;
 738         ArrayRef<const real>          radius_;
 739         int                           flags_;
 740         mutable AnalysisNeighborhood  nb_;
 741 };
 742
 743 SurfaceAreaCalculator::SurfaceAreaCalculator()
 744     : impl_(new Impl())
 745 {
 746 }
 747
 748 SurfaceAreaCalculator::~SurfaceAreaCalculator()
 749 {
 750 }
 751
 752 void SurfaceAreaCalculator::setDotCount(int dotCount)
 753 {
 754     impl_->unitSphereDots_ = make_unsp(dotCount, 4);
 755 }
 756
 757 void SurfaceAreaCalculator::setRadii(const ArrayRef<const real> &radius)
 758 {
 759     impl_->radius_ = radius;
 760     if (!radius.empty())
 761     {
 762         const real maxRadius = *std::max_element(radius.begin(), radius.end());
 763         impl_->nb_.setCutoff(2*maxRadius);
 764     }
 765 }
 766
 767 void SurfaceAreaCalculator::setCalculateVolume(bool bVolume)
 768 {
 769     if (bVolume)
 770     {
 771         impl_->flags_ |= FLAG_VOLUME;
 772     }
 773     else
 774     {
 775         impl_->flags_ &= ~FLAG_VOLUME;
 776     }
 777 }
 778
 779 void SurfaceAreaCalculator::setCalculateAtomArea(bool bAtomArea)
 780 {
 781     if (bAtomArea)
 782     {
 783         impl_->flags_ |= FLAG_ATOM_AREA;
 784     }
 785     else
 786     {
 787         impl_->flags_ &= ~FLAG_ATOM_AREA;
 788     }
 789 }
 790
 791 void SurfaceAreaCalculator::setCalculateSurfaceDots(bool bDots)
 792 {
 793     if (bDots)
 794     {
 795         impl_->flags_ |= FLAG_DOTS;
 796     }
 797     else
 798     {
 799         impl_->flags_ &= ~FLAG_DOTS;
 800     }
 801 }
 802
 803 void SurfaceAreaCalculator::calculate(
 804         const rvec *x, const t_pbc *pbc,
 805         int nat, int index[], int flags, real *area, real *volume,
 806         real **at_area, real **lidots, int *n_dots) const
 807 {
 808     flags |= impl_->flags_;
 809     *area  = 0;
 810     if (volume == nullptr)
 811     {
 812         flags &= ~FLAG_VOLUME;
 813     }
 814     else
 815     {
 816         *volume = 0;
 817     }
 818     if (at_area == nullptr)
 819     {
 820         flags &= ~FLAG_ATOM_AREA;
 821     }
 822     else
 823     {
 824         *at_area = nullptr;
 825     }
 826     if (lidots == nullptr)
 827     {
 828         flags &= ~FLAG_DOTS;
 829     }
 830     else
 831     {
 832         *lidots = nullptr;
 833     }
 834     if (n_dots == nullptr)
 835     {
 836         flags &= ~FLAG_DOTS;
 837     }
 838     else
 839     {
 840         *n_dots = 0;
 841     }
 842     nsc_dclm_pbc(x, impl_->radius_, nat,
 843                  &impl_->unitSphereDots_[0], impl_->unitSphereDots_.size()/3,
 844                  flags, area, at_area, volume, lidots, n_dots, index,
 845                  &impl_->nb_, pbc);
 846 }
 847
 848 } // namespace gmx