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],
 197                            xus[1 + 3 * i],
 198                            xus[2 + 3 * i],
 199                            xus[3 * j],
 200                            xus[1 + 3 * j],
 201                            xus[2 + 3 * j],
 202                            tl,
 203                            tess,
 204                            &xus[3 * tn],
 205                            &xus[1 + 3 * tn],
 206                            &xus[2 + 3 * tn]);
 207                     tn++;
 208                 }
 209             }
 210         }
 211         /* calculate tessalation of icosaeder faces */
 212         for (i = 0; i < 10; i++)
 213         {
 214             for (j = i + 1; j < 11; j++)
 215             {
 216                 x = xus[3 * i] - xus[3 * j];
 217                 y = xus[1 + 3 * i] - xus[1 + 3 * j];
 218                 z = xus[2 + 3 * i] - xus[2 + 3 * j];
 219                 d = x * x + y * y + z * z;
 220                 if (std::fabs(a - d) > DP_TOL)
 221                 {
 222                     continue;
 223                 }
 224
 225                 for (k = j + 1; k < 12; k++)
 226                 {
 227                     x = xus[3 * i] - xus[3 * k];
 228                     y = xus[1 + 3 * i] - xus[1 + 3 * k];
 229                     z = xus[2 + 3 * i] - xus[2 + 3 * k];
 230                     d = x * x + y * y + z * z;
 231                     if (std::fabs(a - d) > DP_TOL)
 232                     {
 233                         continue;
 234                     }
 235                     x = xus[3 * j] - xus[3 * k];
 236                     y = xus[1 + 3 * j] - xus[1 + 3 * k];
 237                     z = xus[2 + 3 * j] - xus[2 + 3 * k];
 238                     d = x * x + y * y + z * z;
 239                     if (std::fabs(a - d) > DP_TOL)
 240                     {
 241                         continue;
 242                     }
 243                     for (tl = 1; tl < tess - 1; tl++)
 244                     {
 245                         divarc(xus[3 * j],
 246                                xus[1 + 3 * j],
 247                                xus[2 + 3 * j],
 248                                xus[3 * i],
 249                                xus[1 + 3 * i],
 250                                xus[2 + 3 * i],
 251                                tl,
 252                                tess,
 253                                &xji,
 254                                &yji,
 255                                &zji);
 256                         divarc(xus[3 * k],
 257                                xus[1 + 3 * k],
 258                                xus[2 + 3 * k],
 259                                xus[3 * i],
 260                                xus[1 + 3 * i],
 261                                xus[2 + 3 * i],
 262                                tl,
 263                                tess,
 264                                &xki,
 265                                &yki,
 266                                &zki);
 267
 268                         for (tl2 = 1; tl2 < tess - tl; tl2++)
 269                         {
 270                             divarc(xus[3 * i],
 271                                    xus[1 + 3 * i],
 272                                    xus[2 + 3 * i],
 273                                    xus[3 * j],
 274                                    xus[1 + 3 * j],
 275                                    xus[2 + 3 * j],
 276                                    tl2,
 277                                    tess,
 278                                    &xij,
 279                                    &yij,
 280                                    &zij);
 281                             divarc(xus[3 * k],
 282                                    xus[1 + 3 * k],
 283                                    xus[2 + 3 * k],
 284                                    xus[3 * j],
 285                                    xus[1 + 3 * j],
 286                                    xus[2 + 3 * j],
 287                                    tl2,
 288                                    tess,
 289                                    &xkj,
 290                                    &ykj,
 291                                    &zkj);
 292                             divarc(xus[3 * i],
 293                                    xus[1 + 3 * i],
 294                                    xus[2 + 3 * i],
 295                                    xus[3 * k],
 296                                    xus[1 + 3 * k],
 297                                    xus[2 + 3 * k],
 298                                    tess - tl - tl2,
 299                                    tess,
 300                                    &xik,
 301                                    &yik,
 302                                    &zik);
 303                             divarc(xus[3 * j],
 304                                    xus[1 + 3 * j],
 305                                    xus[2 + 3 * j],
 306                                    xus[3 * k],
 307                                    xus[1 + 3 * k],
 308                                    xus[2 + 3 * k],
 309                                    tess - tl - tl2,
 310                                    tess,
 311                                    &xjk,
 312                                    &yjk,
 313                                    &zjk);
 314                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess - tl, &x, &y, &z);
 315                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess - tl2, &x2, &y2, &z2);
 316                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl + tl2, &x3, &y3, &z3);
 317                             GMX_ASSERT(tn < ndot, "Inconsistent precomputed surface dot count");
 318                             x               = x + x2 + x3;
 319                             y               = y + y2 + y3;
 320                             z               = z + z2 + z3;
 321                             d               = std::sqrt(x * x + y * y + z * z);
 322                             xus[3 * tn]     = x / d;
 323                             xus[1 + 3 * tn] = y / d;
 324                             xus[2 + 3 * tn] = z / d;
 325                             tn++;
 326                         } /* cycle tl2 */
 327                     }     /* cycle tl */
 328                 }         /* cycle k */
 329             }             /* cycle j */
 330         }                 /* cycle i */
 331         GMX_ASSERT(tn == ndot, "Inconsistent precomputed surface dot count");
 332     } /* end of if (tess > 1) */
 333
 334     return xus;
 335 } /* end of routine ico_dot_arc */
 336
 337 /* densit...required dots per unit sphere */
 338 static std::vector<real> ico_dot_dod(int densit)
 339 {
 340     /* dot distribution on a unit sphere based on an icosaeder *
 341      * great circle average refining of icosahedral face       */
 342
 343     int  i, j, k, tl, tl2, tn, tess, j1, j2;
 344     real a, d, x, y, z, x2, y2, z2, x3, y3, z3, ai_d, adod;
 345     real xij, yij, zij, xji, yji, zji, xik, yik, zik, xki, yki, zki, xjk, yjk, zjk, xkj, ykj, zkj;
 346
 347     /* calculate tesselation level */
 348     a              = std::sqrt(((static_cast<real>(densit)) - 2.) / 30.);
 349     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     tn = 12;
 357     /* square of the edge of an icosaeder */
 358     a = rh * rh * 2. * (1. - cos(TORAD(72.)));
 359     /* dodecaeder vertices */
 360     for (i = 0; i < 10; i++)
 361     {
 362         for (j = i + 1; j < 11; j++)
 363         {
 364             x = xus[3 * i] - xus[3 * j];
 365             y = xus[1 + 3 * i] - xus[1 + 3 * j];
 366             z = xus[2 + 3 * i] - xus[2 + 3 * j];
 367             d = x * x + y * y + z * z;
 368             if (std::fabs(a - d) > DP_TOL)
 369             {
 370                 continue;
 371             }
 372             for (k = j + 1; k < 12; k++)
 373             {
 374                 x = xus[3 * i] - xus[3 * k];
 375                 y = xus[1 + 3 * i] - xus[1 + 3 * k];
 376                 z = xus[2 + 3 * i] - xus[2 + 3 * k];
 377                 d = x * x + y * y + z * z;
 378                 if (std::fabs(a - d) > DP_TOL)
 379                 {
 380                     continue;
 381                 }
 382                 x = xus[3 * j] - xus[3 * k];
 383                 y = xus[1 + 3 * j] - xus[1 + 3 * k];
 384                 z = xus[2 + 3 * j] - xus[2 + 3 * k];
 385                 d = x * x + y * y + z * z;
 386                 if (std::fabs(a - d) > DP_TOL)
 387                 {
 388                     continue;
 389                 }
 390                 x               = xus[3 * i] + xus[3 * j] + xus[3 * k];
 391                 y               = xus[1 + 3 * i] + xus[1 + 3 * j] + xus[1 + 3 * k];
 392                 z               = xus[2 + 3 * i] + xus[2 + 3 * j] + xus[2 + 3 * k];
 393                 d               = std::sqrt(x * x + y * y + z * z);
 394                 xus[3 * tn]     = x / d;
 395                 xus[1 + 3 * tn] = y / d;
 396                 xus[2 + 3 * tn] = z / d;
 397                 tn++;
 398             }
 399         }
 400     }
 401
 402     if (tess > 1)
 403     {
 404         tn = 32;
 405         /* square of the edge of an dodecaeder */
 406         adod = 4. * (cos(TORAD(108.)) - cos(TORAD(120.))) / (1. - cos(TORAD(120.)));
 407         /* square of the distance of two adjacent vertices of ico- and dodecaeder */
 408         ai_d = 2. * (1. - std::sqrt(1. - a / 3.));
 409
 410         /* calculate tessalation of mixed edges */
 411         for (i = 0; i < 31; i++)
 412         {
 413             j1 = 12;
 414             j2 = 32;
 415             a  = ai_d;
 416             if (i >= 12)
 417             {
 418                 j1 = i + 1;
 419                 a  = adod;
 420             }
 421             for (j = j1; j < j2; j++)
 422             {
 423                 x = xus[3 * i] - xus[3 * j];
 424                 y = xus[1 + 3 * i] - xus[1 + 3 * j];
 425                 z = xus[2 + 3 * i] - xus[2 + 3 * j];
 426                 d = x * x + y * y + z * z;
 427                 if (std::fabs(a - d) > DP_TOL)
 428                 {
 429                     continue;
 430                 }
 431                 for (tl = 1; tl < tess; tl++)
 432                 {
 433                     GMX_ASSERT(tn < ndot, "Inconsistent precomputed surface dot count");
 434                     divarc(xus[3 * i],
 435                            xus[1 + 3 * i],
 436                            xus[2 + 3 * i],
 437                            xus[3 * j],
 438                            xus[1 + 3 * j],
 439                            xus[2 + 3 * j],
 440                            tl,
 441                            tess,
 442                            &xus[3 * tn],
 443                            &xus[1 + 3 * tn],
 444                            &xus[2 + 3 * tn]);
 445                     tn++;
 446                 }
 447             }
 448         }
 449         /* calculate tessalation of pentakisdodecahedron faces */
 450         for (i = 0; i < 12; i++)
 451         {
 452             for (j = 12; j < 31; j++)
 453             {
 454                 x = xus[3 * i] - xus[3 * j];
 455                 y = xus[1 + 3 * i] - xus[1 + 3 * j];
 456                 z = xus[2 + 3 * i] - xus[2 + 3 * j];
 457                 d = x * x + y * y + z * z;
 458                 if (std::fabs(ai_d - d) > DP_TOL)
 459                 {
 460                     continue;
 461                 }
 462
 463                 for (k = j + 1; k < 32; k++)
 464                 {
 465                     x = xus[3 * i] - xus[3 * k];
 466                     y = xus[1 + 3 * i] - xus[1 + 3 * k];
 467                     z = xus[2 + 3 * i] - xus[2 + 3 * k];
 468                     d = x * x + y * y + z * z;
 469                     if (std::fabs(ai_d - d) > DP_TOL)
 470                     {
 471                         continue;
 472                     }
 473                     x = xus[3 * j] - xus[3 * k];
 474                     y = xus[1 + 3 * j] - xus[1 + 3 * k];
 475                     z = xus[2 + 3 * j] - xus[2 + 3 * k];
 476                     d = x * x + y * y + z * z;
 477                     if (std::fabs(adod - d) > DP_TOL)
 478                     {
 479                         continue;
 480                     }
 481                     for (tl = 1; tl < tess - 1; tl++)
 482                     {
 483                         divarc(xus[3 * j],
 484                                xus[1 + 3 * j],
 485                                xus[2 + 3 * j],
 486                                xus[3 * i],
 487                                xus[1 + 3 * i],
 488                                xus[2 + 3 * i],
 489                                tl,
 490                                tess,
 491                                &xji,
 492                                &yji,
 493                                &zji);
 494                         divarc(xus[3 * k],
 495                                xus[1 + 3 * k],
 496                                xus[2 + 3 * k],
 497                                xus[3 * i],
 498                                xus[1 + 3 * i],
 499                                xus[2 + 3 * i],
 500                                tl,
 501                                tess,
 502                                &xki,
 503                                &yki,
 504                                &zki);
 505
 506                         for (tl2 = 1; tl2 < tess - tl; tl2++)
 507                         {
 508                             divarc(xus[3 * i],
 509                                    xus[1 + 3 * i],
 510                                    xus[2 + 3 * i],
 511                                    xus[3 * j],
 512                                    xus[1 + 3 * j],
 513                                    xus[2 + 3 * j],
 514                                    tl2,
 515                                    tess,
 516                                    &xij,
 517                                    &yij,
 518                                    &zij);
 519                             divarc(xus[3 * k],
 520                                    xus[1 + 3 * k],
 521                                    xus[2 + 3 * k],
 522                                    xus[3 * j],
 523                                    xus[1 + 3 * j],
 524                                    xus[2 + 3 * j],
 525                                    tl2,
 526                                    tess,
 527                                    &xkj,
 528                                    &ykj,
 529                                    &zkj);
 530                             divarc(xus[3 * i],
 531                                    xus[1 + 3 * i],
 532                                    xus[2 + 3 * i],
 533                                    xus[3 * k],
 534                                    xus[1 + 3 * k],
 535                                    xus[2 + 3 * k],
 536                                    tess - tl - tl2,
 537                                    tess,
 538                                    &xik,
 539                                    &yik,
 540                                    &zik);
 541                             divarc(xus[3 * j],
 542                                    xus[1 + 3 * j],
 543                                    xus[2 + 3 * j],
 544                                    xus[3 * k],
 545                                    xus[1 + 3 * k],
 546                                    xus[2 + 3 * k],
 547                                    tess - tl - tl2,
 548                                    tess,
 549                                    &xjk,
 550                                    &yjk,
 551                                    &zjk);
 552                             divarc(xki, yki, zki, xji, yji, zji, tl2, tess - tl, &x, &y, &z);
 553                             divarc(xkj, ykj, zkj, xij, yij, zij, tl, tess - tl2, &x2, &y2, &z2);
 554                             divarc(xjk, yjk, zjk, xik, yik, zik, tl, tl + tl2, &x3, &y3, &z3);
 555                             GMX_ASSERT(tn < ndot, "Inconsistent precomputed surface dot count");
 556                             x               = x + x2 + x3;
 557                             y               = y + y2 + y3;
 558                             z               = z + z2 + z3;
 559                             d               = std::sqrt(x * x + y * y + z * z);
 560                             xus[3 * tn]     = x / d;
 561                             xus[1 + 3 * tn] = y / d;
 562                             xus[2 + 3 * tn] = z / d;
 563                             tn++;
 564                         } /* cycle tl2 */
 565                     }     /* cycle tl */
 566                 }         /* cycle k */
 567             }             /* cycle j */
 568         }                 /* cycle i */
 569         GMX_ASSERT(tn == ndot, "Inconsistent precomputed surface dot count");
 570     } /* end of if (tess > 1) */
 571
 572     return xus;
 573 } /* end of routine ico_dot_dod */
 574
 575 static int unsp_type(int densit)
 576 {
 577     int i1, i2;
 578     i1 = 1;
 579     while (10 * i1 * i1 + 2 < densit)
 580     {
 581         i1++;
 582     }
 583     i2 = 1;
 584     while (30 * i2 * i2 + 2 < densit)
 585     {
 586         i2++;
 587     }
 588     if (10 * i1 * i1 - 2 < 30 * i2 * i2 - 2)
 589     {
 590         return UNSP_ICO_ARC;
 591     }
 592     else
 593     {
 594         return UNSP_ICO_DOD;
 595     }
 596 }
 597
 598 static std::vector<real> make_unsp(int densit, int cubus)
 599 {
 600     int *ico_wk, *ico_pt;
 601     int  ico_cube, ico_cube_cb, i, j, k, l, ijk, tn, tl, tl2;
 602     int* work;
 603     real x, y, z;
 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         i = 1;
 630         while (i * i * i * 2 < ndot)
 631         {
 632             i++;
 633         }
 634         ico_cube = std::max(i - 1, 0);
 635     }
 636     ico_cube_cb         = ico_cube * ico_cube * ico_cube;
 637     const real del_cube = 2. / (static_cast<real>(ico_cube));
 638     snew(work, ndot);
 639     for (l = 0; l < ndot; l++)
 640     {
 641         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         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         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         ijk     = i + j * ico_cube + k * ico_cube * ico_cube;
 657         work[l] = ijk;
 658     }
 659
 660     snew(ico_wk, 2 * ico_cube_cb + 1);
 661
 662     ico_pt = ico_wk + ico_cube_cb;
 663     for (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     tn = 0;
 670     for (i = 0; i < ico_cube; i++)
 671     {
 672         for (j = 0; j < ico_cube; j++)
 673         {
 674             for (k = 0; k < ico_cube; k++)
 675             {
 676                 tl              = 0;
 677                 tl2             = tn;
 678                 ijk             = i + ico_cube * j + ico_cube * ico_cube * k;
 679                 *(ico_pt + ijk) = tn;
 680                 for (l = tl2; l < ndot; l++)
 681                 {
 682                     if (ijk == work[l])
 683                     {
 684                         x               = xus[3 * l];
 685                         y               = xus[1 + 3 * l];
 686                         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 + ijk) = tl;
 701             } /* cycle k */
 702         }     /* cycle j */
 703     }         /* cycle i */
 704     sfree(ico_wk);
 705     sfree(work);
 706
 707     return xus;
 708 }
 709
 710 static void nsc_dclm_pbc(const rvec*                 coords,
 711                          const ArrayRef<const real>& radius,
 712                          int                         nat,
 713                          const real*                 xus,
 714                          int                         n_dot,
 715                          int                         mode,
 716                          real*                       value_of_area,
 717                          real**                      at_area,
 718                          real*                       value_of_vol,
 719                          real**                      lidots,
 720                          int*                        nu_dots,
 721                          int                         index[],
 722                          AnalysisNeighborhood*       nb,
 723                          const t_pbc*                pbc)
 724 {
 725     const real dotarea = FOURPI / static_cast<real>(n_dot);
 726
 727     if (debug)
 728     {
 729         fprintf(debug, "nsc_dclm: n_dot=%5d %9.3f\n", n_dot, dotarea);
 730     }
 731
 732     /* start with neighbour list */
 733     /* calculate neighbour list with the box algorithm */
 734     if (nat == 0)
 735     {
 736         return;
 737     }
 738     real  area = 0.0, vol = 0.0;
 739     real *dots = nullptr, *atom_area = nullptr;
 740     int   lfnr = 0, maxdots = 0;
 741     if (mode & FLAG_VOLUME)
 742     {
 743         vol = 0.;
 744     }
 745     if (mode & FLAG_DOTS)
 746     {
 747         maxdots = (3 * n_dot * nat) / 10;
 748         snew(dots, maxdots);
 749         lfnr = 0;
 750     }
 751     if (mode & FLAG_ATOM_AREA)
 752     {
 753         snew(atom_area, nat);
 754     }
 755
 756     // Compute the center of the molecule for volume calculation.
 757     // In principle, the center should not influence the results, but that is
 758     // only true at the limit of infinite dot density, so this makes the
 759     // results translation-invariant.
 760     // With PBC, if the molecule is broken across the boundary, the computation
 761     // is broken in other ways as well, so it does not need to be considered
 762     // here.
 763     real xs = 0.0, ys = 0.0, zs = 0.0;
 764     for (int i = 0; i < nat; ++i)
 765     {
 766         const int iat = index[i];
 767         xs += coords[iat][XX];
 768         ys += coords[iat][YY];
 769         zs += coords[iat][ZZ];
 770     }
 771     xs /= nat;
 772     ys /= nat;
 773     zs /= nat;
 774
 775     AnalysisNeighborhoodPositions pos(coords, radius.size());
 776     pos.indexed(constArrayRefFromArray(index, nat));
 777     AnalysisNeighborhoodSearch nbsearch(nb->initSearch(pbc, pos));
 778
 779     std::vector<int> wkdot(n_dot);
 780
 781     for (int i = 0; i < nat; ++i)
 782     {
 783         const int                      iat  = index[i];
 784         const real                     ai   = radius[iat];
 785         const real                     aisq = ai * ai;
 786         AnalysisNeighborhoodPairSearch pairSearch(nbsearch.startPairSearch(coords[iat]));
 787         AnalysisNeighborhoodPair       pair;
 788         std::fill(wkdot.begin(), wkdot.end(), 1);
 789         int currDotCount = n_dot;
 790         while (currDotCount > 0 && pairSearch.findNextPair(&pair))
 791         {
 792             const int  jat = index[pair.refIndex()];
 793             const real aj  = radius[jat];
 794             const real d2  = pair.distance2();
 795             if (iat == jat || d2 > gmx::square(ai + aj))
 796             {
 797                 continue;
 798             }
 799             const rvec& dx     = pair.dx();
 800             const real  refdot = (d2 + aisq - aj * aj) / (2 * ai);
 801             // TODO: Consider whether micro-optimizations from the old
 802             // implementation would be useful, compared to the complexity that
 803             // they bring: instead of this direct loop, the neighbors were
 804             // stored into a temporary array, the loop order was
 805             // reversed (first over dots, then over neighbors), and for each
 806             // dot, it was first checked whether the same neighbor that
 807             // resulted in marking the previous dot covered would also cover
 808             // this dot. This presumably plays together with sorting of the
 809             // surface dots (done in make_unsp) to avoid some of the looping.
 810             // Alternatively, we could keep a skip list here to avoid
 811             // repeatedly looping over dots that have already marked as
 812             // covered.
 813             for (int j = 0; j < n_dot; ++j)
 814             {
 815                 if (wkdot[j] && iprod(&xus[3 * j], dx) > refdot)
 816                 {
 817                     --currDotCount;
 818                     wkdot[j] = 0;
 819                 }
 820             }
 821         }
 822
 823         const real a = aisq * dotarea * currDotCount;
 824         area         = area + a;
 825         if (mode & FLAG_ATOM_AREA)
 826         {
 827             atom_area[i] = a;
 828         }
 829         const real xi = coords[iat][XX];
 830         const real yi = coords[iat][YY];
 831         const real zi = coords[iat][ZZ];
 832         if (mode & FLAG_DOTS)
 833         {
 834             for (int l = 0; l < n_dot; l++)
 835             {
 836                 if (wkdot[l])
 837                 {
 838                     lfnr++;
 839                     if (maxdots <= 3 * lfnr + 1)
 840                     {
 841                         maxdots = maxdots + n_dot * 3;
 842                         srenew(dots, maxdots);
 843                     }
 844                     dots[3 * lfnr - 3] = ai * xus[3 * l] + xi;
 845                     dots[3 * lfnr - 2] = ai * xus[1 + 3 * l] + yi;
 846                     dots[3 * lfnr - 1] = ai * xus[2 + 3 * l] + zi;
 847                 }
 848             }
 849         }
 850         if (mode & FLAG_VOLUME)
 851         {
 852             real dx = 0.0, dy = 0.0, dz = 0.0;
 853             for (int l = 0; l < n_dot; l++)
 854             {
 855                 if (wkdot[l])
 856                 {
 857                     dx = dx + xus[3 * l];
 858                     dy = dy + xus[1 + 3 * l];
 859                     dz = dz + xus[2 + 3 * l];
 860                 }
 861             }
 862             vol = vol + aisq * (dx * (xi - xs) + dy * (yi - ys) + dz * (zi - zs) + ai * currDotCount);
 863         }
 864     }
 865
 866     if (mode & FLAG_VOLUME)
 867     {
 868         *value_of_vol = vol * FOURPI / (3. * n_dot);
 869     }
 870     if (mode & FLAG_DOTS)
 871     {
 872         GMX_RELEASE_ASSERT(nu_dots != nullptr, "Must have valid nu_dots pointer");
 873         *nu_dots = lfnr;
 874         GMX_RELEASE_ASSERT(lidots != nullptr, "Must have valid lidots pointer");
 875         *lidots = dots;
 876     }
 877     if (mode & FLAG_ATOM_AREA)
 878     {
 879         GMX_RELEASE_ASSERT(at_area != nullptr, "Must have valid at_area pointer");
 880         *at_area = atom_area;
 881     }
 882     *value_of_area = area;
 883
 884     if (debug)
 885     {
 886         fprintf(debug, "area=%8.3f\n", area);
 887     }
 888 }
 889
 890 namespace gmx
 891 {
 892
 893 class SurfaceAreaCalculator::Impl
 894 {
 895 public:
 896     Impl() : flags_(0) {}
 897
 898     std::vector<real>            unitSphereDots_;
 899     ArrayRef<const real>         radius_;
 900     int                          flags_;
 901     mutable AnalysisNeighborhood nb_;
 902 };
 903
 904 SurfaceAreaCalculator::SurfaceAreaCalculator() : impl_(new Impl()) {}
 905
 906 SurfaceAreaCalculator::~SurfaceAreaCalculator() {}
 907
 908 void SurfaceAreaCalculator::setDotCount(int dotCount)
 909 {
 910     impl_->unitSphereDots_ = make_unsp(dotCount, 4);
 911 }
 912
 913 void SurfaceAreaCalculator::setRadii(const ArrayRef<const real>& radius)
 914 {
 915     impl_->radius_ = radius;
 916     if (!radius.empty())
 917     {
 918         const real maxRadius = *std::max_element(radius.begin(), radius.end());
 919         impl_->nb_.setCutoff(2 * maxRadius);
 920     }
 921 }
 922
 923 void SurfaceAreaCalculator::setCalculateVolume(bool bVolume)
 924 {
 925     if (bVolume)
 926     {
 927         impl_->flags_ |= FLAG_VOLUME;
 928     }
 929     else
 930     {
 931         impl_->flags_ &= ~FLAG_VOLUME;
 932     }
 933 }
 934
 935 void SurfaceAreaCalculator::setCalculateAtomArea(bool bAtomArea)
 936 {
 937     if (bAtomArea)
 938     {
 939         impl_->flags_ |= FLAG_ATOM_AREA;
 940     }
 941     else
 942     {
 943         impl_->flags_ &= ~FLAG_ATOM_AREA;
 944     }
 945 }
 946
 947 void SurfaceAreaCalculator::setCalculateSurfaceDots(bool bDots)
 948 {
 949     if (bDots)
 950     {
 951         impl_->flags_ |= FLAG_DOTS;
 952     }
 953     else
 954     {
 955         impl_->flags_ &= ~FLAG_DOTS;
 956     }
 957 }
 958
 959 void SurfaceAreaCalculator::calculate(const rvec*  x,
 960                                       const t_pbc* pbc,
 961                                       int          nat,
 962                                       int          index[],
 963                                       int          flags,
 964                                       real*        area,
 965                                       real*        volume,
 966                                       real**       at_area,
 967                                       real**       lidots,
 968                                       int*         n_dots) const
 969 {
 970     flags |= impl_->flags_;
 971     *area = 0;
 972     if (volume == nullptr)
 973     {
 974         flags &= ~FLAG_VOLUME;
 975     }
 976     else
 977     {
 978         *volume = 0;
 979     }
 980     if (at_area == nullptr)
 981     {
 982         flags &= ~FLAG_ATOM_AREA;
 983     }
 984     else
 985     {
 986         *at_area = nullptr;
 987     }
 988     if (lidots == nullptr)
 989     {
 990         flags &= ~FLAG_DOTS;
 991     }
 992     else
 993     {
 994         *lidots = nullptr;
 995     }
 996     if (n_dots == nullptr)
 997     {
 998         flags &= ~FLAG_DOTS;
 999     }
1000     else
1001     {
1002         *n_dots = 0;
1003     }
1004     nsc_dclm_pbc(x,
1005                  impl_->radius_,
1006                  nat,
1007                  &impl_->unitSphereDots_[0],
1008                  impl_->unitSphereDots_.size() / 3,
1009                  flags,
1010                  area,
1011                  at_area,
1012                  volume,
1013                  lidots,
1014                  n_dots,
1015                  index,
1016                  &impl_->nb_,
1017                  pbc);
1018 }
1019
1020 } // namespace gmx