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