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