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